高次元データ分析の方法ーRによる統計的モデリングとモデル統合　朝倉書店　安藤知寛先生

```{r}

#罰則付き線形回帰のモデル統合　p.11

n=100
p=10
X=matrix(runif(n*p,-2,2),nrow=n,ncol=p)
b=c(2,6,-3,1,4,rep(0,5))
y=X%*%b+rnorm(n,0,1)

sign=c(1,1,1,2,2,2,3,3,3,3)


ite=100
eta=0.001
lam=1
beta=rep(1,ncol(X))
sigPL=1

for(l in 1:ite){
AIC_vec=c()
loglik_vec=c()
pre=array(0,dim=c(n,max(sign)))
for(j in 1:max(sign)){
beta[sign==j]=beta[sign==j]+eta*(t(X[,sign==j])%*%(y-X[,sign==j]%*%beta[sign==j])/sigPL-n*lam*beta[sign==j])
sigPL=sum((y-X[,sign==j]%*%beta[sign==j])^2)/n
loglik=-n*log(2*pi*sigPL)/2-sum((y-X[,sign==j]%*%beta[sign==j])^2)/(2*sigPL)-n*lam*sum(beta[sign==j]^2)/2
loglik_vec=c(loglik_vec,loglik)
pre[,j]=X[,sign==j]%*%beta[sign==j]
AIC_vec=-2*loglik_vec+2*length(sign[sign==j])
}
#print(loglik_vec)
print(AIC_vec)
}

W=exp(-AIC_vec/2)/sum(exp(-AIC_vec/2))
MUs=pre%*%W

```

```{r}

#p.11

data("airquality")

data=na.omit(airquality)

Y=data$Ozone

X=as.matrix(data[,colnames(data)!="Ozone"])

n=nrow(X)

ite=100

beta=rep(1,ncol(X))

lam=0.01

for(l in 1:ite){

sig_PL=sum((Y-X%*%beta)^2)/n  
    
beta=solve(t(X)%*%X+n*lam*sig_PL*diag(1,ncol(X)))%*%t(X)%*%Y  

loglik=-n*log(2*pi*sig_PL)/2-sum((Y-X%*%beta)^2)/(2*sig_PL)-n*lam*sum(beta^2)/2  
  
print(loglik)

}

```

```{r}

#p.17

data("anscombe")

R=as.matrix(anscombe[,c(1,4)])

X=as.matrix(cbind(rep(1,nrow(anscombe)),anscombe[,5:8]))

m=ncol(X)

n=nrow(X)
  
B=solve(t(X)%*%X)%*%t(X)%*%R 

sig=t(R-X%*%B)%*%(R-X%*%B)/n

loglik=-n*m*log(2*pi)/2-n*det(sig)/2-sum(diag(solve(sig)%*%(t(R-X%*%B)%*%(R-X%*%B))))/2



```

```{r}
#p.98 3.2.7 分位点ファクター回帰

#サンプルデータの作成
N=200
P=400
r=4
L=matrix(rnorm(P*r,0,1),nrow=P,ncol=r)
F=matrix(rnorm(N*r,0,1),nrow=N,ncol=r)
E=matrix(rnorm(N*P,0,3),nrow=N,ncol=P)
X=F%*%t(L)+E
Beta=c(-2,1,5,1)
z=F[,1]*Beta[1]+F[,2]*Beta[2]+F[,3]*Beta[3]+F[,4]*Beta[4]

VARe=1
y=z+rnorm(N,0,sqrt(VARe))

#パラメータの指定
tau=0.01
IC=10^10
k=4 #次元圧縮する次元数


#漸近主成分による次元圧縮
VEC=eigen(X%*%t(X)/(P*N))$vectors
factor=sqrt(N)*VEC[,1:k]
E=cbind(1,factor)#---次元圧縮された新たな座標行列
lam=t(factor)%*%X/N

#分位点回帰

#最小化したいコスト
f=function(B){
er=y-E%*%B
rho=er*(tau-ifelse(er<0,1,0))
return(rho)
}

#再急降下法による推定
b=rep(1,ncol(E))
ite=10000
eta=0.001
for(l in 1:ite){
error=y-E%*%b

db1=-(tau-1)*E
db2=-tau*E
sig_error=ifelse(error<0,1,-1)
if(sum(sig_error==1)>0 & sum(sig_error==-1)>0){
db1=db1[sig_error==1,]
db2=db2[sig_error==-1,]
db=rbind(db1,db2)
}else{
if(sum(sig_error==1)>0){
db=db1[sig_error==1,] 
}else{
db=db2[sig_error==-1,]  
}}

db=apply(db,2,sum)
b=b-eta*db
print(sum(f(b))) 

}

```

