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

R R データ解析 ロジスティック回帰
  • 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：男性１，女性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を自動的に捨てている。

  • 例２：制約を入れてみる．まず，切片に制約をいれてみる。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*2

出力結果

> 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と同じ結果が初期値に依存せずにえれた．

    • 例３　制約：「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*3

}

#決定係数のスタート値
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

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

    • 例４　平均０制約：「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*4
pp3<-sum(x3)/length(x3)
pp4<-sum(x4)/length(x4)

return*5

}

#決定係数のスタート値
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*6
pp3<-sum(x3)/length(x3)
pp4<-sum(x4)/length(x4)

return*7

}

#決定係数のスタート値
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: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"

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

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

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

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

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