ベイズ統計学入門　渡部洋先生　 福村出版

渡部　洋著　「ベイズ統計学入門」の一部をプログラムしてみました。

```{r}

#渡部　洋著　「ベイズ統計学入門」


library(dplyr)
library(pforeach)

#IQの設定値
ave_IQ=100;var_IQ=15

#現在の標本に関する標本数と平均
n=40;ave_y=104.3

#過去10年間の平均ＩＱと分散ＩＱ
past_10_years_ave=105.6;past_10_years_var=3.88

#現在の標本をもとに母平均の平均値と分散を計算
con_ave=((n/var_IQ^2)*ave_y+past_10_years_ave/past_10_years_var)/(n/var_IQ^2+1/past_10_years_var)

con_var=1/(n/var_IQ^2+1/past_10_years_var)

#実際に現在のＩＱ値の予測値サンプルを発生させ、区間推定を行う。
N=100;
IQ=rnorm(N,mean=con_ave,sd=sqrt(con_var+var_IQ))
plot(IQ)

alpha=0.05
under=con_ave+qnorm(alpha)*sqrt(con_var+var_IQ)
upper=con_ave+qnorm(1-alpha)*sqrt(con_var+var_IQ)
print(paste0("estimated interval is (",under,",",upper,")"))

#母平均のサンプルを発生させ、その区間推定を行う。
N=100;
ave=rnorm(N,mean=con_ave,sd=sqrt(con_var))
plot(ave)

alpha=0.05
under_ave=con_ave+qnorm(alpha)*sqrt(con_var)
upper_ave=con_ave+qnorm(1-alpha)*sqrt(con_var)
print(paste0("estimated interval is (",under_ave,",",upper_ave,")"))

#現在ある標本データを基準にした、母平均ＩＱが100を超える確率
value=(100-con_ave)/sqrt(past_10_years_var)
if(value<0){
  prob=1-pnorm(value)
}else{
  prob=pnorm(value)
}
print(paste0("母平均ＩＱが100を超える確率は ",prob))

#現在ある標本データを基準にした、次の予測ＩＱ値が100を超える確率
value=(100-con_ave)/sqrt(var_IQ^2+con_var)
if(value<0){
  prob=1-pnorm(value)
}else{
  prob=pnorm(value)
}
print(paste0("予測値ＩＱが100を超える確率は ",prob))



```


```{r}
library(dplyr)
library(pforeach)

data=data.frame(students=1:5,train=c(1,2,2,3,4),test=c(5,10,12,16,18)) %>% mutate(train_prob=train/nrow(data),test_prob=test/(nrow(data)*4)) %>% mutate(train_g=asin(sqrt((train+3/8)/(nrow(data)+3/4))),test_g=asin(sqrt((test+3/8)/(nrow(data)*4+3/4))))

PROB=data.frame(students=1:5,g=data$test_g,var=1/(4*nrow(data)*4+2),prior_prob_ave=data$train_g,prior_prob_var=1/(4*1+2),posterior_prob_ave=(data$test_g*(4*4*nrow(data)+2)+data$train_g*(4*1+2))/(4*4*nrow(data)+2+4*1+2),posterior_prob_var=1/(4*4*nrow(data)+2+4*1+2))

#母平均が他の生徒より大きくなる確率
i=1;j=3 

p=1-pnorm((PROB$posterior_prob_ave[j]-PROB$posterior_prob_ave[i])/sqrt(2*PROB$posterior_prob_var[1]))

N=1000
val1=rnorm(N,mean=PROB$posterior_prob_ave[j],sd=sqrt(PROB$posterior_prob_var[1]))
val2=rnorm(N,mean=PROB$posterior_prob_ave[i],sd=sqrt(PROB$posterior_prob_var[1]))

max_val=max(val1,val2);min_val=min(val1,val2)
plot(val1,xlab="sample",ylab="value",xlim=c(0,(N+1)),ylim=c(min_val,max_val),main=paste0("green;",j,", red;",i),col=3)
par(new=T)
plot(val2,xlab="sample",ylab="value",xlim=c(0,(N+1)),ylim=c(min_val,max_val),main=paste0("green;",j,", red;",i),col=2)

print(paste0("サンプル確率＝",length(val1[val2>val1])/N))
print(paste0("実際の確率＝",p))

mean_g=mean(data$test_g)

#事前分布が互いに独立にN(u,1/6)に従うとした場合
normal_data=data.frame(g=data$test_g,posterior_mean=(data$test_g*(4*4*nrow(data)+2)+mean_g*(4*1+2))/(4*1+2+2+4*4*nrow(data)))

```

```{r}
library(dplyr)
library(pforeach)

#逆カイ二乗分布

#標本の大きさ
n=16

#標本分散
s2=28.09

alpha=0.05

inv_chisq=function(u){
  lambda=n*s2;df=n
  y<-(sqrt(lambda^df)/(sqrt(2^df)*gamma(df/2)))*(u^(-(df/2+1)))*exp(-lambda/(2*u))
  return(y)
}
N=100;sep=100
inv_chisq_data=data.frame(value=seq(0,N,by=1/sep))
inv_chisq_data=inv_chisq_data%>% dplyr::mutate(inv_chisq=inv_chisq(inv_chisq_data$value))
inv_chisq_data[is.na(inv_chisq_data)>0]=0

plot(inv_chisq_data)

inv_chisq_data=inv_chisq_data %>% mutate(p_inv_chisq=0)

for(j in 1:nrow(inv_chisq_data)){
  value=inv_chisq_data$value[j]
  if(value>0){
  inv_chisq_data$p_inv_chisq[j]=integrate(inv_chisq,0,inv_chisq_data$value[j])$value
}
}

inv_chisq_data$value[abs(inv_chisq_data$p_inv_chisq-(1-alpha))==min(abs(inv_chisq_data$p_inv_chisq-(1-alpha)))]


```