ノンパラメトリックベイズ　講談社　佐藤一誠先生

青い本の機械学習プロフェッショナルシリーズ「ノンパラメトリックベイズ」の第2,3章をプログラムしました。
少しおかしいところを徐々に修正していきます。

```{r}

#ノンパラメトリックベイズ


library(dplyr)
library(pforeach)

n=1000;N=50;p=0.6

sample=rbinom(n,size=N,prob=p)

sample_data=data.frame(num=sort(unique(sample)),count=0)

for(j in 1:nrow(sample_data)){
  sample_data$count[j]=length(sample[sample==sample_data$num[j]])
}

sample_data=sample_data %>% mutate(prob=count/sum(count))
plot(sample_data$num,sample_data$count,type="o",xlab="num",ylab="prob")

#正規分布で近似

normalization=data.frame(u=sort(unique(sample)),KL=0)

for(j in 1:nrow(normalization)){
  u=normalization$u[j]
  normalization$KL[j]=mean(sample_data$prob*log(sample_data$prob))-mean(sample_data$prob*log((1/sqrt(2*pi))*exp(-((sample_data$num-u)^2)/2)))
}

#最尤推定量(平均のみ：u) p(x|u)

#KL情報量の一番小さい値をとる平均値uを推定量とする。
u=normalization$u[normalization$KL==min(normalization$KL)]


plot(sample_data$num,sample_data$prob,type="o",xlim=c(min(sample_data$num),max(sample_data$num)),ylim=c(min(sample_data$prob),max(dnorm(sample_data$num,mean=u,sd=1))),xlab="num",ylab="prob",col=2,main="green;predict, red;sample")
par(new=T)
plot(sample_data$num,dnorm(sample_data$num,mean=u,sd=1),type="o",xlim=c(min(sample_data$num),max(sample_data$num)),ylim=c(min(sample_data$prob),max(dnorm(sample_data$num,mean=u,sd=1))),xlab="num",ylab="prob",col=3,main="green;predict, red;sample")

#最尤推定量(平均と分散:u,Q) p(x|u,Q)

#正規分布の最尤推定量はu=mean(sample),Q=var(sample)であることより、これらの値がKL情報量を最小にする値になる。

u=mean(sample);Q=var(sample)

KL=mean(sample_data$prob*log(sample_data$prob))-mean(sample_data$prob*log(dnorm(sample_data$num,mean=u,sd=sqrt(Q))))

plot(sample_data$num,sample_data$prob,type="o",xlim=c(min(sample_data$num),max(sample_data$num)),ylim=c(min(sample_data$prob),max(dnorm(sample_data$num,mean=u,sd=sqrt(Q)),sample_data$prob)),xlab="num",ylab="prob",col=2)
par(new=T)
plot(sample_data$num,dnorm(sample_data$num,mean=u,sd=sqrt(Q)),type="o",xlim=c(min(sample_data$num),max(sample_data$num)),ylim=c(min(sample_data$prob),max(dnorm(sample_data$num,mean=u,sd=sqrt(Q)),sample_data$prob)),xlab="num",ylab="prob",col=3)

#p(u|h)がポアソン分布に従うとすると、hはそのパラメータ値となる。今、u=argmax{log(p(u|h))+sum(log(p(x|u,Q)))}_{x:sample}を考えると、sum(log(p(x|u,Q)))を最大にするu,Qは正規分布に従うとしてすでに求まっている。もしこの値を最大にするhを計算するならば、それはh=mean(sample)である。(ポアソン分布の最尤法)

#log(p(u|h))+sum(log(p(x|u,Q)))

data=data.frame(num=sample_data$num,value=0)
for(j in 1:nrow(data)){
  num=sample_data$num[j]
  data$value[j]=sum(log(dnorm(sample_data$num,mean=u,sd=sqrt(Q))))+log((exp(-num)*(num^u))/gamma(u+1))
}

plot(data$num,data$value,type="o",main="param:average",xlab="parameter values",ylab="values")

h=round(mean(sample))

print(paste0("パラメータＱ(分散)を既知としたとき、log(p(u|h))+sum(log(p(x|u,Q)))を最大にするパラメータuは",h,"で、その値は",data$value[data$num==h]))

#式(2.18)において、Ｑ(分散)が既知の場合(predict1)；sum(p(X|u,Q)*p(u|{x:sample}))

sample_data=sample_data %>% dplyr::mutate(predict1=0,predict2=0)

for(k in 1:nrow(sample_data)){

num=sample_data$num[k]  

prob=0

for(j in 1:nrow(sample_data)){

prob=prob+dnorm(num,mean=sample_data$num[j],sd=sqrt(Q))*dpois(sample_data$num[j],lambda=h)

}

sample_data$predict1[k]=prob

}

#p(x*|{x;sample})をプロット

P_x=data.frame(num=seq(0,100,1),prob=0)

for(k in 1:nrow(P_x)){

num=P_x$num[k]  

prob=0

for(j in 1:nrow(P_x)){

prob=prob+dnorm(num,mean=P_x$num[j],sd=sqrt(Q))*dpois(P_x$num[j],lambda=h)

}

P_x$prob[k]=prob

}

plot(P_x$num,P_x$prob,col=4,xlab="num",ylab="prob",main="p(x*|{x;sample})")



#式(2.19)において、Ｑ(分散)が既知の場合(predict2)；sum(p(X|u,Q)/S), u~p(u|{x:sample}))

S=1000

for(j in 1:nrow(sample_data)){
  
  random_num=rpois(S,lambda=h)
  
  num=sample_data$num[j]
  
  prob=0
  
  for(k in 1:S){
    
  prob=prob+dnorm(num,mean=random_num[k],sd=sqrt(Q))/S  
  }
 
sample_data$predict2[j]=prob  
   
}

#式(2.15)

prob_data=data.frame(num=data$num,prob=0)

for(j in 1:nrow(prob_data)){

prob_data$prob[j]=prod(dnorm(data$num,mean=prob_data$num[j],sd=sqrt(Q)))*dpois(prob_data$num[j],lambda=h)/prod(dpois(data$num,lambda=h))

}




```

