Eu gostaria de gerar dados aleatórios a partir de uma distribuição normal restrita usando R.
Por exemplo, eu posso querer simular uma variável de uma distribuição normal com mean=3, sd= 2
e quaisquer valores maiores que 5 são reamostrados da mesma distribuição normal.
Assim, para a função geral, eu poderia fazer o seguinte.
rnorm(n=100, mean=3, sd=2)
Eu tive alguns pensamentos:
- Iterar uma
ifelse
função com um loop que se repete até que todos os valores sejam restritos a ficar dentro dos limites. - Simule muito mais valores do que o necessário e obtenha o primeiro
n
que satisfaça a restrição. - Evite simuladores de variáveis normais vetorizadas e, em vez disso, use um loop for com um do while dentro para simular cada observação uma de cada vez e loop, quando necessário.
Todos os itens acima parecem um pouco desajeitados.
Questão
- O que é uma maneira simples de simular uma variável normal aleatória restrita em R a partir do normal com média = 3, dp = 2 e max = 5?
- De maneira mais geral, qual é uma boa maneira geral de incorporar restrições nas variáveis simuladas em R?
r
normal-distribution
simulation
truncation
Jeromy Anglim
fonte
fonte
Respostas:
Isso é chamado de distribuição normal truncada:
http://en.wikipedia.org/wiki/Truncated_normal_distribution
Christian Robert escreveu sobre uma abordagem para fazê-lo em várias situações (usando diferentes dependendo de onde estavam os pontos de truncamento) aqui:
Robert, CP (1995) "Simulação de variáveis normais truncadas",
Statistics and Computing, Volume 5, Edição 2, Junho, pp 121-125
Artigo disponível em http://arxiv.org/abs/0907.4010
Isso discute várias idéias diferentes para diferentes pontos de truncamento. Não é a única maneira de abordá-los por qualquer meio, mas geralmente apresenta um desempenho muito bom. Se você deseja fazer muitos normais truncados diferentes com vários pontos de truncamento, seria uma abordagem razoável. Como você observou,
msm::tnorm
é baseado na abordagem de Robert, enquantotruncnorm::truncnorm
implementa o amostrador de aceitação / rejeição de Geweke (1991); isso está relacionado à abordagem no trabalho de Robert. Observe quemsm::tnorm
inclui as funções densidade, cdf e quantil (cdf inverso) daR
maneira usual .Uma referência mais antiga com uma abordagem é o livro de Luc Devroye ; desde que saiu do papel, ele recuperou os direitos autorais e os disponibilizou para download.
Também há discussão de uma implementação no código R aqui (e no Rccp em outra resposta para a mesma pergunta, mas o código R é realmente mais rápido). O código R simples lá gera 50000 normais truncados em 6 milissegundos, embora esse normal truncado em particular apenas elimine as caudas extremas, portanto, um truncamento mais substantivo significaria que os resultados seriam mais lentos. Ele implementa a idéia de gerar "muitos", calculando quantos deve gerar para ter quase certeza de obter o suficiente.
Se eu precisasse apenas de um tipo específico de normal truncado muitas vezes, provavelmente tentaria adaptar uma versão do método zigurate, ou algo semelhante, ao problema.
Na verdade, parece que Nicolas Chopin já fez exatamente isso, então eu não sou a única pessoa que ocorreu:
http://arxiv.org/abs/1201.6140
Ele discute vários outros algoritmos e compara o tempo para três versões de seu algoritmo com outros algoritmos para gerar 10 ^ 8 normais aleatórios para vários pontos de truncamento.
Talvez sem surpresa, seu algoritmo acaba sendo relativamente rápido.
Edit: Um que não tenho certeza é mencionado aqui (mas talvez esteja em um dos links) é transformar (via cdf normal inverso) um uniforme truncado - mas o uniforme pode ser truncado simplesmente gerando um uniforme dentro dos limites de truncamento . Se o cdf normal inverso for rápido, isso é rápido e fácil e funciona bem para uma ampla variedade de pontos de truncamento.
fonte
Seguindo as referências do @ glen_b e focando exclusivamente na implementação de R. Existem algumas funções projetadas para amostrar a partir de uma distribuição normal truncada:
rtruncnorm(100, a=-Inf, b=5, mean=3, sd=2)
notruncnorm
pacotertnorm(100, 3, 2, upper=5)
nomsm
pacotefonte