一个R软件新的做Bayesian分析的包

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

名字叫“rethinking” ，可以作为RStan的前端使用，自动计算DIC和WAIC，非常方便。

rethinking

---

      This R package accompanies a course and book on Bayesian data analysis. It contains tools for conducting both MAP estimation and Hamiltonian Monte Carlo (through RStan). These tools force the user to specify the model as a list of explicit distributional assumptions.
For example, a simple Gaussian model could be specified with this list of formulas:

1. f <- alist(
   
2.     y ~ dnorm( mu , sigma ),
   
3.     mu ~ dnorm( 0 , 10 ),
   
4.     sigma ~ dcauchy( 0 , 1 )
   
5. )

复制代码

      The first formula in the list is the likelihood; the second is the prior for mu; the third is the prior for sigma (implicitly a half-Cauchy, due to positive constraint on sigma).


MAP estimation

---

      Then to use maximum a posteriori (MAP) fitting:

1. library(rethinking)
   
2. fit <- map(
   
3.     f ,
   
4.     data=list(y=c(-1,1)) ,
   
5.     start=list(mu=0,sigma=1)
   
6. )

复制代码

      The object fit holds the result.




Hamiltonian Monte Carlo estimation

---

            The same formula list can be compiled into a Stan (mc-stan.org) model:

1. fit.stan <- map2stan(
   
2.     f ,
   
3.     data=list(y=c(-1,1)) ,
   
4.     start=list(mu=0,sigma=1)
   
5. )

复制代码

      The start list is optional, provided a prior is defined for every parameter. In that case, map2stan will automatically sample from each prior to get starting values for the chains. The chain runs automatically, provided rstan is installed. The Stan code can be accessed by using stancode(fit.stan):

1. data{
   
2.     int<lower=1> N;
   
3.     real y[N];
   
4. }
   
5. parameters{
   
6.     real mu;
   
7.     real<lower=0> sigma;
   
8. }
   
9. model{
   
10.     mu ~ normal( 0 , 10 );
    
11.     sigma ~ cauchy( 0 , 1 );
    
12.     y ~ normal( mu , sigma );
    
13. }
    
14. generated quantities{
    
15.     real dev;
    
16.     dev <- 0;
    
17.     dev <- dev + (-2)*normal_log( y , mu , sigma );
    
18. }

复制代码

   
Multilevel model formulas

---

      While map is limited to fixed effects models for the most part, map2stan can specify multilevel models, even quite complex ones. For example, a simple varying intercepts model looks like:

1. f2 <- alist(
   
2.     y ~ dnorm( mu , sigma ),
   
3.     mu <- a + aj,
   
4.     aj[group] ~ dnorm( 0 , sigma_group ),
   
5.     a ~ dnorm( 0 , 10 ),
   
6.     sigma ~ dcauchy( 0 , 1 ),
   
7.     sigma_group ~ dcauchy( 0 , 1 )
   
8. )

复制代码

And with varying slopes as well:

1. f3 <- alist(
   
2.     y ~ dnorm( mu , sigma ),
   
3.     mu <- a + aj + (b + bj)*x,
   
4.     c(aj,bj)[group] ~ dmvnorm( 0 , Sigma_group ),
   
5.     a ~ dnorm( 0 , 10 ),
   
6.     b ~ dnorm( 0 , 1 ),
   
7.     sigma ~ dcauchy( 0 , 1 ),
   
8.     Sigma_group ~ inv_wishart( 3 , diag(2) )
   
9. )

复制代码

   Nice covariance priors

---

      And map2stan supports decomposition of covariance matrices into vectors of standard deviations and a correlation matrix, such that priors can be specified independently for each:

1. f4 <- alist(
   
2.     y ~ dnorm( mu , sigma ),
   
3.     mu <- a + aj + (b + bj)*x,
   
4.     c(aj,bj)[group] ~ dmvnorm2( 0 , sigma_group , Rho_group ),
   
5.     a ~ dnorm( 0 , 10 ),
   
6.     b ~ dnorm( 0 , 1 ),
   
7.     sigma ~ dcauchy( 0 , 1 ),
   
8.     sigma_group ~ dcauchy( 0 , 1 ),
   
9.     Rho_group ~ dlkjcorr(2)
   
10. )

复制代码

    Semi-automated Bayesian imputation

---

          It is possible to code simple Bayesian imputations this way. For example, let's simulate a simple regression with missing predictor values:

1. N <- 100
   
2. N_miss <- 10
   
3. x <- rnorm( N )
   
4. y <- rnorm( N , 2*x , 1 )
   
5. x[ sample(1:N,size=N_miss) ] <- NA

复制代码

That removes 10 x values. Then the map2stan formula list just defines a distribution for x:

1. f5 <- alist(
   
2.     y ~ dnorm( mu , sigma ),
   
3.     mu <- a + b*x,
   
4.     x ~ dnorm( mu_x, sigma_x ),
   
5.     a ~ dnorm( 0 , 100 ),
   
6.     b ~ dnorm( 0  , 10 ),
   
7.     mu_x ~ dnorm( 0 , 100 ),
   
8.     sigma_x ~ dcauchy(0,2),
   
9.     sigma ~ dcauchy(0,2)
   
10. )
    
11. m5 <- map2stan( f5 , data=list(y=y,x=x) )

复制代码

       What map2stan does is notice the missing values, see the distribution assigned to the variable with the missing values, build the Stan code that uses a mix of observed and estimated x values in the regression. See the stancode(m) for details of the implementation.

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏6 回帖

关键词：Bayesian分析 Bayesian Bayes baye bay rethink Bayesian分析

Nicolle 发表于 2015-5-10 03:54

谢谢分享，欢迎多多分享此类好帖。

Nicolle 发表于 2015-5-10 03:54

多谢分享！

多谢分享！

可以的，library名字也很有意思

多谢楼主分享~~

好人~~~~~~