```{r}

#SCAD推定量

#サンプルデータの作成
N=200
P=400
r=4
L=matrix(rnorm(P*r,0,1),nrow=P,ncol=r)
F=matrix(rnorm(N*r,0,1),nrow=N,ncol=r)
E=matrix(rnorm(N*P,0,3),nrow=N,ncol=P)
X=F%*%t(L)+E
Beta=c(-2,1,5,1)
z=F[,1]*Beta[1]+F[,2]*Beta[2]+F[,3]*Beta[3]+F[,4]*Beta[4]

VARe=1
y=z+rnorm(N,0,sqrt(VARe))

#パラメータの指定
k=4 #次元圧縮する次元数


#漸近主成分による次元圧縮
VEC=eigen(X%*%t(X)/(P*N))$vectors
factor=sqrt(N)*VEC[,1:k]
E=cbind(1,factor)#---次元圧縮された新たな座標行列
lam=t(factor)%*%X/N


#再急降下法による推定

lambda=1;r=3.7

scad=function(B){
val1=lambda*abs(B)  
val2=(r*lambda*abs(B)-(B^2+lambda^2)*0.5)/(r-1)
val3=rep((r^2-1)*lambda^2/(2*(r-1)),length(B))  
val_sig=ifelse(abs(B)<=lambda,1,ifelse(abs(B)<=r*lambda,2,3))
values=matrix(c(val1,val2,val3),nrow=length(B))
val=c()
for(j in 1:length(B)){
val=c(val,values[j,val_sig[j]])  
}
return(val)  
}

b=rep(1,ncol(E))
b0=b
ite=10^3
eta=0.0001
cost_vec=c()
for(l in 1:ite){
  
if(l==1){  
db=-2*t(E)%*%(y-E%*%b)
}else{
db=-2*t(E)%*%(y-E%*%b)+2*b*(scad(b)-scad(b0))/((b-b0)^2)  
}  
b=b-eta*db
cost=sum((y-E%*%b)^2)+N*sum(scad(b))
cost_vec=c(cost_vec,cost)
print(cost)
}

plot(cost_vec)

```

```{r}

#高次元パネルデータ分析　p.103
#---ちょっと勉強中---

rm(list=ls())

#scad 推定量
lambda=1;r=3.7

scad=function(B){
val1=lambda*abs(B)  
val2=(r*lambda*abs(B)-(B^2+lambda^2)*0.5)/(r-1)
val3=rep((r^2-1)*lambda^2/(2*(r-1)),length(B))  
val_sig=ifelse(abs(B)<=lambda,1,ifelse(abs(B)<=r*lambda,2,3))
values=matrix(c(val1,val2,val3),nrow=length(B))
val=c()
for(j in 1:length(B)){
val=c(val,values[j,val_sig[j]])  
}
return(val)  
}

#サンプルデータ
data("iris")
X=as.matrix(iris[,1:3]);Y=iris[,4]
t=nrow(X)/3
class=as.integer(iris$Species)

#パラメータ

N=ncol(X)
W=matrix(c(Y),ncol=max(class))
beta=t(array(10,dim=c(max(class),N)))
eigens=eigen(W%*%t(W)/(N*t))  
dim=length(eigens$values[eigens$values>10^(-7)])
f=sqrt(t)*eigens$vectors
lam=t(f)%*%W/t
pre=f%*%lam
beta0=beta+1
cost=0;cost_pre=0
for(j in 1:max(class)){
cost=cost+sum((Y[class==j]-X[class==j,]%*%beta[,j]-pre[,j])^2)+t*sum(scad(beta[,j]))
cost_pre=cost_pre+sum((Y[class==j]-X[class==j,]%*%beta0[,j]-pre[,j])^2)+t*sum(scad(beta0[,j]))
}
pre_beta=beta0  
ite=300
eta=0.0001

for(l in 1:ite){
  
while(sum(abs(cost-cost_pre))>10^(-5))({ 
cost_pre=cost
cost=0
W=matrix(c(Y),ncol=max(class))
for(j in 1:max(class)){
if(sum(beta==beta0)==prod(dim(beta))){  
db=-2*t(X[class==j,])%*%(W[,j]-X[class==j,]%*%beta[,j])
}else{
db=-2*t(X[class==j,])%*%(W[,j]-pre[,j]-X[class==j,]%*%beta[,j])+2*(beta[,j])*(scad(beta[,j])-scad(beta0[,j]))/((beta[,j]-beta0[,j])^2)  
}  
beta[,j]=beta[,j]-eta*db
W[,j]=W[,j]-X[class==j,]%*%beta[,j]
cost=cost+sum((Y[class==j]-X[class==j,]%*%beta[,j]-pre[,j])^2)+t*sum(scad(beta[,j]))
}
print(cost)
})
cost_pre=cost
eigens=eigen(W%*%t(W)/(N*t))  
dim=length(eigens$values[eigens$values>10^(-7)])
f=sqrt(t)*eigens$vectors
lam=t(f)%*%W/t
pre=f%*%lam
cost=sum((W-f%*%lam)^2)
for(j in 1:max(class)){
cost=cost+sum((Y[class==j]-X[class==j,]%*%beta[,j]-pre[,j])^2)+t*sum(scad(beta[,j]))
}
print(cost)
}



```

