Rでgroup lasso の実験 一般化線形モデルの正則化

x1<-sample(1:100)
x2<-sample(1:100)
y<-2*x1+0.01*x2+rnorm(100)
p<-data.frame(y,x1,x2)
lm(y~x1+x2,data=p)

f<-function(r){
  sum((y-r[1]*x1-r[2]*x2)^2)
}

optim(c(1,3),f)

##Lasso
f2<-function(r){
  sum((y-r[1]*x1-r[2]*x2)^2)+20000*sum(abs(r))
}

optim(c(1,3),f2)$par

#Ridge
library(dummies)
x<-rbinom(100,size=1,prob=0.5)
y<-100*x+rnorm(100)
xx<-dummy(x)
p<-data.frame(y,x1)
lm(y~.,p)
x1<-xx[,1]
x2<-xx[,2]
f<-function(r){
  sum((y-r[1]*x1-r[2]*x-r[3])^2)+0000*sum(abs(r))
}

optim(c(1,3,1),f)$par

f2<-function(r){
  sum((y-r[1]*x1-r[2]*x2-r[3])^2)+100*sum(abs(r))
}
optim(c(1,3,1),f2)$par

f3<-function(r){
  sum((y-r[1]*x1-r[2]*x2-r[3])^2)+100*sum((r)^2)
}
optim(c(1,3,1),f3)$par

#group lasso
f4<-function(r){
  sum((y-r[1]*x1-r[2]*x2-r[3])^2)+100*((sum((r[1:2])^2))^0.5+(r[3]^2)^0.5)
}
optim(c(1,3,1),f4)$par

一般化線形モデル

正則化
―逸脱度/データ数＋λ×正則化項？
によると、
http://bayes.sigmath.es.osaka-u.ac.jp/ftanaka/T/modeling/2016chap12_2.pdf
対数尤度＋正則化項でいいっぽい。

逸脱度（二乗誤差みたいなやつ）＝-2×対数尤度比?

https://mumu.jpn.ph/forest/computer/2017/08/20/11019/

###GLM check

x<-rbinom(1000,size=1,prob=0.5)
mu<-exp(log(100)*x+log(3))

y<-rgamma(1000,shape=mu,scale =0.5)
xx<-dummy(x)
#p<-data.frame(y,x1)
x1<-xx[,1]
x2<-xx[,2]
#group lasso
f5<-function(r){
    ll<-exp(r[1]*x1+r[2]*x2+r[3])
    l<- -log(dgamma(y,shape=ll,scale=(r[4])))
    l[l==Inf]<-10^8
    mean(l)+0.01*((sum((r[1:2])^2))^0.5+(r[3]^2)^0.5+(r[4]^2)^0.5)
}

ans<-optim(c(1,3,1,1),f5)$par
ans
#optim(c(1,3,1),f5)$par
y_th<-rgamma(1000,shape=exp(ans[1]*x1+ans[2]*x2+ans[3]),scale=ans[4])
qqplot(y_th,y)
plot(density(y))
points(density(y_th),col=2)

###Genraization###

x1<-(rbinom(1000,size=1,prob=0.5))
x2<-(rbinom(1000,size=2,prob=0.5))
#mu<-exp(log(60)*x1-log(30)*x2+log(3))
mu<-exp(log(60)*x1-0*x2+log(3))


xx<-dummy.data.frame(data.frame(factor(x1),factor(x2)))

y<-rgamma(length(mu),shape=mu,scale =0.5)
#xx<-dummy(x)
g<-c(1,1,2,2,2)
gamma_group_lasso<-function(y,xx,g,lambda=0.01){
  colno<-length(xx[1,]) 
  ff<-function(r){
    r2<-r[1:((colno)+2)]
    ll<-exp(as.matrix(xx)%*%(r[1:colno])+r[colno+1])
    #colno<-length(xx[1,])  
    #ll<-exp(r[1]*x1+r[2]*x2+r[3])
    l<- -log(dgamma(y,shape=ll,scale=(r[colno+2])))
    l[l==Inf]<-10^8
    g2<-c(g,length(g),length(g)+1)
    mean(l)+lambda*mean(tapply(r2,g2,function(x){mean(x^2)^0.5}))
  }
  rr<-rep(1,colno+2)
  optim(rr,ff,method="CG") #共役勾配法のがうまくいくっぽい
}

gamma_group_lasso(y,xx,g,30)