Amostragem da distribuição bivariada com densidade conhecida usando MCMC

9

Tentei simular a partir de uma densidade bivariada usando algoritmos Metropolis em R e não tive sorte. A densidade pode ser expressa como p ( y | x ) p ( x ) , onde p ( x ) é a distribuição de Singh-Maddalap(x,y)p(y|x)p(x)p(x)

p(x)=aqxa1ba(1+(xb)a)1+q

com os parâmetros de , q , b , e p ( y | x ) é o log-normal, com log-média como uma fracção de x , e log-sd uma constante. Para testar se minha amostra é a que eu quero, observei a densidade marginal de x , que deve ser p ( x ) . Eu tentei algoritmos diferentes do Metropolis dos pacotes R MCMCpack, mcmc e dream. Eu descartei a queima, usei o desbaste, usei amostras com tamanho de até milhões, mas a densidade marginal resultante nunca foi a que eu forneci.aqbp(y|x)xxp(x)

Aqui está a edição final do meu código que usei:

logvrls <- function(x,el,sdlog,a,scl,q.arg) {
    if(x[2]>0) {
         dlnorm(x[1],meanlog=el*log(x[2]),sdlog=sdlog,log=TRUE)+
         dsinmad(x[2],a=a,scale=scl,q.arg=q.arg,log=TRUE)
    }
    else -Inf    
}

a <- 1.35
q <- 3.3
scale <- 10/gamma(1 + 1/a)/gamma(q - 1/a)*  gamma(q) 

Initvrls <- function(pars,nseq,meanlog,sdlog,a,scale,q) {
    cbind(rlnorm(nseq,meanlog,sdlog),rsinmad(nseq,a,scale,q))
}

library(dream)
aa <- dream(logvrls,
        func.type="logposterior.density",
        pars=list(c(0,Inf),c(0,Inf)),
        FUN.pars=list(el=0.2,sdlog=0.2,a=a,scl=scale,q.arg=q),
        INIT=Initvrls,
        INIT.pars=list(meanlog=1,sdlog=0.1,a=a,scale=scale,q=q),
        control=list(nseq=3,thin.t=10)
        )

Eu me conformei com o pacote dos sonhos, desde amostras até a convergência. Eu testei se tenho os resultados corretos de três maneiras. Usando estatística KS, comparando quantis e estimando os parâmetros da distribuição Singh-Maddala com máxima probabilidade da amostra resultante:

ks.test(as.numeric(aa$Seq[[2]][,2]),psinmad,a=a,scale=scale,q.arg=q)

lsinmad <- function(x,sample)
    sum(dsinmad(sample,a=x[1],scale=x[2],q.arg=x[3],log=TRUE))
 optim(c(2,20,2),lsinmad,method="BFGS",sample=aa$Seq[[1]][,2])

 qq <- eq(0.025,.975,by=0.025)   
 tst <- cbind(qq,
              sapply(aa$Seq,function(l)round(quantile(l[,2],qq),3)),
              round(qsinmad(qq,a,scale,q),3))
 colnames(tst) <- c("Quantile","S1","S2","S3","True")

 library(ggplot2)
 qplot(x=Quantile,y=value,
       data=melt(data.frame(tst),id=1), 
       colour=variable,group=variable,geom="line")

Quando observo os resultados dessas comparações, a estatística KS quase sempre rejeita a hipótese nula de que a amostra é da distribuição de Singh-Maddala com os parâmetros fornecidos. Às vezes, os parâmetros estimados com probabilidade máxima aproximam-se dos valores reais, mas geralmente estão muito longe da zona de conforto, para aceitar que o procedimento de amostragem foi bem-sucedido. Idem para os quantis, quantis empíricos não estão muito longe, mas muito longe.

Minha pergunta é o que estou fazendo de errado? Minhas próprias hipóteses:

  1. O MCMC não é apropriado para este tipo de amostragem
  2. O MCMC não pode convergir devido a razões teóricas (a função de distribuição não satisfaz as propriedades necessárias, sejam elas quais forem)
  3. Eu não uso o algoritmo Metropolis corretamente
  4. Meus testes de distribuição não estão corretos, pois não tenho amostra independente.
mpiktas
fonte
No link de distribuição Singh-Maddala , o pdf possui dois parâmetros - {c, k}, mas a função R dsinmadusa três parâmetros ou estou faltando alguma coisa.
Csgillespie
Desculpe, o link da wikipedia cita a fórmula errada, parecia ok à primeira vista, quando eu estava compondo a pergunta. Não encontrei um link pronto, apenas coloquei a fórmula na pergunta.
mpiktas

Respostas:

3

Acho que a ordem está correta, mas os rótulos atribuídos a p (x) ep (y | x) estavam errados. O problema original indica que p (y | x) é log-normal ep (x) é Singh-Maddala. Então é

  1. Gere um X de um Singh-Maddala e

  2. gerar um Y a partir de um log-normal com uma média que é uma fração do X gerado.

Jan Galkowski
fonte
3

Na verdade, você não deve fazer o MCMC, pois seu problema é muito mais simples. Experimente este algoritmo:

Etapa 1: gerar um X a partir do log normal

Etapa 2: mantendo esse X fixo, gere um Y a partir do Singh Maddala.

Voilà! Amostra pronta !!!

Mohit
fonte
Presumo que você quis dizer os passos invertidos. Mas se isso é tão simples, por que precisamos da amostragem de Gibbs?
mpiktas
11
Não, eu quis dizer as etapas 1 e 2 na ordem em que escrevi. Afinal, a distribuição de y é especificada condicionalmente em X, portanto, você deve gerar um X antes de Y. Quanto à amostragem de Gibbs, essa é uma solução mais complicada para problemas mais complicados. O seu, como você o descreve, é bastante direto, IMHO.
Mohit
11
p(y|x)p(x|y)p(x)