スパース推定法による統計モデリング　川野・松井・廣瀬先生　共立出版

```r

#p.15 最小角回帰

data("EuStockMarkets")
X=as.matrix(EuStockMarkets[,1:3])
Y=as.numeric(EuStockMarkets[,4])
n=nrow(X)


lam=1

beta=solve(t(X)%*%X)%*%t(X)%*%Y
ite=10^4
eta=10^(-0)

for(l in 1:ite){
for(j in 1:length(beta)){
rj=Y-X[,-j]%*%beta[-j]
beta[j]=beta[j]*(1-eta)+eta*(sum(rj*X[,j])-n*lam*ifelse(beta[j]>0,1,ifelse(beta[j]<0,-1,0)))/sum(X[,j]^2)
}
L=sum((t(X)%*%(Y-X%*%beta)-n*lam*ifelse(beta>0,1,ifelse(beta<0,-1,0)))^2)
print(L)
}

```

```r

#p.32 Elastic Net
#式(2.4)を準ニュートン法で解く

data("EuStockMarkets")
X=as.matrix(EuStockMarkets[,1:3])
Y=as.numeric(EuStockMarkets[,4])
X=cbind(rep(1,nrow(X)),X)
n=nrow(X)

#Broyden Quasi-Newton Algorithm p.213
a=1
lam=1

f=function(x){
S=ifelse(x>0,n*lam*a,ifelse(x<0,-n*lam*a,0))
vec=-t(X)%*%(Y-X%*%x)+n*lam*(1-a)*x+S
return(c(vec))
}

beta=solve(t(X)%*%X)%*%t(X)%*%Y
ite=1000
eta=0.1
H=diag(f(beta))

for(l in 1:ite){
if(max(abs(f(beta)))>10^(-4)){  
beta_pre=beta 
beta=beta-eta*H%*%f(beta)  
s=beta-beta_pre
y=f(beta)-f(beta_pre)
H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H
print((f(beta)))
L=sum((Y-X%*%beta)^2)+sum(a*abs(beta)+(1-a)*beta^2/2)/lam
#print(L)
}}  


```

```r

#p.32 Elastic Net
#(2.4)を計算した結果を用いてパラメータ最適化

data("EuStockMarkets")
X=as.matrix(EuStockMarkets[,1:3])
Y=as.numeric(EuStockMarkets[,4])
X=cbind(rep(1,nrow(X)),X)
n=nrow(X)

a=1
lam=10

beta=solve(t(X)%*%X)%*%t(X)%*%Y
ite=10^4
eta=10^(-0)

for(l in 1:ite){
for(j in 1:length(beta)){
rj=Y-X[,-j]%*%beta[-j]
beta[j]=beta[j]*(1-eta)+eta*(sum(rj*X[,j])-n*lam*a*ifelse(beta[j]>0,1,ifelse(beta[j]<0,-1,0)))/(sum(X[,j]^2)+n*lam*(1-a))
}
L=sum((Y-X%*%beta)^2)+sum(a*abs(beta)+(1-a)*beta^2/2)*lam
print(L)
}


```

```r

#p.49 クラスターlasso
#勾配降下法

data("EuStockMarkets")
X=as.matrix(EuStockMarkets[,1:3])
Y=as.numeric(EuStockMarkets[,4])
X=cbind(1,X)
n=nrow(X)

n=2000
X=cbind(rnorm(n,3,1),rnorm(n,18,4),rnorm(n,9,1),rnorm(n,3,1),rnorm(n,6,2))
X=cbind(1,X)
Y=X[,2]*2+X[,4]*3+9*(X[,5])+6*(X[,6])+rnorm(n,1,1)

lam1=0 #
lam2=0.8 #

d=function(x,value){
return(ifelse(x>0,value,ifelse(x<0,-value,0)))
}

f=function(x){
vec=-t(X)%*%(Y-X%*%x)
if(lam1>0){
vec[-1]=vec[-1]+d(x,n*lam1)[-1]
}
if(lam2>0){
for(i in 2:length(x)){
vec[i]=vec[i]+sum(d(x[i]-x[-i],n*lam2))  
}}
return(c(vec))
}

beta=solve(t(X)%*%X)%*%t(X)%*%Y
ite=3*10^4
eta=0.0000001

for(l in 1:ite){
beta=beta-eta*f(beta)
L=sum((Y-X%*%beta)^2)/(2*n)+sum(abs(beta[-1]))*lam1
for(i in 2:(length(beta)-1)){
L=L+sum(abs(beta[i]-beta[(i+1):length(beta)]))*lam2  
}
print(L)
}

```