```{r}

library(dplyr)
library(pforeach)
library(DirichletReg)
library(MCMCpack)

#Dirichlet分布やBeta分布はある一点の密度が１を超える分布関数です。

#dirichlet dis(density function)

dirichlet=function(x,a){
  y<-(gamma(sum(a))/(prod(gamma(a))))*prod(x^(a-1))
  return(y)
}

  #sample(x=c(0.1,0.9),a=c(2,3))

  print(paste0("自身で作成したもの;",dirichlet(c(0.1,0.9),c(2,3))))
  
  #use package
  
  print(paste0("パッケージを使ったもの;",ddirichlet(c(0.1,0.9),c(2,3))))
  
  
#beta dis(density function)

beta=function(x){
  a=c(3,4)
  y<-(x^(a[1]-1))*(1-x)^(a[2]-1)/(1/(gamma(sum(a))/(prod(gamma(a)))))
  return(y)
}

plot(seq(0,1,by=0.01),beta(seq(0,1,by=0.01)),xlab="x",ylab="dencity",type="o",main="Beta")


print(paste0("全体の確率は ",integrate(beta,0,1)$value))

multi=function(x,n){
  y<-(gamma(sum(n)+1)/prod(gamma(n)+1))*prod(x^(n))
  return(y)
}

#p.24 サンプルの作成

dice=6 #diceの目の数

N=100 #標本数

dice=sample(1:dice,N,replace=TRUE) #z_{1:N}

#z_{1:N}のサンプルから、目の出た頻度を計算

dice_table=data.frame(values=1:max(dice),count=0,freq=0)

for(j in 1:nrow(dice_table)){
  dice_table$count[j]=length(dice[dice==dice_table$values[j]])
  dice_table$freq[j]=length(dice[dice==dice_table$values[j]])/N
}

#p(q|z_{1:N},a)をp_dataのdensity_probとして計算

#Dirichlet分布のalpha(a)を仮定

a=sample(length(unique(dice_table$values)))

ite=1000
p_data=pforeach(k=1:ite,.c=rbind)({
  
  
#q=(q1,q2,...)をランダムに作成(0<=qi<=1,i=1,2,...)
N=100;sample=sample(1:max(dice),N,replace=TRUE)

sample_data=data.frame(num=1:max(dice),freq=0)
for(j in 1:nrow(sample_data)){
  sample_data$freq[j]=length(sample[sample==sample_data$num[j]])/N
}
n=dice_table$count

result=data.frame(num=k,q=t((sample_data$freq)),density_prob=0)

result$density_prob[1]=(gamma(sum(a))/(prod(gamma(a))))*prod(sample_data$freq^(n+a-1))

result
})

#z_{1:N}のデータとalpha(a)をもとに、最もらしいq=(q1,q2,..)をもとめると；

best_data=p_data %>% filter(density_prob==max(density_prob))

print(paste0("もっともらしいq=(q1,q2,...)のあたいは;"))
print(c(best_data[1,(1+1):(1+length(a))]))             
print(paste0("で確率は;",best_data[1,ncol(best_data)]))

#p(z_{N+1}=k|z_{1:N},a)を計算する

P_data=data.frame(num=sort(unique(dice_table$values)),prob=(dice_table$count+a)/sum(dice_table+a))

#この確率をもとにさいころの目の乱数を発生
samples=sample(sort(unique(dice_table$values)),N,replace=TRUE,prob=P_data$prob)

plot(dice,type="p",col=2,xlab="num",ylab="A roll of the dice",main="red:samples,green:predict")
par(new=T)
plot(samples,type="p",col=3,xlab="num",ylab="A roll of the dice",main="red:samples,green:predict")
```

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

