ベイズ統計モデリング (統計ライブラリー) 　安藤知寛先生　朝倉書店

```{r}

#ベイズ線形回帰

library(dplyr)

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])

Y=as.numeric(data$y)

beta_MLE=solve(t(X)%*%X)%*%t(X)%*%Y

beta0=c(3,0.2)

A=diag(1,2)

A_hat=solve(t(X)%*%X+A)

beta_hat=solve(t(X)%*%X+A)%*%(t(X)%*%X%*%beta_MLE+A%*%beta0)

inv_gamma=function(lambda0,v0,x){
  
  return(((lambda0/2)^(v0/2))*exp(-lambda0/(2*x))*(x)^(-(v0/2+1))/gamma(v0/2))
  
  
}

#事前分布p(sigma)

lambda=3;v=2

inv_gamma_data=data.frame(x=seq(0.1,10,0.1),density=0)

for(j in 1:nrow(inv_gamma_data)){
  
inv_gamma_data$density[j]=inv_gamma(lambda,v,inv_gamma_data$x[j])  
  
}

plot(inv_gamma_data$x,inv_gamma_data$density,type="o")

v_hat=nrow(data)+v

lambda_hat=lambda+t(Y-X%*%beta_hat)%*%(Y-X%*%beta_hat)+t(beta0-beta_MLE)%*%solve(solve(t(X)%*%X)+solve(A))%*%(beta0-beta_MLE)

#尤度関数p(beta|sigma)

library(mvtnorm)

rmvnorm(10,mean=beta0,sigma=A)

#周辺事後分布p(beta|Y,X)

p_beta_Y_X=function(b){
  
return((1+t(b-beta_hat)%*%solve(A_hat*as.numeric(lambda_hat/v_hat))%*%(b-beta_hat)/v_hat)^(-(v_hat+ncol(X))/2))  
  
  
}


mu=X%*%beta_hat;sigma_star=as.numeric(lambda_hat/v_hat)*(diag(1,nrow(X))+X%*%solve(A_hat)%*%t(X))



```

```{r}

#ベイズ多項式回帰

library(dplyr)

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

X=as.matrix(data[,colnames(data) %in% c("x1")])

dim=3

XX=array(0,dim=c(nrow(X),dim))

for(j in 1:dim){
  
XX[,j]=X^j  
  
}

X=XX

Y=as.numeric(data$y)

beta_MLE=solve(t(X)%*%X)%*%t(X)%*%Y

beta0=c(3,0.2,0.007)

A=diag(1,dim)

A_hat=solve(t(X)%*%X+A)

beta_hat=solve(t(X)%*%X+A)%*%(t(X)%*%X%*%beta_MLE+A%*%beta0)

inv_gamma=function(lambda0,v0,x){
  
  return(((lambda0/2)^(v0/2))*exp(-lambda0/(2*x))*(x)^(-(v0/2+1))/gamma(v0/2))
  
  
}

#事前分布p(sigma)

lambda=3;v=2

inv_gamma_data=data.frame(x=seq(0.1,10,0.1),density=0)

for(j in 1:nrow(inv_gamma_data)){
  
inv_gamma_data$density[j]=inv_gamma(lambda,v,inv_gamma_data$x[j])  
  
}

plot(inv_gamma_data$x,inv_gamma_data$density,type="o")

v_hat=nrow(data)+v

lambda_hat=lambda+t(Y-X%*%beta_hat)%*%(Y-X%*%beta_hat)+t(beta0-beta_MLE)%*%solve(solve(t(X)%*%X)+solve(A))%*%(beta0-beta_MLE)

#尤度関数p(beta|sigma)

library(mvtnorm)

rmvnorm(10,mean=beta0,sigma=A)

#周辺事後分布p(beta|Y,X)

p_beta_Y_X=function(b){
  
return((1+t(b-beta_hat)%*%solve(A_hat*as.numeric(lambda_hat/v_hat))%*%(b-beta_hat)/v_hat)^(-(v_hat+ncol(X))/2))  
  
  
}


mu=X%*%beta_hat;sigma_star=as.numeric(lambda_hat/v_hat)*(diag(1,nrow(X))+X%*%solve(A_hat)%*%t(X))


```



