Hatena::ブログ(Diary)

arupaka-_-arupakaの日記

2018-09-10

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###

2017-03-16

Rでパラメータ制約付き非線形最適化とパラメータ制約付の最小二乗法とロジスティック回帰(カテゴリカルデータ、質的データ)

Rsolnpパッケージをもちいる。

Rで非線形制約条件付非線形最適化を行う方法 by Rsolnp

http://d.hatena.ne.jp/teramonagi/20091217/1261048574を参考に

  • とりあえず最小二乗法でためす.

library(Rsolnp)

#答えと問題のデータを作成

#y=5*x1-8*x2-3*x3

a1<-5

a2<--8

a3<--3

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

x3<-rnorm(1000,mean=0,sd=1)


y<-a1*x1+a2*x2+a3*x3+1*rnorm(1000,mean=0,sd=1)


###最小二乗法の目的関数

eq1 <- function(x)

{

a1<-x[1]

a2<-x[2]

a3<-x[3]

return(sum*1^2))

}

#係数a1+a2=0の非線形制約

eq2 <- function(x)

{

return((x[1]+x[2])^1)

# return(x_[1]^4)

}

#初期パラメータ

x0 <- c(1,-1,1)

#最適化 a1+a2=0の制約で最適化. eqB=0 は制約式=0という意味

solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

#solution <- solnp(x0,fun = objectiveFunc)

#結果出力

print(solution)

a1<-5

a2<--8

a3<--3

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

x3<-rnorm(1000,mean=0,sd=1)


y<-a1*x1+a2*x2+a3*x3+1*rnorm(1000,mean=0,sd=1)


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

#対数尤度

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+a[4]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}


#決定係数のスタート値

x0 <- c(1,1,1,0)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution <- solnp(x0,fun = log_li)

#結果出力

print(solution$par)

出力

print(solution$par)

[1] 4.50399000 -6.65304329 -2.54410502 -0.01038315


一応、標準のglmと比較

d<-data.frame(z,x1,x2,x3)

glm(z~.,d,family=binomial(link="logit"))

#Call: glm(formula = z ~ ., family = binomial(link = "logit"), data = d)

#Coefficients:

#(Intercept) x1 x2 x3

# -0.01038 4.50399 -6.65304 -2.54410

#

#Degrees of Freedom: 999 Total (i.e. Null); 996 Residual

#Null Deviance: 1386

#Residual Deviance: 294 AIC: 302

#>

だいたい同じ結果.

  • カテゴリカルデータ(質的データ)のロジスティック回帰(パラメータ制約なし,不定形)
    • 不定形のため初期値により最適解がかわる。
      • x3とx4が線形従属している。x3:男性1,女性0, x4:女性1, 男性0

a1<-5

a2<--8

a3<--3

a4<-2

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

##Man flg

x3<-rbinom(1000,prob=0.1,size=1)

##Woman flg

x4<-1-x3

y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+x4*a[4]+a[5]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}

#初期値を変えて、不定性を確認

x0a <- c(1,1,1,1,0)

x0b<-c(3,4,2,2,1.764)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution1 <- solnp(pars=x0a,fun = log_li)

solution2 <- solnp(pars=x0b,fun = log_li)

#結果出力

print(solution1$par)

print(solution2$par)

#glmと比較

d<-data.frame(z,x1,x2,x3,x4)

glm(z~.,d,family=binomial(link="logit"))

    • 出力結果

> #結果出力

> print(solution1$par)

[1] 4.352484439 -7.153849857 -3.101372925 1.684972724 -0.002192511

> print(solution2$par)

[1] 4.352489 -7.153857 -1.877847 2.908487 -1.225704

>

> d<-data.frame(z,x1,x2,x3,x4)

> glm(z~.,d,family=binomial(link="logit"))

Call: glm(formula = z ~ ., family = binomial(link = "logit"), data = d)

Coefficients:

(Intercept) x1 x2 x3 x4

1.683 4.352 -7.154 -4.786 NA

Degrees of Freedom: 999 Total (i.e. Null); 996 Residual

Null Deviance: 1378

Residual Deviance: 312.5 AIC: 320.5

初期値によって解が違うことがわかる.glmあx4を自動的に捨てている。


  • 例2:制約を入れてみる.まず,切片に制約をいれてみる。gslの結果にあわせる

#####Logistic regression with Categolical data and constraints####

####Example 1: Fixed thresholds

a1<-5

a2<--8

a3<--3

a4<-2

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

##Man flg

x3<-rbinom(1000,prob=0.1,size=1)

##Woman flg

x4<-1-x3

y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+x4*a[4]+a[5]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}

constraint <- function(x)

{

return((x[5]-1.683)^1)

# return(x_[1]^4)

}

#決定係数のスタート値

x0a <- c(1,1,1,1,0)

