Stanサンプリング

• 参考サイト

---
title: "StanDistribution"
author: "ryamada"
date: "2017年6月10日"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(rstan)
```


# Stanの各種分布の練習

```{r}
.stan_code1 = '
data {
  int n_obs;
  vector [n_obs] x;
}

parameters {
  real mu;
  real<lower=0> sigma;
}

model {
  x ~ normal(mu, sigma);
}'
```
```{r,echo=FALSE}
#Compiling
MODEL1 = stan_model(model_code = .stan_code1)
```

```{r,echo=FALSE}
.stan_code2 = '
data {
  int n_obs;
  vector [n_obs] x;
}

parameters {
  real mu;
  real<lower=0> sigma;
}

model {
  x ~ cauchy(mu, sigma);
}'
#Compiling
MODEL2 = stan_model(model_code = .stan_code2)
```
```{r,echo=FALSE}
.stan_code3 = '
data {
  int n_obs;
  vector [n_obs] x;
}

parameters {
  real a;
  real b;
}

model {
  x ~ gamma(a,b);
}'
#Compiling
MODEL3 = stan_model(model_code = .stan_code3)
```
```{r,echo=FALSE}
.stan_code4 = '
data {
  int n_obs;
  int x[n_obs];
}

parameters {
  real <lower=0> a;
}

model {
  x ~ poisson(a);
}'
#Compiling
MODEL4 = stan_model(model_code = .stan_code4)
```

```{r,echo=FALSE}
.stan_code5 = '
data {
  int n_obs;
  int x[n_obs];
  int ns[n_obs];
}

parameters {
  real <lower=0> p;
}

model {
  x ~ binomial(ns, p);
}'
#Compiling
MODEL5 = stan_model(model_code = .stan_code5)
```
```{r, echo=FALSE}
.stan_code6 = '
data {
  int n_obs;
  int K;
  int<lower=1,upper=K> Z[n_obs];  # data: category index
  vector<lower=0>[K] alpha;
}

parameters {
  simplex[K] theta; 
}

model {
  for (i in 1:n_obs)
   Z[i] ~ categorical(theta);
  theta ~ dirichlet(alpha); 
}'
#Compiling
MODEL6 = stan_model(model_code = .stan_code6)

```
```{r,echo=FALSE}
.data1 = list(x=rnorm(10000, 10, 3))
.data1$n_obs = length(.data1$x)
#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit1 = sampling(MODEL1, data=.data1)


EXT1 = extract(.fit1)
```
```{r,echo=FALSE}
.data2 = list(x=rcauchy(10000, 10, 3))
.data2$n_obs = length(.data2$x)
#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit2 = sampling(MODEL2, data=.data2)


EXT2 = extract(.fit2)
```
```{r,echo=FALSE}
.data3 = list(x=rgamma(10000, 2,5))
.data3$n_obs = length(.data3$x)
#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit3 = sampling(MODEL3, data=.data3)


EXT3 = extract(.fit3)
```
```{r,echo=FALSE}
.data4 = list(x=rpois(10000, 2.5))
.data4$n_obs = length(.data4$x)
#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit4 = sampling(MODEL4, data=.data4)


EXT4 = extract(.fit4)
```

```{r,echo=FALSE}
.data5 = list(ns=rpois(10000, 2.5))
p = 0.3
.data5$x = rbinom(length(.data5$ns),.data5$ns,p)
.data5$n_obs = length(.data5$x)

#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit5 = sampling(MODEL5, data=.data5)


EXT5 = extract(.fit5)
```
```{r,echo=FALSE}
K <- 5
library(MCMCpack)
r <- rdirichlet(1,rep(1,K))
.data6 = list(Z=sample(1:K,10000,replace=TRUE,r))
.data6$K = K
.data6$n_obs = length(.data6$Z)
.data6$alpha = rep(1,K)

#stan.out <- stan(file="SNPcovariate.stan",data=dat)
.fit6 = sampling(MODEL6, data=.data6)


EXT6 = extract(.fit6)
```
```{r}
print(.fit1)
summary(.fit1)
plot(.fit1)
pairs(.fit1)
rstan::traceplot(.fit1)
rstan::stan_trace(.fit1)
rstan::stan_hist(.fit1)
rstan::stan_dens(.fit1)
```
```{r}
print(.fit2)
summary(.fit2)
plot(.fit2)
pairs(.fit2)
rstan::traceplot(.fit2)
rstan::stan_trace(.fit2)
rstan::stan_hist(.fit2)
rstan::stan_dens(.fit2)
```
```{r}
print(.fit3)
summary(.fit3)
plot(.fit3)
pairs(.fit3)
rstan::traceplot(.fit3)
rstan::stan_trace(.fit3)
rstan::stan_hist(.fit3)
rstan::stan_dens(.fit3)
```
```{r}
print(.fit4)
summary(.fit4)
plot(.fit4)
pairs(.fit4)
rstan::traceplot(.fit4)
rstan::stan_trace(.fit4)
rstan::stan_hist(.fit4)
rstan::stan_dens(.fit4)
```
```{r}
print(.fit5)
summary(.fit5)
plot(.fit5)
pairs(.fit5)
rstan::traceplot(.fit5)
rstan::stan_trace(.fit5)
rstan::stan_hist(.fit5)
rstan::stan_dens(.fit5)
```

```{r}
print(.fit6)
summary(.fit6)
plot(.fit6)
pairs(.fit6)
rstan::traceplot(.fit6)
rstan::stan_trace(.fit6)
rstan::stan_hist(.fit6)
rstan::stan_dens(.fit6)
```