#gamma-gauss dis model

#ガウス分布の事後分布を求める。平均事後分布と分散の事後分布、平均と分散の事後分布を求める。

#sample datas
N=10
n=100

sample=rbinom(n,size=N,prob=0.6)
x=sample
uniq_num=sort(unique(sample))
sample_data=data.frame(num=uniq_num,count=0)

for(j in 1:nrow(sample_data)){
  sample_data$count[j]=length(sample[sample==uniq_num[j]])
}


plot(sample_data$num,sample_data$count,type="o")

#p(u|x,Q,u0,Q0), Q:分散,u:平均

#事前分布のparameter, u0とQ0を仮定

u0=6;Q0=10 #二項分布の平均が6なので6とするときれいに映る

#分散値(共分散行列で1次元の場合で、p.26の3.4.1より固定値なので指定する。)

Q=10

#p(u|sample_data$num,Q,u0,Q0)を計算--(3.20)

#u_ave=(母平均分散/(標本分散+母平均分散))*(標本平均)＋(標本分散/(標本分散+母平均分散))*(母平均値)=(母平均分散割合)*(標本平均)+(標本平均分散割合)*(母平均値)=(分散割合を確率とした)平均値の平均

u_ave=(Q0/((Q/n)+Q0))*mean(x)+((Q/n)/((Q/n)+Q0))*u0

u_var=(1/Q0+n/Q)^(-1)

#p(u|sample_data$num,Q,u0,Q0)の９０％区間推定

alpha=0.05

upper=u_ave+qnorm(1-alpha)*sqrt(u_var)/sqrt(n)
under=u_ave+qnorm(alpha)*sqrt(u_var)/sqrt(n)

print(paste0("u の90%区間は(",under," , ",upper,")"))


#p(x_n+1|x,Q,u0,Q0)を計算--(3.23)

x_ave=u_ave

x_var=Q+u_var #標本分散＋母平均分散

#p(x_n+1|x,Q,u0,Q0)の９０％区間推定

alpha=0.05

upper=x_ave+qnorm(1-alpha)*sqrt(x_var)/sqrt(n)
under=x_ave+qnorm(alpha)*sqrt(x_var)/sqrt(n)

print(paste0("x_n+1 の90%区間は(",under," , ",upper,")"))

#実際のサンプルの値と予測値の比較

Sample=data.frame(num=1:length(sample),sample=sample,predict_x=0)

for(j in 1:nrow(Sample)){
  Sample$predict_x[j]=rnorm(1,mean=x_ave,sd=sqrt(x_var))
}
plot(Sample$num,Sample$sample,type="p",xlab="num",ylab="value",xlim=c(min(Sample$num),max(Sample$num)),ylim=c(min(Sample$sample,Sample$predict_x),max(Sample$sample,Sample$predict_x)),col=2,main="red:samples, green:predict_x")
par(new=T)
plot(Sample$num,Sample$predict_x,type="p",xlab="num",ylab="value",xlim=c(min(Sample$num),max(Sample$num)),ylim=c(min(Sample$sample,Sample$predict_x),max(Sample$sample,Sample$predict_x)),col=3,main="red:samples, green:predict_x")
```

```{r}