```r

#p.49 クラスターlasso
#式(2.4)を準ニュートン法で解く

data("EuStockMarkets")
X=as.matrix(EuStockMarkets[,1:3])
Y=as.numeric(EuStockMarkets[,4])
n=nrow(X)

n=2000
X=cbind(rnorm(n,3,1),rnorm(n,18,4),rnorm(n,9,1),rnorm(n,3,1))
Y=X[,1]*2+X[,3]*3+rnorm(n,1,3)

#Broyden Quasi-Newton Algorithm p.213
lam1=10^(1)
lam2=10^(5)

d=function(x,value){
return(ifelse(x>0,value,ifelse(x<0,-value,0)))
}

f=function(x){
vec=-t(X)%*%(Y-X%*%x)+d(x,n/lam1)
for(i in 1:length(x)){
vec[i]=vec[i]+sum(d(x[i]-x[-i],n/lam2))  
}
return(c(vec))
}

beta=solve(t(X)%*%X)%*%t(X)%*%Y
ite=10^3
eta=0.05
H=diag(f(beta))

for(l in 1:ite){
if(max(abs(f(beta)))>10^(-7)){  
beta_pre=beta 
beta=beta-eta*H%*%f(beta)  
s=beta-beta_pre
y=f(beta)-f(beta_pre)
H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H
#print((f(beta)))
L=sum((Y-X%*%beta)^2)/(2*n)+sum(abs(beta))/lam1
for(i in 1:(length(beta)-1)){
L=L+sum(abs(beta[i]-beta[(i+1):length(beta)]))/lam2  
}
print(L)
}}  

```

```r

#p.70 一般化lassoに対するアルゴリズム

data("EuStockMarkets")

X=as.matrix(EuStockMarkets)
Y=X[,4]
X=cbind(rep(1,nrow(X)),X[,1:3])
D=diag(1,ncol(X))


n=nrow(X)
p=1
lam=0.5
beta=rep(0.01,ncol(X))
r=rep(0.01,ncol(X))
u=rep(10,ncol(X))
ite=1000

for(l in 1:ite){
beta=solve(t(X)%*%X+n*p*t(D)%*%D)%*%(t(X)%*%Y+n*p*t(D)%*%(r-u/p))  
r=sign(D%*%beta+u/p)*ifelse(abs(D%*%beta+u/p)-lam/p>0,abs(D%*%beta+u/p)-lam/p,0)
u=u+p*(D%*%beta-r)
L=sum((Y-X%*%beta)^2)/(2*n)+lam*sum(r^2)+sum(u*(D%*%beta-r))+p*sum((D%*%beta-r)^2)/2
print(L) 
}


```

```r

#p.85 logistic lasso
data("Seatbelts")
df=data.frame(Seatbelts)

X=df[,colnames(df) %in% c("DriversKilled","drivers","front","rear","kms","vankilled")]
X=as.matrix(cbind(1,X))
Y=df$law
n=nrow(X)

beta=rep(0.000001,ncol(X))
eta=1
ite=300
lam=1

for(l in 1:ite){
pix=1-1/(1+exp(X%*%beta)) 
pix[pix==1]=0.99999999999
pix[pix==0]=0.00000000001
wi=pix*(1-pix)  
zi=X%*%beta+(Y-pix)/wi
beta[1]=beta[1]*(1-eta)+eta*sum(wi*(zi-X[,-1]%*%beta[-1]))/sum(wi)
for(j in 2:length(beta)){
zij=X[,-j]%*%beta[-j]  
value1=sum(wi*(zi-zij)*X[,j]) 
value2=n*lam
beta[j]=beta[j]*(1-eta)+eta*sign(value1)*(abs(value1)-value2)/sum(wi*X[,j]^2)
}
L=sum(Y*log(pix)+(1-Y)*log(1-pix))-lam*sum(abs(beta[-1]))
print(L)
}


```