x0b<-c(3,4,2,2,1.764)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution1 <- solnp(pars=x0a,fun = log_li,eqfun = constraint,eqB =0)

solution2 <- solnp(pars=x0b,fun = log_li,eqfun = constraint,eqB =0)

#結果出力

print(solution1$par)

print(solution2$par)

d<-data.frame(z,x1,x2,x3,x4)

glm(z~.,d,family=binomial(link="logit"))

出力結果

> print(solution1$par)

[1] 4.2505029 -6.8385070 -4.4021475 -0.0818387 1.6830000

> print(solution2$par)

[1] 4.25048803 -6.83826458 -4.40227303 -0.08189154 1.68300000

> d<-data.frame(z,x1,x2,x3,x4)

> glm(z~.,d,family=binomial(link="logit"))

Call: glm(formula = z ~ ., family = binomial(link = "logit"), data = d)

Coefficients:

(Intercept) x1 x2 x3 x4

1.601 4.250 -6.838 -4.320 NA

予想通りglmと同じ結果が初期値に依存せずにえれた.

    • 例3 制約:「a3=-a4」を入れてみる(男性=-女性).

a1<-5

a2<--8

a3<--3

a4<-2

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

##Man flg

x3<-rbinom(1000,prob=0.1,size=1)

##Woman flg

x4<-1-x3

y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+x4*a[4]+a[5]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}

constraint <- function(a)

{

return*2

}

#決定係数のスタート値

x0a <- c(1,1,1,1,0)

x0b<-c(3,4,2,2,1.764)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution1 <- solnp(pars=x0a,fun = log_li,eqfun = constraint,eqB =0)

solution2 <- solnp(pars=x0b,fun = log_li,eqfun = constraint,eqB =0)

#結果出力

print(solution1$par)

print(solution2$par)

#出力 does not depends on the innitail value

#> print(solution1$par)

#[1] 4.4599154 -6.9335231 -2.5641968 2.5641968 -0.9094139

#> print(solution2$par)

#[1] 4.4599725 -6.9335964 -2.5644191 2.5644191 -0.9096364

予想通り解が初期値に依存しない.

    • 例4 平均0制約:「a3*p3+a4*p4=0」を入れてみる(ランダムサンプルして平均をとると男女効果は0になる).
      • すべてにこの制約をいれれば、切片が平均になる。

a1<-5

a2<--8

a3<--3

a4<-2

x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

##Man flg

x3<-rbinom(1000,prob=0.1,size=1)

##Woman flg

x4<-1-x3

y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+x4*a[4]+a[5]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}

constraint <- function(a)

{

#return*3

pp3<-sum(x3)/length(x3)

pp4<-sum(x4)/length(x4)

return*4

}

#決定係数のスタート値

x0a <- c(1,1,1,1,0)

x0b<-c(3,4,2,2,1.764)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution1 <- solnp(pars=x0a,fun = log_li,eqfun = constraint,eqB =0)

solution2 <- solnp(pars=x0b,fun = log_li,eqfun = constraint,eqB =0)

#結果出力

print(solution1$par)

print(solution2$par)


出力

> #結果出力

> print(solution1$par)

[1] 4.3169926 -6.6144392 -3.4385721 0.3693704 1.1827882

> print(solution2$par)

[1] 4.3704808 -6.6988868 -3.4753002 0.3733157 1.2258578


  • 係数ax+ayの平均が0制約. 回帰の分散が係数の説明力より十分小さいとき,係数をy=1/(1+exp(-mu))

で変換したものがyの平均確率(p=(p1+p2+..+pN)/N)に対応する.

#####Logistic regression with Categolical data and constraints####

####Example 3: p3*a3+p4*a4=0 Mean zero constraionts

####The case of var[a1x1+a1x2] << mu0, p=1(1+exp(-mu0));p~mean(z)

a1<-5

a2<--8

a3<--3

a4<-2

a1<-0.0309115

a2<--0.07597877

a3<--0.1681380

a4<-0.137884

#a1<-0

#a2<-0

#a3<-0

#a4<-0

a5<-0.4594540


x1<-rnorm(1000,mean=0,sd=1)

x2<-rnorm(1000,mean=0,sd=1)

##Man flg

x3<-rbinom(1000,prob=0.1,size=1)

##Woman flg

x4<-1-x3

x1<-scale(x1)

x2<-scale(x2)


#y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)+a5

y<-a1*x1+a2*x2+a3*x3+a4*x4+1*rnorm(1000,mean=0,sd=1)+a5


prob<-1/(1+exp(-y))

z<-rbinom(length(prob),prob=prob,size=1)

log_li<-function(a){

y<-x1*a[1]+x2*a[2]+x3*a[3]+x4*a[4]+a[5]

p<-1/(1+exp(-y))

-sum(z*log(p)+(1-z)*log(1-p))

}

