ベイズ信号処理  関原謙介先生　共立出版

```{r}

#線形正規モデルにおけるEMアルゴリズム

#p.39 勾配上昇法

data=data.frame(y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

Y=data$y

H=as.matrix(data[,!(colnames(data) %in% "y")])

lam=matrix(sample(1:9,nrow(H)^2,replace=T),ncol=nrow(H))+diag(rep(1,nrow(H)))

pthi=matrix(sample(1:9,ncol(H)^2,replace=T),ncol=ncol(H))+diag(rep(1,ncol(H)))

ite=10^4

eta=0.01

for(l in 1:ite){
  
pthi_pre=pthi;lam_pre=lam

gam=pthi+t(H)%*%lam%*%H
  
x_bar=solve(gam)%*%t(H)%*%lam%*%Y

pthi=pthi+eta*(solve(pthi)-(x_bar%*%t(x_bar)+solve(gam)))/2

lam=lam+eta*(solve(lam)-(Y-H%*%x_bar)%*%t(Y-H%*%x_bar)+H%*%solve(gam)%*%t(H))/2

loglik=-log(det(solve(lam)+H%*%solve(pthi)%*%t(H))+1)/2-t(c(Y))%*%solve(solve(lam)+H%*%solve(pthi)%*%t(H))%*%t(t(c(Y)))/2

print(as.numeric(loglik))

}



```

```{r}

#4.5 l2正則化ミニマムノルム解のベイズ的導出

# 勾配上昇法

data=data.frame(y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

Y=data$y

H=as.matrix(data[,!(colnames(data) %in% "y")])

alpha=beta=1

ite=2*10^5

eta=0.01

for(l in 1:ite){
  
gam=diag(rep(beta,ncol(H)))+beta*t(H)%*%H
  
x_bar=beta*solve(gam)%*%t(H)%*%Y

alpha=alpha+eta*(nrow(H)/(2*alpha)-(sum(x_bar^2)+sum(diag(solve(gam))))/2)

beta=beta+eta*(ncol(H)/(2*beta)-(sum((Y-H%*%x_bar)^2)+sum(diag(t(H)%*%H%*%solve(gam))))/2)

loglik=log(alpha)*nrow(H)/2-(sum(x_bar^2)+sum(diag(solve(gam))))*alpha/2+ncol(H)*log(beta)/2-beta*(sum((Y-H%*%x_bar)^2)+sum(diag(t(H)%*%H%*%solve(gam))))/2

#print(as.numeric(loglik))

loglik2=-det(diag(1/beta,nrow(H))+H%*%t(H)/alpha)/2-t(c(Y))%*%solve(diag(1/beta,nrow(H))+H%*%t(H)/alpha)%*%t(t(c(Y)))/2

print(as.numeric(loglik2))

}



```

```{r}

#p.40　ちょっと気になるなあ

#4.5 l2正則化ミニマムノルム解のベイズ的導出

data=data.frame(y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

Y=data$y

H=as.matrix(data[,!(colnames(data) %in% "y")])

alpha=beta=1

ite=10^4

for(l in 1:ite){
  
gam=diag(rep(beta,ncol(H)))+beta*t(H)%*%H
  
x_bar=beta*solve(gam)%*%t(H)%*%Y

alpha=sum(x_bar^2)+sum(diag(solve(gam)))/nrow(H)

beta=(sum((Y-H%*%x_bar)^2)+sum(diag(t(H)%*%H%*%solve(gam))))/ncol(H)

alpha=1/alpha;beta=1/beta

loglik=log(alpha)*nrow(H)/2-(sum(x_bar^2)+sum(diag(solve(gam))))*alpha/2+ncol(H)*log(beta)/2-beta*(sum(Y^2)-t(c(x_bar))%*%t(H)%*%Y-t(Y)%*%H%*%x_bar+t(c(x_bar))%*%t(H)%*%H%*%t(t(c(x_bar)))+sum(diag(t(H)%*%H%*%solve(gam))))/2

print(as.numeric(-loglik))

}



```

```{r}

#p.49

data=data.frame(y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

Y=data$y

H=(as.matrix(data[,!(colnames(data) %in% "y")]))

beta=10;alpha=rep(1,ncol(H))

ite=1000

for(l in 1:ite){
  
beta_pre=beta;alpha_pre=alpha


inv_pthi=diag(c(1/alpha))

gam=solve(inv_pthi)+beta*t(H)%*%H

x_bar=beta*solve(gam)%*%t(H)%*%Y

alpha=((x_bar^2)+(diag(solve(gam))))

alpha=1/alpha

beta=(sum(Y^2)-2*t(c(x_bar))%*%t(H)%*%Y+t(c(x_bar))%*%t(H)%*%H%*%t(t(c(x_bar))))/ncol(H)

beta=1/as.numeric(beta)

D=beta*(sum(Y^2)-2*t(c(x_bar))%*%t(H)%*%Y)+t(c(x_bar))%*%gam%*%t(t(c(x_bar)))

log_sigma_y=log(det(gam))-ncol(H)*log(beta)-log(det(solve(inv_pthi)))

sigma_y=diag(rep(1/beta,nrow(H)))+(H)%*%inv_pthi%*%t(H)

loglik=-(1/2)*(log_sigma_y+t(Y)%*%solve(sigma_y)%*%Y)
  
print(as.numeric(loglik))

}



```

```{r}



#p.55 5.4　ハイパーパラメータalphaの更新

data=data.frame(y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

Y=data$y

H=(as.matrix(data[,!(colnames(data) %in% "y")]))

beta=1;alpha=rep(1,ncol(H))

ite=1000

for(l in 1:ite){

beta_pre=beta;alpha_pre=alpha

inv_pthi=diag(c(1/alpha))

gam=solve(inv_pthi)+beta*t(H)%*%H

x_bar=beta*solve(gam)%*%t(H)%*%Y

alpha=diag(beta*solve(gam)%*%t(H)%*%H)/(x_bar^2)

sig_y=diag(1/beta,nrow(H))+H%*%inv_pthi%*%t(H)

f=as.numeric(log(det(sig_y))+t(c(Y))%*%solve(sig_y)%*%t(t(c(Y))))

print(f)

}




```

```{r}

#p.65 ベイズ因子分析

data(EuStockMarkets)

data=data.frame(EuStockMarkets)

K=nrow(data);M=ncol(data)

Y=data/1000

L=3

A=array(0.1,dim=c(M,L))

lambda=diag(rep(0.1,M))

ite=30

for(l in 1:ite){

gam=t(A)%*%lambda%*%A+diag(rep(1,L))

uk=solve(gam)%*%t(A)%*%lambda%*%as.matrix(t(Y))

Ruu=uk%*%t(uk)+K*solve(gam)

Ryu=t(as.matrix(Y))%*%t(uk)

A_hat=Ryu%*%solve(Ruu)

Ryy=t(as.matrix(Y))%*%as.matrix(Y)

lambda=diag(diag(Ryy-A%*%t(Ryu))/K)

sigma_y=solve(lambda)+A%*%t(A)

loglik=-K*log(det(sigma_y))/2

for(j in 1:K){
  
loglik=loglik-as.numeric(t(as.numeric(Y[j,]))%*%solve(sigma_y)%*%t(t(as.numeric(Y[j,]))))/2  
  
}

print(loglik)

}

```