library(dplyr)
library(pforeach)
library(DirichletReg)
library(MCMCpack)

#3.4.2 逆ガンマ分布

#平均値を固定、分散が確率変数の場合

#sampleを作成

N=100

x=rnorm(N,mean=2,sd=10)

#分散Qが確率変数として、事前分布が逆ガンマ分布に従う

#p(Q|x_{1:N},u,a,b)を計算する

#事前分布p(Q|a,b)

Inv_gamma=function(Q,a,b){
  
  y<-((b^a)/gamma(a))*(Q^(-a-1))*exp(-b/Q)
  
  return(y)
  
}

#p(x_{1:N}|u,Q)

normal=function(x,u,Q){
  
  y<-prod((1/sqrt(2*pi*Q))*(exp(-((x-u)^2)/(2*Q))))
  
  return(y)
  
}

#データの値から、normalの平均と分散を予想

#最尤推定量

u=mean(x);Q=var(x)

plot(x,xlab="num",ylab="values",type="o",col=2,ylim=c(-30,30))

par(new=T)

plot(rnorm(length(x),mean=u,sd=sqrt(Q)),xlab="num",ylab="values",type="o",col=3,ylim=c(-30,30))

#Q=b/(a+1)になるようにa,bを求める(対数尤度法)

#a,bを仮定して事前分布を計算

a=1;b=2*Q

plot(seq(0,200,by=0.1),Inv_gamma(seq(0,200,by=0.1),a,b),type="o",xlab="X",ylab="d_inv_gamma",col=2)

inv_gamma=data.frame(X=seq(0,200,by=0.1),inv_gamma=Inv_gamma(seq(0,200,by=0.1),a,b))

inv_gamma[is.na(inv_gamma)>0]=0

max_density=max(inv_gamma$inv_gamma)

max_density_X=inv_gamma$X[inv_gamma$inv_gamma==max(inv_gamma$inv_gamma)]

print(paste0("このことから、Q=",max_density_X,"と予想"))

print(paste0("p(Q|x_{1:N},u,a,b)=",normal(x,u,Q)*Inv_gamma(Q,a,b)))

Prob_data=data.frame(Q=seq(Q-20,Q+20,by=0.1),prob=0)

for(j in 1:nrow(Prob_data)){
  
  q=Prob_data$Q[j]
  
  Prob_data$prob[j]=normal(x,u,q)*Inv_gamma(q,a,b)
  
}

#p(x_{N+1}|x_{1:N},u,a,b)--(3.30)

#t=1/Q(精度)を事前分布のパラメータとしてもとめる。なぜなら、

inv_gamma_samples<-rinvgamma(200,shape=10,scale=5)

gamma_samples<-rgamma(200,shape=10,scale=1/5)

plot(1/gamma_samples,xlab="num",ylab="values",type="o",col=2,main="red:1/gamma, green:inv_gamma",ylim=c(min(inv_gamma_samples,1/gamma_samples),max(inv_gamma_samples,1/gamma_samples)))

par(new=T)

plot(inv_gamma_samples,xlab="num",ylab="values",type="o",col=3,main="red:1/gamma, green:inv_gamma",ylim=c(min(inv_gamma_samples,1/gamma_samples),max(inv_gamma_samples,1/gamma_samples)))

print(paste0("1/gamma の乱数値の平均は",mean(1/gamma_samples),"で分散は",var(1/gamma_samples)))

print(paste0("inv_gamma の乱数値の平均は",mean(inv_gamma_samples),"で分散は",var(inv_gamma_samples)))

t=1/Q;

x=seq(0,1000,by=0.1)

plot(x,dgamma(x,shape=a,scale=b),xlab="x",ylab="d_gamma",type="o",col=3)


```

```{r}
#gamma分布調査
library(dplyr)
library(MCMCpack)
library(pforeach)
x=rpois(n,lambda=2)
a=-42;b=-100.695;n=100