constraint <- function(a)

{

#return*5

pp3<-sum(x3)/length(x3)

pp4<-sum(x4)/length(x4)

return*6

}

#決定係数のスタート値

x0a <- c(1,1,1,1,0)

x0b<-c(3,4,2,2,1.764)

#最適化

#solution <- solnp(x0,fun = eq1,eqfun = eq2,eqB =0)

solution1 <- solnp(pars=x0a,fun = log_li,eqfun = constraint,eqB =0)

solution2 <- solnp(pars=x0b,fun = log_li,eqfun = constraint,eqB =0)

#結果出力

print(solution1$par)

print(solution2$par)

#> 1/(1+exp(-0.45194))

#[1] 0.6111004

#> mean(z)

#[1] 0.61

#>



http://documentation.statsoft.com/portals/0/formula%20guide/Logistic%20Regression%20Formula%20Guide.pdf

http://home.hiroshima-u.ac.jp/tkurita/lecture/statface/node13.html

回帰

http://ifs.nog.cc/gucchi24.hp.infoseek.co.jp/MRA2.htm

http://norimune.net/625

*1:y-(a1*x1+a2*x2+a3*x3

*2:a[3]+a[4]

*3:a[3]+a[4]

*4:a[3]*pp3+a[4]*pp4

*5:a[3]+a[4]

*6:a[3]*pp3+a[4]*pp4

2015-07-31

ggplot2で両対数グラフ その3. いくつかの系列をプロット

ggplot2で両対数グラフ その3. いくつかの系列をプロット

 d<-data.frame(rr_rd_mean_0[[1]],rr_rd_sd[[1]],rr_rd_mean_0[[42]],rr_rd_sd[[42]],rr_rd_mean_0[[79]],rr_rd_sd[[79]])
 colnames(d)<-c("x1","y1","x2","y2","x3","y3")
 g<-ggplot(data=d)+theme_bw(base_size=24)+scale_x_log10(labels=trans_format("log10",math_format(10^.x)))+scale_y_log10(labels=trans_format("log10",math_format(10^.x)))
 g<-g+geom_point(aes(x=x3,y=y3),size=I(1),colour=I(3),shape=I(3),alpha=0.5)
 g<-g+geom_point(aes(x=x2,y=y2),size=I(1),colour=I(2),shape=I(2),alpha=0.5)
 
 g<-g+geom_point(aes(x=x1,y=y1),size=I(1),colour=I(4),shape=I(1),alpha=0.5)
 
 g<-g+geom_line(aes(x=x1,y=sqrt(2*x1/len1[1]+0.05^2*x1^2)),size=1,linetype=2)
 g<-g+geom_line(aes(x=x1,y=sqrt(2*x1/len1[42]+0.05^2*x1^2)),colour=1,size=1,linetype=2)
 g<-g+geom_line(aes(x=x1,y=sqrt(2*x1/len1[79]+0.05^2*x1^2)),colour=1,size=1,linetype=2)
 g<-g+xlab("Mean")+ylab("S.D.")
 print(g)

2015-06-02

ggplot 対数グラフに補助線をつける例

ggplot 対数グラフに補助線をつける例

サイズ指定をIでやるのはポイント。。

p<-qplot(vf+vf/2,iqr1,log="xy",geom="point",shape=I(17))+theme_bw(base_size=32)+scale_x_log10(labels=trans_format("log10",math_format(10^.x)))+scale_y_log10(labels=trans_format("log10",math_format(10^.x)))
x<-exp(seq(log(600),log(10^7),length.out=300))
p<-p+geom_line(aes(x,0.16*(x/400)^-0.3),colour=I(2),linetype=I(2),size=2)
p<-p+xlab("Count 2009")+ylab("IQR of r; 2010-2009")
print(p)

2015-05-07

ggplot2のqplotで簡単にきれなグラフ

事前にlibrary(ggplot2)はしておく

x<-1:10
y<-2*x
qplot(x,y,geom="line",xlab="date",ylab="count")+theme_bw(base_size=30)

plot(x,y,type="b")は,geom=c("line","point")

qplot(cc[4:length(cc)],exp(m1),xlab="Mean of start",ylab="Mean of 6 year ratio",log="x",geom=c("line","point"))+theme_bw(base_size=30)

対数プロットに補助線

v<-qplot(sort(r0b),length(r0b):1/length(r0b),log="xy",type="l",xlab="Counts",ylab="CDF")+scale_x_log10(labels=trans_format("log10",math_format(10^.x)))+scale_y_log10(labels=trans_format("log10",math_format(10^.x)))+theme_bw(base_size=25)
v<-v+geom_line(aes(x=(1:100),y=(0.2*(1:100)^-1)),colour=2,size=2,linetype=2)
print(v)