```{r}

#p.107 モデル統合法

#情報量基準による線形回帰のモデル統合

n=100
p=10
X=matrix(runif(n*p,-2,2),nrow=n,ncol=p)
b=c(2,6,-3,1,4,rep(0,5))
y=X%*%b+rnorm(n,0,1)

#各モデルの推定
M<-3
sign=c(1,1,1,2,2,2,3,3,3,3)


AIC=c()
MU=matrix(0,nrow=n,ncol=M)

for(j in 1:M){

Xk=X[,sign==j]    
Bk=solve(t(Xk)%*%Xk)%*%t(Xk)%*%y
Sk=sqrt((sum((Xk%*%Bk-y)^2))/n)
AICk=n*log(2*pi*Sk^2)+n+2*(length(Bk)+1)
AIC=c(AIC,AICk)
MU[,j]=Xk%*%Bk

}

W=exp(-AIC/2)/sum(exp(-AIC/2))
MUS=MU%*%W

```

```{r}

#4.2.2 ロジスティック回帰モデルの統合　p.110

Y=c(2,rep(0,5),1,0,1,2,0,0,1,0,1,0,0,0,1,0,0,0,0)
X=cbind(rep(1,23),c(53,73,70,67,78,75,57,68,70,75,79,69,58,76,70,70,81,76,63,68,72,68,66))

prob=function(s,b){
if(s==1){
p=1-1/(1+exp(X%*%b))  
}else{
if(s==2){
p=pnorm(X%*%b)
}else{
p=1-exp(-exp(X%*%b))  
}}  
return(p) 
}

beta_box=matrix(0,nrow=3,ncol=ncol(X))
AIC_vec=c()
MU=array(0,dim=c(length(Y),3))
for(j in 1:3){
beta=rep(0.01,2)
ite=100
eta=0.01
for(l in 1:ite){
W=diag(c(prob(j,beta)*(1-prob(j,beta))))
guzai=Y-6*prob(j,beta)
beta=beta+eta*solve(t(X)%*%W%*%X)%*%t(X)%*%guzai
loglik=sum(Y*log(ifelse(prob(j,beta)>0,prob(j,beta),1)))+sum((6-Y)*log(ifelse(1-prob(j,beta)>0,1-prob(j,beta),1)))
AIC=-2*loglik+2*ncol(X)
print(loglik)
}  
AIC_vec=c(AIC_vec,AIC)
beta_box[j,]=beta  
MU[,j]=6*prob(j,beta)
}

w=exp(-AIC_vec/2)/sum(exp(-AIC_vec/2))
MUs=MU%*%w

```