```{r}

#自己相関を持つ回帰モデルのベイズ推定

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])

Y=as.numeric(data$y)

p=ncol(X);n=nrow(X)

epsilon=rnorm(n)

E=array(0,dim=c(n-p,p))

epsil=epsilon[(p+1):n]

for(j in 1:nrow(E)){
  
epsilon_sub=epsilon[j:(j+p-1)]  
  
E[j,]=epsilon_sub  

}


beta0=c(1,1)

A=diag(1,p)

ps=solve(t(E)%*%E)%*%t(E)%*%epsil

for(j in (p+1):n){
  
val=rep(0,p) 

val_y=0
  
for(l in 1:p){
  
val=val+X[(j-l),]*ps[l]  

val_y=val_y+ps[l]*Y[j-l]
  
}  

X[j,]=X[j,]-val
  
Y[j]=Y[j]-val_y
}

X_star=X[(p+1):n,];Y_star=Y[(p+1):n]

beta_star=solve(t(X_star)%*%X_star)%*%t(X_star)%*%Y_star

A_hat=solve(t(X_star)%*%X_star+A)

beta_hat=solve(t(X_star)%*%X_star+A)%*%(t(X_star)%*%X_star%*%beta_star+A%*%beta0)

V_star=solve(t(E)%*%E)

predict=(X_star)%*%beta_hat

cor(predict,Y_star)

res=data.frame(num=c(1:length(Y_star)),Y=Y_star,predict=predict)

plot(res$num,res$Y,col=2,type="o",xlab="num",ylab="values",xlim=c(1,nrow(res)),ylim=c(min(Y_star,predict),max(Y_star,predict)))

par(new=T)

plot(res$num,res$predict,col=3,type="o",xlab="num",ylab="values",xlim=c(1,nrow(res)),ylim=c(min(Y_star,predict),max(Y_star,predict)))


```
```{r}

library(dplyr);library(MCMCpack)

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=sample(1:10,20,replace = T),x4=sample(20:30,20,replace=T))

X=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

Y=as.numeric(data$y)

n0=5;s0=0.01;

B0=diag(100,ncol(X))

beta0=rep(0,ncol(X))

ite=10000

beta_box=array(0,dim=c(ite,ncol(X)))

confidence_interval=array(0,dim=c(ite,2))

for(i in 1:ite){


sigma=rinvgamma(1,n0/2,s0/2)

B_inv=t(t(X[1,]))%*%t(X[1,])/sigma

for(j in 2:nrow(X)){
  
B_inv=B_inv+ t(t(X[2,]))%*%t(X[2,])/sigma 
  
}

B_inv=B_inv+solve(B0)

beta=solve(B0)%*%beta0

for(j in 1:nrow(X)){
  
beta=beta+t(t(X[j,]))*Y[j]/sigma  
  
}

beta=solve(B_inv)%*%beta

beta_box[i,]=c(beta)


}

hist(beta_box[,1])

hist(beta_box[,2])

hist(beta_box[,3])

hist(beta_box[,4])



```

```{r}

#非線形多項ロジスティックモデル

library(dummies)

data(iris)

Species=iris$Species

X=as.matrix(iris[,!(colnames(iris)%in% c("Species"))])

v=0.1;lambda=0.001

#ガウスカーネルに変換

mu=apply(X,2,mean);sig=apply(X,2,var)

#for(j in 1:ncol(X)){

#X[,j]=exp(-(X[,j]-mu[j])^2/(2*v*sig[j]))

#}

X=cbind(rep(1,nrow(X)),X)

Y=dummy((Species))

class=ncol(Y)-1

ite=100000

eta=0.01

w=array(0.01,dim=c(class,ncol(X)))

for(l in 1:ite){

pi=cbind(exp(X%*%t(w)),rep(1,nrow(X)))

pi_sum=apply(pi,1,sum)

pi=pi/c(pi_sum)

dw=t(Y[,-ncol(Y)]-pi[,-ncol(pi)])%*%X-nrow(X)*lambda*w

loglik=sum(Y*log(ifelse(pi>0,pi,1)))

print(loglik)

#w=w+eta*dw

dw2=t(X)%*%diag(pi[,1])%*%(diag(pi[,1])-diag(rep(1,nrow(X))))%*%X-nrow(X)*diag(rep(1,ncol(X)))

for(j in 2:(ncol(pi)-1)){

dw_mat=t(X)%*%diag(pi[,1])%*%diag(pi[,j])%*%X

dw2=cbind(dw2,dw_mat)

}

for(i in 2:(ncol(pi)-1)){

dw3=t(X)%*%diag(pi[,i])%*%diag(pi[,1])%*%X

for(j in 2:(ncol(pi)-1)){

if(i==j){

dw_mat=t(X)%*%diag(pi[,i])%*%(diag(pi[,j])-diag(rep(1,nrow(X))))%*%X-nrow(X)*diag(rep(1,ncol(X)))

}else{

dw_mat=t(X)%*%diag(pi[,i])%*%diag(pi[,j])%*%X

}

dw3=cbind(dw3,dw_mat)

}

dw2=rbind(dw2,dw3)

}

dw_vec=solve(dw2)%*%c(t(dw))

w=w-eta*t(matrix(dw_vec,ncol=class))

}






```