a=a+n/2;b=b+(1/2)*sum((x-mean(x))^2)
ave=a/b
var=a/(b^2)
print(ave)
print(var)
plot(rgamma(n,shape=a,scale=1/b),type="o")
plot(x,type="o")

```

```{r}
#p.31

library(dplyr)
library(pforeach)

#IQデータ（平均；100、分散；15）

Average=100;
Variance=15;
Precision=1/Variance

#実際のサンプルデータ作成


samples=c(rnorm(10,mean=101,sd=10^(1/2)),rnorm(10,mean=98,sd=8^(1/2)))

data=data.frame(num=1:length(samples),samples=samples)

average=mean(samples);variance=var(samples);precision=1/variance

b=Average/(1/Variance);a=(1/Variance)*b

u0=Average

gamma1=function(t,a,b){
  
  y<-(b/gamma(a+1))*(prod(rep(b*t,a-1)))*exp(-b*t)
  
  return(y)
  
}

normal1=function(u,u0,t,p){
  
  y<-((p*t/(2*pi))^(1/2))*exp(-((p*t*(u-u0)^2)/2))
  
  return(y)
  
}

normal2=function(X,u,t){
  
  y<-((t/(2*pi))^(length(X)/2))*exp(-(t/2)*sum((X-u)^2))
  
  return(y)
  
}

gamma1_data=data.frame(t=seq(0,a/b+1,by=0.0001),value=0)

for(j in 1:nrow(gamma1_data)){
  
  gamma1_data$value[j]=gamma1(gamma1_data$t[j],a,b)
  
}

gamma1_data[is.na(gamma1_data)>0]=0

plot(gamma1_data$t,gamma1_data$value,type="l")

est_t=gamma1_data$t[gamma1_data$value==max(gamma1_data$value)]

normal1_data=data.frame(u=seq(Average-10,Average+10,by=0.1),value=0)

for(j in 1:nrow(normal1_data)){
  
  normal1_data$value[j]=normal1(normal1_data$u[j],u0,est_t,1)
  
}


plot(normal1_data$u,normal1_data$value,type="l")

est_u=normal1_data$u[normal1_data$value==max(normal1_data$value)]


#p(u,t|X,u0,p,a,b)を計算(式(3.35)を使用)

prob=function(z,v,num){
  p=1;n=length(samples)
  ave=(n/(n+p))*mean(samples)+(p/(p+n))*Average
  var=v*(p+n)
  A=a+n/2;
  B=b+(1/2)*sum((samples-mean(samples))^2)+(n*p*(mean(samples)-Average)^2)/(2*(p+n))
  B=1/B
  normal=(1/sqrt(2*pi*var))*exp(-((z-ave)^2)/(2*var))
  gamma=((B^(-1))/(gamma(A+1)))*prod(rep(v/B,A-1))*exp(-v/B)
  
  if(num==1){
  y<-normal*gamma
  print(paste0("平均が",average,"で,精度が",precision,"である確率p(u,t|X,u0,p,a,b)を導出"))
  
  }else{
    if(num==2){
      y<-normal
  print(paste0("平均が",average,"で,精度が",precision,"である確率p(u|X,t,u0,p)を導出"))    
    }else{
     y<-gamma
    print(paste0("精度が",precision,"である確率p(t|X,u0,p,a,b)を導出"))
      } }
  return(y)
}

print(prob(average,precision,1))
#p(t|X,u0,p,a,b)を計算(式(3.37))
print(prob(average,precision,0))
#p(u|X,t,u0,p)を計算(式(3.36))
print(prob(average,precision,2))
print(paste0("基準値は平均",Average,";精度",Precision))

```

```{r}
#x=u,y=t,z=p(u,t|X,u0,p,a,b)を軸としてプロット