```{r}

#p.123 ベイズアプローチによるモデル統合法

#線形回帰モデル

#データ発生
n=100
x=runif(n,-2,2)
b=c(3,2,-4,3,2,0,0,0)
X=cbind(1,x,x^2,x^3,x^4,x^5,x^6,x^7)
y=X%*%b+rnorm(n,0,1)

#各モデルの推定、BIC
M=7
ML=c() #各モデルの周辺尤度
MU=matrix(0,nrow=n,ncol=M)　#各モデルの条件付期待値


for(k in 1:M){
  
#事前パラメータ  
p=k+1
Ak=(10^(-3))*diag(1,p)
nu0=10^(-10)
lam0=10^(-10)

#パラメータ推定
Xk=X[,1:(k+1)]
Bkmle=solve(t(Xk)%*%Xk)%*%t(Xk)%*%y
Bk=solve(t(Xk)%*%Xk+solve(Ak))%*%t(Xk)%*%y
nu=nu0+n
lam=lam0+t(y-Xk%*%Bkmle)%*%(y-Xk%*%Bkmle)+t(Bkmle)%*%solve(solve(t(Xk)%*%Xk)+solve(Ak))%*%Bkmle
Sk=as.numeric((lam/2)/(nu/2+1))

#周辺尤度の計算
#(4.12)の分子
PriorB=-(p/2)*log(2*pi*Sk)+(1/2)*log(det(Ak))-t(Bk)%*%Ak%*%Bk/(2*Sk)
PriorS=(nu0/2)*log(lam0/2)-log(gamma(nu0/2))-(nu0/2+1)*log(Sk)-lam0/(2*Sk)
Dk=solve(t(Xk)%*%Xk+Ak)
#(4.12)の分母
PostB=-(p/2)*log(2*pi*Sk)+(1/2)*log(det(Dk))-t(Bk)%*%Dk%*%Bk/(2*Sk)
PostS=(nu/2)*log(lam/2)-log(gamma(nu/2))-(nu/2+1)*log(Sk)-lam/(2*Sk)
loglik=-(n/2)*log(2*pi*Sk)-as.numeric(t(y-Xk%*%Bk)%*%(y-Xk%*%Bk))

#すぐに0に近づく
ML=c(ML,(loglik+PriorB+PriorS)-(PostB+PostS))
MU[,k]=Xk%*%Bk
}

#重みベクトルの計算
PriorM=c(2:(M+1))/sum(c(2:(M+1)))
W=PriorM*ML/sum(PriorM*ML)

M=MU%*%W

```

```{r}



#ベイズアプローチによるモデル統合　p.129

library(dummies)

data(iris)

X=as.matrix(iris[,-ncol(iris)]);Y=dummy(iris[,ncol(iris)])

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

w=rep(1,(ncol(Y)-1)*ncol(X))

n=nrow(X)


ite=100000

eta=0.001

lambda=1

for(l in 1:ite){
  
w_mat=matrix(w,ncol=ncol(Y)-1)  
  
pis=exp(X%*%w_mat)/(1+apply(exp(X%*%w_mat),1,sum))

mat=Y[,-ncol(Y)]-pis

dw=w_mat

for(j in 1:(ncol(Y)-1)){
  
X2=X*c(mat[,j]);X2=apply(X2,2,sum)-n*lambda*w_mat[,j]  
  
dw[,j]=c(X2)  

}

dw2=array(0,dim=c(ncol(X)*(ncol(Y)-1),ncol(X)*(ncol(Y)-1)))

for(i in 1:(ncol(Y)-1)){
for(j in 1:(ncol(Y)-1)){ 
  
if(i==j){
  
dw2[((i-1)*ncol(X)+1):(i*ncol(X)),((i-1)*ncol(X)+1):(i*ncol(X))]=-t(X*c(pis[,i]))%*%(X*c(1-pis[,j]))-n*diag(lambda,nrow(w_mat))
  
}else{
  
dw2[((i-1)*ncol(X)+1):(i*ncol(X)),((j-1)*ncol(X)+1):(j*ncol(X))]=t(X*c(pis[,i]))%*%(X*c(pis[,j]))  
  
}

}}

w=w-eta*solve(dw2)%*%c(dw)

loglik=sum(Y[,-ncol(Y)]*log(ifelse(pis>0,pis,1)))-n*lambda*sum(w^2)/2

print(loglik)
  
}

Snbeta=-dw2/n

```