Amostrador de Gibbs – probit

Hoje criei (com sucesso) meu primeiro simulador do amostrador de Gibbs para uma regressão probit e entendi tudo que fiz. A amostra era pequena e na única cadeia que rodei, coloquei como valores iniciais as estimativas de máxima verossimilhança fornecidas pelo R na função GLM (family=binomial(link=probit)). Probit é muito mais fácil de amostrar (priori conjugada) do que logística, que durante muito tempo foi preferida por mim (é mais fácil de interpretar e plotar o gráfico no R).

Como porém probit é mais fácil na Bayesiana e inclusive para utilizar um modelo hierárquico (a generalização do amostrador de Gibbs probit simples pro hierárquico é relativamente simples), que é o que vou usar na tese por causa dos dados longitudinais, estou mudando da logística para a probit.

A dica que eu dou para quem quiser usar o R para fazer seu amostrador de Gibbs é usar a função solve, do pacote Matrix. Solve é usado para resolver sistemas lineares, mas se você usar o solve só para uma matriz ele retorna a inversa dela. E nós precisamos inverter umas matrizes no amostrador de Gibbs em regressão. Aliás, nesse ponto, eu tenho uma dúvida.

Todos os livros recomendam centrar e padronizar as variáveis, pois essa parametrização aumenta a eficiência do amostrador de Gibbs, o que acredito. Porém, valores próximos de zero (que ocorerrá com padronização) tornam a matriz em geral mais mal condicionada e inclusive dificultam a utilização de Pivots para transformar a matriz e depois calcular a inversa. Essa função solve provavelmente usa a decomposição LU para inverter a matriz (já que ela chama o famoso pacote Lapack), já que é computacional caro aplicar aquele negócio que aprendemos no colégio de cofator, calcular determinante etc. Então, será que não pode dar problema nas inversões das Matrizes essa padronização?

Enfim, o caminho agora é aprender mais um pouco de Gibbs sampler no R, partir pro Metropolis-Hastings e, ao final, sair do R e fazer tudo isso no C. Pelo que vi, o mais complicado mesmo de fazer em C será carregar os dados e gerar números aleatórios de uma normal multivarida. Essa parte de inversão de Matriz tem biliotecas que fazem, como as do Lapack (que ainda naõ aprendi a usar em C). Outra complicação será chamar o C do R, para não precisar ter que carregar os dados no C e depois exportar pro R.

De todo modo, terminado esse aprendizado, estarei dominando o ferramental básico da inferência Bayesiana. Saberei fazer minhas próprias simulações de MCMC de forma eficiente computacionalmente. O passo seguinte tenho ainda de decidir: aprender técnicas mais avançadas de MCMC (como o Monte carlo Hamiltoniano?), estatística (modelos, testes etc.) ou matemática (EDO, EDP, Topologia, Jogos etc.). Sem falar nos estudos substantivos.

ps.: Essa foi a boa notícia de hoje, para contrabalançar o péssimo resultado do corinthians.

pstu: se alguém quiser o código em R, eu mando, é só pedir. Tá mal comentado ainda, mas posso ajeitar.

Anúncios

Sobre Manoel Galdino

Corinthiano, Bayesiano e Doutor em ciência Política pela USP.
Esse post foi publicado em estatística, Manoel Galdino e marcado , , , , , , , , . Guardar link permanente.

7 respostas para Amostrador de Gibbs – probit

  1. Eduardo Leoni disse:

    Já deu uma fuçada no Scythe? http://scythe.wustl.edu/ Eu cheguei a implementar um probit usando esta “library” do C++ no passado remoto. Tentei depois fazer um modelo probit hierárquico mas meus talentos de programador não se revelaram muito bons. (Hoje acho que meu problema era falta de conhecimento estatístico mesmo … mas esta é outra história.)

  2. rodolpho disse:

    Eu quero.

  3. Opa, não conhecia não. Vou dar uma olhada. Valeu.

  4. Mariane disse:

    Manoel, se puder me enviar o algoritmo eu agradeceria….to com certa dificuldade de montar o meu.

  5. Gabriel disse:

    Manoel, por gentileza, gostaria muito de ter os códigos do amostrador de Gibbs. Desde já, muito obrigado! mg2054@yahoo.com.br

  6. Michelle disse:

    Olá!
    Gostaria, por favor, do código em R do Amostrador de Gibbs.
    Muito obrigada.

  7. Débora disse:

    Olá,
    Você poderia,por favor, me enviar teu código em R do amostrador de Gibbs?
    Obrigada

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s