z_data=data.frame(x=0,y=0,z=0)
width=3;sep=0.1
x=seq(Average-width,Average+width,by=sep)
width=0.02;sep=0.001
for(j in 1:length(x)){
y=seq(Precision-width,Precision+width,by=sep)
Z_data=data.frame(x=x[j],y=y,z=0)

for(l in 1:nrow(Z_data)){
   p=1;n=length(samples);z=Z_data$x[l];v=Z_data$y[l]
  ave=(n/(n+p))*mean(samples)+(p/(p+n))*Average
  var=v*(p+n)
  A=a+n/2;
  B=b+(1/2)*sum((samples-mean(samples))^2)+(n*p*(mean(samples)-Average)^2)/(2*(p+n))
  B=1/B
  normal=(1/sqrt(2*pi*var))*exp(-((z-ave)^2)/(2*var))
  gamma=((B^(-1))/(gamma(A+1)))*prod(rep(v/B,A-1))*exp(-v/B)
  Z_data$z[l]=normal*gamma
  Z_data[is.na(Z_data)]=0
}
z_data=rbind(z_data,Z_data)
}
library(rgl)

plot3d(z_data$x,z_data$y,z_data$z,type="l",xlab="average",ylab="precision",zlab="p(u,t|X,u0,p,a,b)",xlim=c(min(x),max(x)),ylim=c(min(y),max(y)),zlim=c(min(z_data$x),max(z_data$z)))

max_data=z_data %>% filter(z==max(z))
```

```{r}
#y=t,z=p(t|X,u0,p,a,b)を軸としてプロット

z_data=data.frame(y=0,z=0)
width=0.02;sep=0.001


y=seq(Precision-width,Precision+width,by=sep)
Z_data=data.frame(y=y,z=0)

for(l in 1:nrow(Z_data)){
   p=1;n=length(samples);z=1;v=Z_data$y[l]
  A=a+n/2;
  B=b+(1/2)*sum((samples-mean(samples))^2)+(n*p*(mean(samples)-Average)^2)/(2*(p+n))
  B=1/B
gamma=((B^(-1))/(gamma(A+1)))*prod(rep(v/B,A-1))*exp(-v/B)
  Z_data$z[l]=gamma
  Z_data[is.na(Z_data)]=0
}
z_data=rbind(z_data,Z_data)

plot(z_data$y,z_data$z,xlab="precision",ylab="p(t|X,u0,p,a,b)",xlim=c(min(y),max(y)),ylim=c(min(z_data$z),max(z_data$z)),type="p")
```
```{r}

#p(x_n+1|X,u0,p,a,b)=st(x|u_n,a_n,(1+(1/(n+p)))*b_n)を計算--式(3.39)

st_dis=function(X){
  n=length(samples)
  u_n=(n/(n+p))*mean(samples)+(p/(p+n))*u0
  a_n=2*a+n
  b_n=2*b+sum((samples-mean(samples))^2)+((n*p)/(n+p))*(mean(samples)-u0)^2
  u_n=u_n;v_n=a_n;b_n=b_n*(1+(1/(n+p)))
  
  y<-(gamma((v_n/2)+1/2)/gamma(v_n/2))*sqrt(b_n/(pi*v_n))*(1+(b_n*((X-u_n)^2))/v_n)^(-(v_n+1)/2)
  return(y)
  
  
}

width=0.5;sep=0.01
Z_data=data.frame(x=seq(Average-width,Average+width,by=sep),z=0)


for(m in 1:nrow(Z_data)){
  Z_data$z[m]=st_dis(Z_data$x[m])
}


plot(Z_data$x,Z_data$z,xlim=c(min(Z_data$x),max(Z_data$x)),ylim=c(min(Z_data$z),max(Z_data$z)),type="o",xlab="次回の予測値;x",ylab="student t分布の確率密度")


```


```{r}

# p.35 

dice_num=3

N=10

#サンプルの作成
z=sample(1:dice_num,N,replace=TRUE)

n=c()
for(j in 1:dice_num){
  
  n=c(n,length(z[z==j]))
}

dice_table=data.frame(num=1:dice_num,count=n)

alpha=c()
for(j in 1:dice_num){
  alpha=c(alpha,sample(1:sum(n),1,replace=TRUE))
}

#p(z_{1:n} | alpha )を計算

prob=(gamma(sum(alpha))/gamma(sum(n)+sum(alpha)))*(prod(gamma(n+alpha)/gamma(alpha)))





```

```{r}

#ディリクレ分布の最尤推定量の計算

dice_num=5

#観測データ

N=3
samples=20


data=data.frame(num=0,p=t(rep(0,dice_num)))

