Amostragem inversa de CDF para uma distribuição mista

9

A versão curta fora do contexto

Seja uma variável aleatória com CDF y

F(){θ y = 0 θ+(1θ)×CDFlog-normal(;μ,σ) y > 0

Digamos que eu queira simular empates de usando o método CDF inverso. Isso é possível? Esta função não tem exatamente um inverso. Então, novamente, há amostragem de transformação inversa para distribuição de mistura de duas distribuições normais, o que sugere que há uma maneira conhecida de aplicar amostragem de transformação inversa aqui.y

Estou ciente do método de duas etapas, mas não sei como aplicá-lo à minha situação (veja abaixo).


A versão longa com fundo

Eu ajustei o seguinte modelo para uma resposta com valor vetorial, , usando o MCMC (especificamente, Stan):yi=(y1,,yK)i

θkilogit1(αkxi),μkiβkxiσk22F(){θ y = 0 θ+(1θ)×CDFlog-normal(;μ,σ) y > 0ukF(yk),zkΦ1(uk)zN(0,R)×kf(yk)(α,β,σ,R)priors

onde indexa observações, é uma matriz de correlação e é um vetor de preditores / regressores / características.N R xiNRx

Ou seja, meu modelo é um modelo de regressão em que a distribuição condicional da resposta é assumida como uma cópula gaussiana com marginais log-normais inflados a zero. Eu postei sobre esse modelo antes; Acontece que Song, Li e Yuan (2009, gated ) o desenvolveram e o chamam de vetor GLM, ou VGLM. A seguir, a especificação deles é o mais aproximada possível: MeuF K G m z q R Γ

f(y;μ,φ,Γ)=c{G1(y1),,Gm(ym)|Γ}i=1mg(yi;μi,φi)c(u|Γ)=|Γ|1/2exp(12qT(ImΓ1)q)q=(q1,,qm)T,qi=Φ1(ui)
FKcorresponde ao , meu corresponde ao e meu corresponde ao ; os detalhes estão na página 62 (página 3 do arquivo PDF), mas são idênticos ao que escrevi aqui.GmzqRΓ

A peça inflada com zero segue aproximadamente as especificações de Liu e Chan (2010, sem porta ).

Agora, eu gostaria de simular dados dos parâmetros estimados, mas estou um pouco confuso sobre como fazê-lo. Primeiro, pensei em simular diretamente (no código R):y

for (i in 1:N) {
    for (k in 1:K) {
        Y_hat <- rbinom(1, 1, 1 - theta[i, k])
        if (Y_hat == 1)
            Y_hat <- rlnorm(1, mu[i, k], sigma[k])
    }
}

que não usa em tudo. Eu gostaria de tentar realmente usar a matriz de correlação que estimei.R

Minha próxima idéia foi pegar e convertê-los de volta para . Isso também parece coincidir com as respostas em Gerando amostras de Copula na amostragem R e bivariada para distribuição expressa no teorema da cópula de Sklar? . Mas que diabos é o meu aqui? A amostragem de transformação inversa para distribuição de mistura de duas distribuições normais faz parecer que isso é possível, mas não tenho idéia de como fazê-lo.e F - 1zyF1

shadowtalker
fonte
@ Xi'an é uma cópula gaussiana, para estimar a dependência entre os componentes . y
shadowtalker
11
O encadeamento que você menciona sobre a amostragem de misturas de normais aplica-se diretamente ao seu problema sem modificação essencial: em vez de usar os CDFs inversos dos normais, use os CDFs inversos de seus dois componentes. O CDF inverso do átomo em é a função constante, sempre igual a . 0y=00
whuber
@whuber Estou confuso sobre como usar os CDFs inversos dos dois componentes: o que eu desenho, o que eu desenho e o que eu conecto cada coisa?
shadowtalker
11
@ Xi'an explica bem que, em sua resposta à pergunta da mistura normal: você usa uma variável uniforme para selecionar o componente da mistura e depois extrai um valor desse componente (da maneira que desejar). No seu caso, é extremamente fácil desenhar um valor a partir do primeiro componente: é sempre ! Para desenhar um valor do segundo componente, use qualquer gerador de números aleatórios lognormal que você desejar. Em cada caso, você acaba com um número: não há "plug-in" para realizar; todo o objetivo da geração aleatória de números é obter esse número. 0
whuber
@whuber a nova resposta esclareceu para mim. Obrigado a ambos.
shadowtalker

Respostas:

5

A resposta para a versão longa com fundo:

Essa resposta para a versão longa aborda um pouco outra questão e, como parecemos ter dificuldades em formular o modelo e o problema, escolhi reformulá-lo aqui, espero que corretamente.

Para , o objetivo é simular vetores modo que, condicional a uma covariável , com . Portanto, se alguém quiser simular dados desse modelo, poderá proceder da seguinte maneira:1iIyi=(y1i,,yKi)xi

yki={0 with probability logit1(αkxi)log(σkzki+βkxi) with probability 1logit1(αkxi)
zi=(z1i,,zKi)NK(0,R)

Para ,1iI

  1. Gerezi=(z1i,,zKi)NK(0,R)
  2. Gereu1i,,uKiiidU(0,1)
  3. Derivar parayki=I{uki>logit1(αkxi)}log{σkzki+βkxi}1kK

Se alguém está interessado na geração a partir de posterior de dada a , esse é um problema mais difícil, embora possível pela amostragem de Gibbs ou ABC.(α,β,μ,σ,R)yki

Xi'an
fonte
11
Eu sabia que estava faltando alguma coisa. "Tudo é óbvio em retrospectiva." Minha intenção: estou interessado no valor de , então sim, estou interessado em desenhar a partir da junção posterior dos parâmetros. Eu quero que o simulado é para ver se os ajustes do modelo. yF(yi|xi)y
shadowtalker
11
Como o segundo problema é muito mais difícil? Já estimei o modelo e tenho desenhos posteriores. Podemos continuar no bate-papo, se você quiser, para não bagunçar os comentários aqui.
shadowtalker
11
Oh, em geral, sim. Felizmente, tenho Stan e o Sampler No-U-Turn fazendo o trabalho duro para mim lá.
shadowtalker
7

A resposta para a versão curta fora de contexto:

"Inverter" um cdf que não é invertível no sentido matemático (como sua distribuição mista) é viável, conforme descrito na maioria dos livros didáticos de Monte Carlo. (Como o nosso , consulte o Lema 2.4.) Se você definir o inverso generalizado então Isso significa que, quando tem um salto de em , para . Em outras palavras, se você desenhar um uniforme e ele acabar menor que , sua geração de X F  é equivalente a  X = F - ( U )  quando  U U ( 0 , 1 )

F(u)=inf{xR; F(x)u}
XF is equivalent to X=F(U) when UU(0,1).
θ y = 0F(y)θy=0u θF(u)=0uθU(0,1)θXé . Senão, quando , você acaba gerando da parte contínua, ou seja, o log-normal no seu caso. Isso significa usar uma segunda geração aleatória uniforme, , independente do primeiro desenho uniforme e definir para obter uma geração log-normal.x=0u>θvy=exp(μ+σΦ1(v))

Isso é quase o que o seu código R

Y_hat <- rbinom(1, 1, theta[i, k]) if (Y_hat == 1) Y_hat <- rlnorm(1, mu[i, k], sigma[k])

está fazendo. Você gera um Bernoulli com probabilidade e, se for igual a , transforma-o em um log-normal. Uma vez que é igual a 1 com probabilidade você deve em vez transformá-lo em uma simulação de log-normal quando é igual a zero, em vez disso, acabar com o código R modificado: 1 θ i kθki1θki

Y_hat <- rbinom(1, 1, theta[i, k])
    if (Y_hat == 0)
        Y_hat <- rlnorm(1, mu[i, k], sigma[k])
Xi'an
fonte
Assim, todos juntos, meu procedimento de simulação seria: 1) desenhar , 2) calcular , depois 3) calcular se e caso contrário. Corrigir? u k = Φzy k = 0 u kθ y k = F - 1 log-normal ( u k )uk=Φ(zk)yk=0ukθyk=Flog-normal1(uk)
shadowtalker
Não, incorreto. Você desenha um primeiro uniforme para decidir entre e log-normal, depois um segundo uniforme, caso tenha decidido por um log-normal. Veja a versão editada da minha resposta. 0
Xian
Mas isso ignora o componente ; daí a minha pergunta. Fiz uma edição esclarecedora e também resolvi o erro no meu pseudocódigo. z
shadowtalker
Minha resposta é para a versão curta e para o código R que você forneceu. Espero que ajude na versão longa, mas sua fórmula para o modelo conjunto ainda está incorreta. Você deve definir o modelo nos sem usar os uniformes ...y
Decian
Como esse modelo está incorreto? Acabei de meu na fórmula fornecida pelo artigo que citei (correspondente a em sua notação). Isso é inválido? G 1 , , G mF1,,FKG1,,Gm
shadowtalker