# Stanサンプリング

```---
title: "StanDistribution"
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)
```
```