```{r}

#ベイズ統計モデリング

#ロジスティック回帰　

library(dplyr)

#目的変量（故障数）
Y=c(2,0,0,0,0,0,1,0,1,2,0,0,1,0,1,0,0,0,1,0,0,0,0)

#確率
Y=Y/6

#特徴量行列（華氏）　
X=cbind(rep(1,length(Y)),c(53,73,70,67,78,75,57,68,70,75,79,69,58,76,70,70,81,76,63,68,72,68,66))

X=as.matrix(X)

alpha=rep(0.01,(ncol(X)))

ite=10000

eta=0.1

alpha_beta_data=array(0,dim=c(ite,ncol(X)))

for(j in 1:ite){

p=1/(1+exp(-X%*%alpha))

V=array(0,dim=c(nrow(X),nrow(X)))

diag(V)=(p*(1-p))

D=diag(1/c(dnorm(X%*%alpha)))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(Y-p)

alpha=alpha+eta*mat

alpha_beta_data[j,]=alpha

}

res=data.frame(predict=p,Y=Y)




#プロビット回帰　

library(dplyr)

#目的変量（故障数）
Y=c(2,0,0,0,0,0,1,0,1,2,0,0,1,0,1,0,0,0,1,0,0,0,0)

#確率
Y=Y/6

#特徴量行列（華氏）　
X=cbind(rep(1,length(Y)),c(53,73,70,67,78,75,57,68,70,75,79,69,58,76,70,70,81,76,63,68,72,68,66))

X=as.matrix(X)

alpha=rep(0.01,(ncol(X)))

ite=10000

eta=0.1

alpha_beta_data=array(0,dim=c(ite,ncol(X)))

for(j in 1:ite){

p=pnorm(X%*%alpha)

V=array(0,dim=c(nrow(X),nrow(X)))

diag(V)=(dnorm(X%*%alpha)^2)/(p*(1-p))

D=diag(1/c(dnorm(X%*%alpha)))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%V%*%D%*%(Y-p)

alpha=alpha+eta*mat

alpha_beta_data[j,]=alpha

}

res=data.frame(predict=p,Y=Y)



#補対数対数回帰　

library(dplyr)

#目的変量（故障数）
Y=c(2,0,0,0,0,0,1,0,1,2,0,0,1,0,1,0,0,0,1,0,0,0,0)

#確率
Y=Y/6

#特徴量行列（華氏）　
X=cbind(rep(1,length(Y)),c(53,73,70,67,78,75,57,68,70,75,79,69,58,76,70,70,81,76,63,68,72,68,66))

X=as.matrix(X)

alpha=rep(0.01,(ncol(X)))

ite=10000

eta=0.1

alpha_beta_data=array(0,dim=c(ite,ncol(X)))

for(j in 1:ite){

p=1-exp(-exp(X%*%alpha))

V=array(0,dim=c(nrow(X),nrow(X)))

diag(V)=((1-p)*log(1-p))^2/(p*(1-p))

D=diag(c(1/(-(1-p)*log(1-p))))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%V%*%D%*%(Y-p)

alpha=alpha+eta*mat

alpha_beta_data[j,]=alpha

}

res=data.frame(predict=p,Y=Y)


```