```r

#p.85 poisson lasso
data("Seatbelts")
df=data.frame(Seatbelts)

X=df[,colnames(df) %in% c("PetrolPrice","kms","drivers")]
X=as.matrix(cbind(1,X))
Y=df$DriversKilled
n=nrow(X)

beta=rep(0.00005,ncol(X))
eta=0.1
ite=1000
lam=0.1

for(l in 1:ite){
wi=exp(X%*%beta) 
zi=X%*%beta+(Y-wi)/wi
beta[1]=beta[1]*(1-eta)+eta*sum(wi*(zi-X[,-1]%*%beta[-1]))/sum(wi)
for(j in 2:length(beta)){
zij=X[,-j]%*%beta[-j]  
value1=sum(wi*(zi-zij)*X[,j]) 
value2=n*lam
beta[j]=beta[j]*(1-eta)+eta*sign(value1)*(abs(value1)-value2)/sum(wi*X[,j]^2)
}
L=sum(Y*log(wi)-wi)/n-lam*sum(abs(beta[-1]))
print(L)
}


```

```r

#p.94 multinominal logistic lasso　
data(iris)

sets=unique(iris$Species)
Y=array(0,dim=c(nrow(iris),length(sets)))
X=as.matrix(cbind(rep(1,nrow(iris)),iris[,-5]))
for(i in 1:nrow(Y)){
no=c(1:length(sets))[sets==iris$Species[i]]
Y[i,no]=1  
}

beta0=rep(0.01,length(sets))
beta=array(0.01,dim=c(length(sets),ncol(X)-1))
n=nrow(iris)

ite=10^4
eta=0.1
lam=0.01

for(l in 1:ite){

b=t(cbind(beta0,beta))
values=exp(X%*%b)
pij=values/apply(values,1,sum)
wij=pij*(1-pij)
Zij=X%*%b+(Y-pij)/wij


beta0=beta0*(1-eta)+eta*apply(wij*(Zij-X[,-1]%*%t(beta)),2,sum)/apply(wij,2,sum)
beta0=beta0-mean(beta0)

for(k in 2:(1+ncol(beta))){
  Zij_childer=t(t(X[,-c(1,k)]%*%t(beta[,-k+1]))+beta0)
  Sval=t(wij*(Zij-Zij_childer))%*%X[,k]
  beta[,(k-1)]=beta[,(k-1)]*(1-eta)+eta*sign(Sval)*ifelse(abs(Sval)-n*lam>0,abs(Sval)-n*lam,0)/(t(wij)%*%(X[,k]^2))
}
beta=t(t(beta)-apply(beta,2,mean))

print(sum(Y*log(pij)+(1-Y)*log(1-pij)))
}

```