for(j in 1:N){
Samples=sample(1:dice_num,samples,replace=TRUE)

data_sub=data.frame(num=j,p=t(rep(0,dice_num)))
for(l in 1:dice_num){
  data_sub[,l+1]=length(Samples[Samples==l])/length(Samples)
}
data=rbind(data,data_sub)
}

data=tail(data,nrow(data)-1)



digamma=function(z){
  n=100000;
  
  Y=0;
  for(k in 0:n){
 Y=Y+(1/(k+z))
}

  y<-log(n)-Y
  return(y)
}

trigamma=function(z){
  
  n=100000
  
  Y=0
  for(j in 0:n){
    Y=Y+1/((j+z)^2)
  }
  
  return(Y)
}



alpha=rep(1000,dice_num)

ite=100

newton=data.frame(ite=0,alpha=t(alpha))
for(j in 1:ite){
  
  
  g=N*(digamma(sum(alpha))-digamma(alpha))+apply(log(data[,(1+1):(ncol(data))]),2,mean)
  
  H=array(0,dim=c(dice_num,dice_num))

  for(k in 1:dice_num){
    for(l in 1:dice_num){
      H[k,l]=N*(trigamma(sum(alpha))-ifelse(k==l,trigamma(alpha[l]),0))
    }
  }
  
  alpha=alpha-solve(H)%*%g
  newton_sub=data.frame(ite=j, alpha=t(alpha))
  newton=rbind(newton,newton_sub)
  
  }






```




```{r}
#gamma分布の最尤推定量

library(dplyr)
library(MCMCpack)
library(pforeach)


#samples
N=1000

x=rgamma(N,shape=30,scale=40)

#x=rgamma(N,lambda=5)

#精度の分布 p(t|x_{1:N},u,a,b)--(3.28)

#gamma分布の最尤推定量はb=mean(x)/aより、aをiteration  methodで求める。

ave_x=mean(x)

sum_log_x=sum(log(x))

digamma=function(z){
  n=100000;y=0;count=0
  while(count<n)({
  y<-y+(1/(z+count))
  count=count+1
  })
  y<-(log(n)-y)
  return(y)
}

trigamma=function(z){
  count=0;y=0
  while(count<100000)({
  y<-y+1/((z+count)^2)
  count=count+1
  })
  
  return(y)
}
a<-0.5/(log(mean(x))-mean(log(x)))
b=mean(x)/a #初期値

log_p_data=data.frame(num=1:20,a=0,b=0,diff_log_p=0)
for(j in 1:nrow(log_p_data)){
  log_p_data$a[j]=a
  log_p_data$b[j]=mean(x)/a
  if(j>1){
  log_p_data$diff_log_p[j]=abs(sum(log(dgamma(x,shape=a,scale=b)))-sum(log(dgamma(x,shape=log_p_data$a[j-1],scale=log_p_data$b[j-1]))))
  }
    
  a=1/(1/a+(mean(log(x))-digamma(a)-log(mean(x))+log(a))/(((1/a)-trigamma(a))*(a^2)))
  
  
}
#a<-0.5/(log(mean(x))-mean(log(x)))
#b=mean(x)/a

print(paste0("gamma分布の最尤推定量はa=",a,",b=",b))

#gamma分布の出力

u=mean(x)

num=seq(0,10000,by=10)

t_prob_data=data.frame(num=num,values=0)


t_prob_data$values=dgamma(num,shape=a,scale=b)

a_h=a-length(x)/2;b_h=b-(1/2)*sum((x-u)^2)
plot(num,t_prob_data$values,type="o",xlab="num",ylab="values",main="最尤推定量で予測された分布",col=2)



```

```{r}

library(dplyr)

#p.39

p0=1;

N=100

sample=rnorm(N,mean=3,sd=5)

#最尤推定量

mu0=mean(sample)

var=var(sample)

tau=1/var

#gamma分布の平均値より、tauの値が一番出やすいように設定する。

a0=tau;b0=1

#p(mu,tau|mu0,p0)

prob_data1=data.frame(average=0,tau=0,prob=0)

separate=0.1

sep_tau=seq(0,1,0.01)

width=10

