Estou tentando declarar uma distribuição anterior para uma meta-análise bayesiana.
Eu tenho as seguintes informações sobre uma variável aleatória:
- Duas observações: 3.0, 3.6
- um cientista que estuda a variável me disse que , e que valores tão altos quanto 6 têm probabilidade diferente de zero.
Eu usei a seguinte abordagem para otimização (o modo de log-N = :
prior <- function(parms, x, alpha) {
a <- abs(plnorm(x[1], parms[1], parms[2]) - (alpha/2))
b <- abs(plnorm(x[2], parms[1], parms[2]) - (1-alpha/2))
mode <- exp(parms[1] - parms[2]^2)
c <- abs(mode-3.3)
return(a + b + c)
}
v = nlm(prior,c(log(3.3),0.14),alpha=0.05,x=c(2.5,7.5))
x <- seq(1,10,0.1)
plot(x, dlnorm(x, v$estimate[1], v$estimate[2]))
abline(v=c(2.5,7.5), lty=2) #95%CI
Na figura, você pode ver a distribuição que isso retorna, mas eu gostaria de encontrar algo mais parecido com as linhas vermelhas que desenhei.
Isso fornece a mesma distribuição modelada usando o lognormal, gama ou normal, e resulta em uma distribuição com e , ou seja:P ( X = 6 ) < 0,01
plnorm(c(5,6), v$estimate[1],v$estimate[2])
Alguém pode sugerir alternativas? Eu preferiria ficar com uma distribuição única do que com uma mistura.
Obrigado!
r
distributions
probability
bayesian
optimization
David LeBauer
fonte
fonte
Respostas:
Se, dada uma resposta ao meu comentário acima, você deseja limitar o intervalo da distribuição, por que não ajustar simplesmente uma distribuição Beta na qual você redimensiona o intervalo de unidade? Em outras palavras, se você souber que o parâmetro de interesse deve estar entre , por que não definir . Onde eu centralizei o intervalo em zero, dividido pela largura para que Y tenha um intervalo de 1 e, em seguida, adicionei volta para que o intervalo de Y seja . (Você pode pensar de qualquer maneira: diretamente de ou de[2,8] Y=X−56+12=X−26 12 [0,1] [2,8]→[0,1] [2,8]→[−12,12]→[0,1] , mas achei que o último poderia ser mais fácil a princípio).
Então, com dois pontos de dados, você poderia ajustar um beta posterior com um beta uniforme antes?
fonte
E a distribuição Kumaraswamy , que possui o seguinte pdf:
a > 0 b > 0 0 < x < 1
fonte
Como a distribuição log-normal possui dois parâmetros, você não pode ajustá-la satisfatoriamente a três restrições que não se ajustam naturalmente a ela. Com quantis extremos de 2,5 e 7,5, o modo é ~ 4 e não há muito que você possa fazer sobre isso. Como a escala dos erros para
a
eb
é muito menor que parac
, um deles será praticamente ignorado durante a otimização.Para um melhor ajuste, você pode escolher uma distribuição de três parâmetros, por exemplo, a distribuição gama generalizada (implementada no
VGAM
pacote) ou adicionar um parâmetro shift à distribuição lognormal (ou gama, ...).Como última observação, como a distribuição que você está procurando claramente não é simétrica, a média das duas observações fornecidas não é o valor correto para o modo. Eu maximizaria a soma das densidades em 3,0 e 3,6, mantendo os quantis extremos em 2,5 e 7,5 - isso é possível se você tiver três parâmetros.
fonte
Você também pode tentar a distribuição triangular. Para ajustar isso, você basicamente especifica um limite inferior (esse seria X = 2), um limite superior (esse seria X = 8) e um valor "mais provável". A página da wikepedia http://en.wikipedia.org/wiki/Triangular_distribution possui mais informações sobre esta distribuição. Se não houver muita fé no valor "mais provável" (como parece ser, antes da observação de qualquer dado), pode ser uma boa ideia colocar uma distribuição prévia não informativa e usar os dois dados pontos para estimar esse valor. Um bom é o prior de jeffrey, que para esse problema seria p (c) = 1 / (pi * sqrt ((c-2) * (c-8))), onde "c" é o "valor mais provável "(consistente com a notação da Wikipedia).
Dado isso anteriormente, é possível calcular a distribuição posterior de c analiticamente ou via simulação. A forma analítica da probabilidade não é particularmente agradável, portanto a simulação parece ser mais atraente. Este exemplo é particularmente adequado para amostragem por rejeição (consulte a página da wiki para obter uma descrição geral da amostragem por rejeição), porque a probabilidade máxima é de 1/3 ^ n, independentemente do valor de c, que fornece o "limite superior". Então, você gera um "candidato" a partir do prior de jeffrey (chame de c_i) e avalia a probabilidade desse candidato L (x1, .., xn | c_i) e divide pela probabilidade maximizada para dar (3 ^ n ) * L (x1, .., xn | c_i). Você gera uma variável aleatória U (0,1) e, se u for menor que (3 ^ n) * L (x1, .., xn | c_i), aceite c_i como um valor amostrado posterior, caso contrário, jogue fora c_i e comece de novo. Repita esse processo até ter amostras aceitas suficientes (100, 500, 1.000 ou mais, dependendo da precisão que você deseja). Em seguida, faça a média da amostra de qualquer função de c em que você esteja interessado (a probabilidade de uma nova observação é um candidato óbvio para sua aplicação).
Uma alternativa para aceitar-rejeitar é usar o valor da probabilidade como um peso (e não gerar u) e, em seguida, continuar com as médias ponderadas usando todos os candidatos, em vez de médias não ponderadas com os candidatos aceitos.
fonte