```r

# p.99 コックス回帰モデル lasso
#座標降下法　確認中

data=data.frame(id=c(1:65),t=c(1,1,2,2,2,3,5,5,6,6,6,6,7,7,7,9,11,11,11,11,11,13,14,15,16,16,17,17,18,19,19,24,25,26,32,35,37,41,42,51,52,54,58,66,67,88,89,92,4,4,7,7,8,12,11,12,13,16,19,19,28,41,53,57,77),X2=c(2.218,1.94,1.519,1.748,1.301,1.544,2.236,1.681,1.362,2.114,1.114,1.415,1.978,1.041,1.176,1.724,1.114,1.23,1.301,1.508,1.079,0.778,1.398,1.602,1.342,1.322,1.23,1.591,1.447,1.079,1.255,1.301,1,1.23,1.322,1.114,1.602,1,1.146,1.568,1,1.255,1.204,1.447,1.322,1.176,1.322,1.431,1.945,1.924,1.114,1.532,1.079,1.146,1.613,1.398,1.663,1.146,1.322,1.322,1.23,1.756,1.114,1.255,1.079),X3=c(9.4,12,9.8,11.3,5.1,6.7,10.1,6.5,9,10.2,9.7,10.4,9.5,5.1,11.4,8.2,14,12,13.2,7.5,9.6,5.5,14.6,10.6,9,8.8,10,11.2,7.5,14.4,7.5,14.6,12.4,11.2,10.6,7,11,10.2,5,7.7,10.1,9,12.1,6.6,12.8,10.6,14,11,10.2,10,12.4,10.2,9.9,11.6,14,8.8,4.9,13,13,10.8,7.3,12.8,12,12.5,14),X4=c(67,38,81,75,57,46,50,74,77,70,60,67,48,61,53,55,61,43,65,70,51,60,66,70,48,62,53,68,65,51,60,56,67,49,46,48,63,69,70,74,60,49,42,59,52,47,63,58,59,49,48,81,57,46,60,66,71,55,59,69,82,72,66,66,60),X5=c(0,0,0,0,0,1,1,0,0,1,0,1,0,1,1,0,0,0,0,0,1,1,0,0,0,1,0,0,1,0,1,1,0,1,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0),X6=c(10,18,15,12,9,10,9,9,8,8,10,8,10,10,13,12,10,9,10,12,9,10,10,11,10,10,9,10,8,15,9,9,10,11,9,10,9,10,9,13,10,10,10,9,10,9,9,11,10,13,10,11,8,7,9,9,9,9,10,10,9,9,11,11,12),delta=c(rep(1,48),rep(0,17)))

t=data$t;n=length(t)
delta=data$delta

X=as.matrix(data[,colnames(data) %in% c("X2","X3","X4","X5","X6")])
t_values=sort(unique(t))

# パラメータの初期値
ite=10^4
beta=rep(0.05,ncol(X))
lam=10^(-2)

for(l in 1:ite){
w=c()
z=c()
for(s in 1:length(t)){
if(t[s]==min(t)){
w=c(w,0)
z=c(z,sum(X[s,]*beta))
}else{  
Cj=t[t<t[s]]
val1=exp(sum(X[s,]*beta))
val2=c()
for(k in 1:length(Cj)){
Rj=t[t>=Cj[k]]  
val2=c(val2,sum(exp(X[t>=Cj[k],]%*%beta)))  
}
w_val=sum((val1*val2-val1^2)/(val2^2))
w=c(w,w_val)
z=c(z,sum(X[s,]*beta)+(delta[s]-sum(val1/val2))/w_val)    
}}  

#dbeta=0
for(j in 1:length(beta)){
zij=X[,-j]%*%beta[-j]
#beta[j]=beta[j]+sign(sum(w*X[,j]*(z-zij)))*ifelse(abs(sum(w*X[,j]*(z-zij)))-n*lam>0,abs(sum(w*X[,j]*(z-zij)))-n*lam,0)/sum(w*X[,j]^2)
beta[j]=(sum(w*X[,j]*(z-zij))-(n*lam/2)*ifelse(beta[j]>0,1,ifelse(beta[j]<0,-1,0)))/sum(w*X[,j]^2)
#dbeta=c(dbeta,abs(-2*sum(w*X[,j]*(z-X%*%beta))/n+lam*ifelse(beta[j]>0,1,ifelse(beta[j]<0,-1,0))))
}

L=-lam*sum(abs(beta))
for(s in 1:(length(t))){
L=L+2*(sum(X[s,]*beta)-log(sum(exp(X[t>=t[s],]%*%beta))))/n
}

#L=sum((w*(z-X%*%beta))^2)/n+lam*sum(abs(beta))
print(L)
}

```