for(j in 1:length(sep_tau)){
  
  prob_data=data.frame(average=seq(mu0-width,mu0+width,1),tau=sep_tau[j],prob=0)
  
  for(k in 1:nrow(prob_data)){
    
    prob_data$prob[k]=dnorm(prob_data$average[k],mean=mu0,sd=sqrt(1/prob_data$tau[k]))*dgamma(prob_data$tau[k],shape=a0,scale=b0)
    
  }
  prob_data[is.na(prob_data)>0]=0
  
  prob_data1=rbind(prob_data,prob_data1)
  
  plot(prob_data$average,prob_data$prob,col=2,type="o",xlab="average",ylab="prob",main=paste0("p(mu,p0|mu0,p0)(tau;",sep_tau[j],")"))
}

prob_data1=tail(prob_data1,nrow(prob_data1)-1)

#p(mu,tau|x_{x;sample},mu0,p0,a0,b0)

prob_data2=data.frame(average=0,tau=0,prob=0)



mu_n=(N*mean(sample))/(N+p0)+(p0*mu0)/(N+p0)

p_n=p0+N

a_n=a0+N/2

b_n=b0+(1/2)*sum((sample-mean(sample))^2)+((N*p0)/(2*(N+p0)))*(mean(sample)-mu0)^2

width=2

sep_tau=seq(0,1,0.01)

for(j in 2:length(sep_tau)){
  
  prob_data=data.frame(average=seq(mu0-width,mu0+width,0.01),tau=sep_tau[j],prob=0)
  
  for(k in 1:nrow(prob_data)){
    
    prob_data$prob[k]=dnorm(prob_data$average[k],mean=mu_n,sd=sqrt(1/(p_n*prob_data$tau[k])))*dgamma(prob_data$tau[k],shape=a_n,scale=b_n)
    
  }
  prob_data[is.na(prob_data)>0]=0
  
  #prob_data[abs(prob_data)<10^(-200)]=0
  
  prob_data2=rbind(prob_data2,prob_data)
  
  plot(prob_data$average,prob_data$prob,col=2,type="o",xlab="average",ylab="prob",main=paste0("p(mu,tau|x_{x;sample},mu0,p0,a0,b0)(tau;",sep_tau[j],")"))
}

prob_data2=tail(prob_data2,nrow(prob_data2)-1)

#p(x_{x;sample}|mu0,p0,a0,b0);サンプルが発生する確率

print(paste0("p(x_{x;sample}|mu0,p0,a0,b0)=",sqrt(p0/p_n)*((b0^a0)/(b_n^a_n))*(gamma(a_n)/gamma(a0))*sqrt(2*pi)^(-N)))

```
```{r}

#CRP(multiple)


library(dplyr);library(rmvnorm)

dim=3

mu0=sample(1:10,dim);sigma0=t(t(mu0))%*%t(mu0);K=1

N=100;alpha=5

mu_vec=c(mu0)

n_vec=c(1)

samples=c(1)

for(i in 1:N){
  
val=sample(c(1:(K+1)),1,prob=c(n_vec/(i-1+alpha),alpha/(i-1+alpha)))  
  
samples=c(samples,val)

if(val>K){
  
K=K+1  

mu_vec=c(mu_vec,rnorm(1,mu0,sigma0))
  
}

Z_data=data.frame(num=samples) %>% group_by(num) %>% summarise(n=n())

n_vec=Z_data$n

print(paste0("number of clusters;",K))

print(n_vec)

}

mu_vec=rmvnorm(length(n_vec),mean=mu0,sigma=sigma0)

```

```{r}
#CRP

library(dplyr)

mu0=5;sigma0=9;K=1

N=100;alpha=5

mu_vec=c(mu0)

n_vec=c(1)

samples=c(1)

for(i in 1:N){
  
val=sample(c(1:(K+1)),1,prob=c(n_vec/(i-1+alpha),alpha/(i-1+alpha)))  
  
samples=c(samples,val)

if(val>K){
  
K=K+1  

mu_vec=c(mu_vec,rnorm(1,mu0,sigma0))
  
}

Z_data=data.frame(num=samples) %>% group_by(num) %>% summarise(n=n())

n_vec=Z_data$n

print(paste0("number of clusters;",K))

print(n_vec)

}


```