jssp05.Rmd

---
title: "確率モデルは面白い"
---



# 確率モデルは面白い

## 7scientist 

```{r seven socientist}
## データの準備
X <- c(-27.020,3.570,8.191,9.808,9.603,9.945,10.056)
sc7 <- list(N=NROW(X),X=X)
```
```{r,results="hide"}
## モデルコンパイル
model <- stan_model("SevenScientist.stan")
```
```{r, echo=FALSE, eval=TRUE}
paste(readLines("SevenScientist.stan"), collapse = "\n") %>% cat()
```
```{r}
## 推定
fit <- sampling(model,sc7,iter=10000)
##  表示
fit
## 描画
plot(fit,pars=c("sig[1]","sig[2]","sig[3]","sig[4]",
                "sig[5]","sig[6]","sig[7]"),show_density=T)

```

## 打ち切りデータ 

```{r censored}
## データの準備
nfails     <- 949  
n          <- 50    # Number of questions
datasasoon <- list(nF=nfails,N=n) 
```

```{r,results="hide"}
## モデルコンパイル
sasoon <- stan_model("censored.stan")
```
```{r, echo=FALSE, eval=TRUE}
paste(readLines("censored.stan"), collapse = "\n") %>% cat()
```
```{r}
fit <- sampling(sasoon,datasasoon)
fit
fit %>% as.array %>% 
  bayesplot::mcmc_dens(pars="theta") + xlim(0,1)

```

### 二項分布を確認しておこう

##### 二項分布の例

```{r}
theta <- 0.4
Ntrial <- 50
size <- 10000
data.frame(x=rbinom(size,Ntrial,theta)) %>% 
  ggplot(aes(x=x))+geom_histogram(binwidth = 0.5)
```

### 今回の例を乱数で作ってみる

```{r}
### 仮にtheta=0.4とする
N <- 10000
data.frame(X=rbinom(N,50,0.4)) %>% 
  table() %>% data.frame() %>% 
  rename(tbl=1) %>% mutate(tbl=as.numeric(levels(tbl)[tbl])) %>% 
  right_join(.,data.frame(tbl=1:50)) %>% 
  mutate(Freq=replace_na(Freq,0)) %>% 
  mutate(Cum = cumsum(Freq)/N) %>% 
  mutate(COL = ifelse(tbl<15,1,ifelse(tbl<25,2,3))) %>% 
  ggplot(aes(x=tbl,y=Cum,fill=as.factor(COL))) + 
  geom_bar(stat="identity",alpha=0.5) +
  xlab("")+ylab("")+ theme(legend.position = "none")
```


## 変化点の検出 

```{r change point}
c <- scan("changepointdata.txt")
n <- length(c)
t <- 1:n
datapoints <- list(C=c,N=n,Time=t)
```
```{r,results="hide"}
ChangeDetection <- stan_model("ChangeDetection.stan")
```
```{r, echo=FALSE, eval=TRUE}
paste(readLines("CHangeDetection.stan"), collapse = "\n") %>% cat()
```
```{r fitting CD}
fit <- sampling(ChangeDetection,datapoints)
print(fit)
```

###  可視化
```{r}
df <- transform(c)
Ms <- rstan::get_posterior_mean(fit,pars="mu")[,5]
point <- round(apply(as.matrix(rstan::extract(fit,pars="tau")$tau),2,median))
df$Mu <- c(rep(Ms[1],point),rep(Ms[2],n-point))
df %>% dplyr::mutate(num=row_number()) %>% ggplot(aes(x=num,y=X_data))+geom_line(alpha=0.5)+
  geom_point(aes(y=Mu),color="blue")

```


## 飛行機を再捕獲する

```{r plane}
x <- 10 # number of captures
k <- 4  # number of recaptures from n
n <- 5  # size of second sample
tmax <- 50 # maximum population size
datastan <- list(X=x,N=n,K=k,TMax=tmax)
```

```{r,results="hide"}
planeModel <- stan_model("plane.stan")
```
```{r, echo=FALSE, eval=TRUE}
paste(readLines("plane.stan"), collapse = "\n") %>% cat()
```
```{r}
fit <- sampling(planeModel,datastan,algorithm="Fixed_param")
fit %>% as.array %>% bayesplot::mcmc_hist(pars="t",binwidth=1)

```