```r
# p.99 コックス回帰モデル lasso
#Newton-raphson

data=data.frame(id=c(1:65),t=c(1,1,2,2,2,3,5,5,6,6,6,6,7,7,7,9,11,11,11,11,11,13,14,15,16,16,17,17,18,19,19,24,25,26,32,35,37,41,42,51,52,54,58,66,67,88,89,92,4,4,7,7,8,12,11,12,13,16,19,19,28,41,53,57,77),X2=c(2.218,1.94,1.519,1.748,1.301,1.544,2.236,1.681,1.362,2.114,1.114,1.415,1.978,1.041,1.176,1.724,1.114,1.23,1.301,1.508,1.079,0.778,1.398,1.602,1.342,1.322,1.23,1.591,1.447,1.079,1.255,1.301,1,1.23,1.322,1.114,1.602,1,1.146,1.568,1,1.255,1.204,1.447,1.322,1.176,1.322,1.431,1.945,1.924,1.114,1.532,1.079,1.146,1.613,1.398,1.663,1.146,1.322,1.322,1.23,1.756,1.114,1.255,1.079),X3=c(9.4,12,9.8,11.3,5.1,6.7,10.1,6.5,9,10.2,9.7,10.4,9.5,5.1,11.4,8.2,14,12,13.2,7.5,9.6,5.5,14.6,10.6,9,8.8,10,11.2,7.5,14.4,7.5,14.6,12.4,11.2,10.6,7,11,10.2,5,7.7,10.1,9,12.1,6.6,12.8,10.6,14,11,10.2,10,12.4,10.2,9.9,11.6,14,8.8,4.9,13,13,10.8,7.3,12.8,12,12.5,14),X4=c(67,38,81,75,57,46,50,74,77,70,60,67,48,61,53,55,61,43,65,70,51,60,66,70,48,62,53,68,65,51,60,56,67,49,46,48,63,69,70,74,60,49,42,59,52,47,63,58,59,49,48,81,57,46,60,66,71,55,59,69,82,72,66,66,60),X5=c(0,0,0,0,0,1,1,0,0,1,0,1,0,1,1,0,0,0,0,0,1,1,0,0,0,1,0,0,1,0,1,1,0,1,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0),X6=c(10,18,15,12,9,10,9,9,8,8,10,8,10,10,13,12,10,9,10,12,9,10,10,11,10,10,9,10,8,15,9,9,10,11,9,10,9,10,9,13,10,10,10,9,10,9,9,11,10,13,10,11,8,7,9,9,9,9,10,10,9,9,11,11,12),delta=c(rep(1,48),rep(0,17)))

t=data$t;n=length(t)
delta=data$delta

X=as.matrix(data[,colnames(data) %in% c("X2","X3","X4","X5","X6")])
t_values=sort(unique(t))

# パラメータの初期値
ite=10^4
eta=0.05
beta=rep(0.001,ncol(X))
lam=1

for(l in 1:ite){
dbeta=-ifelse(beta>0,lam,ifelse(beta<0,-lam,0))
ddbeta=array(0,dim=c(length(beta),length(beta)))
loglik=-lam*sum(abs(beta))
for(i in 1:nrow(X)){
if(sum(t>t[i])>0){  
x=X[t>=t[i],]  
dbeta=dbeta+delta[i]*(X[i,]-t(x)%*%exp(x%*%beta)/sum(exp(x%*%beta)))
H1=t(x)%*%diag(c(exp(x%*%beta)/sum(exp(x%*%beta))))%*%x
H2=t(x*c(exp(x%*%beta)))%*%diag(rep(1/sum(exp(x%*%beta)),nrow(x)))%*%(x*c(exp(x%*%beta)))
ddbeta=ddbeta-delta[i]*(H1-H2)
loglik=loglik+delta[i]*(sum(X[i,]*beta)-log(sum(exp(x%*%beta))))
}}  
beta=beta+eta*solve(ddbeta)%*%dbeta
print(loglik)
}

```

```r

#p.121 SPCA

data("EuStockMarkets")
X=as.matrix(EuStockMarkets)
S=cov(X)
K=ncol(X)
n=nrow(X)

#正則化パラメータ
lamvec=c(2,5,5,7)*10
lam=0.00001
ite=10^4

d=eigen(S)$valuesd
A=eigen(S)$vectors[,1:K]

for(l in 1:ite){
B=array(0,dim=dim(A))
for(i in 1:K){
Y=X%*%A[,i]
beta=solve(t(X)%*%X)%*%t(X)%*%Y
beta_pre=beta+1

while(max(abs(beta-beta_pre))>10^(-4)){
beta_pre=beta
for(j in 1:length(beta)){
rj=Y-X[,-j]%*%beta[-j]
beta[j]=(sum(rj*X[,j])-n*lamvec[i]*ifelse(beta[j]>0,1,ifelse(beta[j]<0,-1,0)))/(sum(X[,j]^2)+2*n*lam)
}
L=sum((Y-X%*%beta)^2)/2+lam*sum(beta^2)+lamvec[i]*sum(abs(beta))
#print(L)
#print(max(abs(beta_pre-beta)))
}
B[,i]=beta
}
A=(svd(t(X)%*%X%*%B)$u)%*%t(svd(t(X)%*%X%*%B)$v)
L=sum(diag(t(X-X%*%B%*%t(A))%*%(X-X%*%B%*%t(A))))+lam*sum(diag(t(B)%*%B))+sum(lamvec*apply(abs(B),2,sum))
print(L)
}


```