Hatena::ブログ(Diary)

arupaka-_-arupakaの日記

2018-11-11

多変量回帰のチェック

多変量回帰で共通した情報を持つ場合は、

共通した情報程度の係数になる。 ただし、AIC Step wise などを使うとノイズによって

不安定になり、片方に集中することも。。 Ridge回帰なら両方の変数を選んでくる。

例えば、

a1=情報1+0.2情報2

a2=情報1+0.3情報3

で、y=情報1+情報2

なら、a1の回帰係数のが少し大きくなる。


もし、片方に回帰係数が寄った場合、ランダム誤差か、それとも、情報1の割合が

大きく、情報3の割合がちいさかったからか? yへの

情報もってるからといって、回帰係数がおおいとは限らない。

回帰係数が真に大きい潜在変数に近いほうが回帰係数が大きくなる。

b1が増えて、b2=0が増えたとき、a1とa2がどれくらい増えるか連立方程式で計算。それをもとに係数が決められる。

library(glmnet)
a1<-rexp(1000)
a2<-rexp(1000)
a3<-rexp(1000)
b1<-0.9*a1+0.05*a2
b2<-0.9*a1+0.06*a2

y<-3*a1+0.1*a2+0.0*a3+1*rnorm(1000,sd=1)

ans<-lm(y~b1+b2)
coef(ans)
step(ans)

ans0<-lm(y~1)

step(ans0,direction=c("forward"),scope=list(upper="y~b1+b2",lower="y~1"))


#Lasso
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 1)
coefficients(fit.glmnet.lasso.cv)

#Ridge
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 0)
coefficients(fit.glmnet.lasso.cv)


a1<-rexp(1000)
a2<-rexp(1000)
a3<-rexp(1000)
b1<-0.9*a1+0.03*a2
b2<-0.9*a1+0.03*a3

y<-3*a1+0.1*a2+0.1*a3+1*rnorm(1000,sd=1)

ans<-lm(y~b1+b2)
coef(ans)
step(ans)

ans0<-lm(y~1)

step(ans0,direction=c("forward"),scope=list(upper="y~b1+b2",lower="y~1"))


#Lasso
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 1)
coefficients(fit.glmnet.lasso.cv)

#Ridge
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 0)
coefficients(fit.glmnet.lasso.cv)




a1<-rexp(1000)
a2<-rexp(1000)
a3<-rexp(1000)
b1<-0.9*a1+0.5*a2
b2<-0.9*a1+0.03*a3

y<-3*a1+0.1*a2+0.1*a3+1*rnorm(1000,sd=1)

ans<-lm(y~b1+b2)
coef(ans)
step(ans)

ans0<-lm(y~1)

step(ans0,direction=c("forward"),scope=list(upper="y~b1+b2",lower="y~1"))


#Lasso
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 1)
coefficients(fit.glmnet.lasso.cv)

#Ridge
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 0)
coefficients(fit.glmnet.lasso.cv)

3 x 1 sparse Matrix of class "dgCMatrix"
                     1
(Intercept) -0.1048879
b1           0.7910052
b2           2.2927133

2016-06-29

Rでカテゴリカルデータ回帰の勉強(質的データ)

Rでカテゴリカルデータでの勉強(質的データ);数量化1類

(1)順番がつけられない質的データを考える。

   それを回帰する場合適切コーディングが求められる。

 *コーディングによって回帰係数がかわり、結果の解釈がかわる。

  普通に数字を割り振ると、順番に意味がでてしまうので望ましくない。

   例)東京1、大阪2、名古屋3:

それぞれにフラグを付けると順番の問題は防げるが、不定問題になるので定数項をいれると係数決められなくなる.

   例)(u1,u2,u3):東京(1,0,0) 大阪(0,1,0), 名古屋(0,0,1)

(v1,v2):男(1,0), 女(0,1)

    例えば、東京400円、大阪300円、名古屋200円とする、、

        モデルがy=a1*u1+a2*u2+a3*u3 ならば、a1=400,a2=300,a3=200となるが、

        モデルがy=a1*u1+a2*u3+a3*u3+C の場合ば、a1,a2,a3,Cの4変数で値は方程式3種類しかでないので不定問題になる(ノイズ0の極端場合で類推している)。   

          

(B)基本:一つのデフォルト値を決めて2進数でコーディング

   例1:(u1,u2):東京(1,0), u2大阪(0,1), 名古屋(0,0)

例2:(v1): 男(1), 女(0):

   モデルをとすして、y=a1*u1+a2*u2+b1*v1+C で回帰すれば不定にならない。

   解釈として、Cは、名古屋+女+定数 を合わせた値となり、a1,a2は名古屋を基準として有意かどうか、v1は女に対して有利がどうかになり、解釈が難しくなる。名古屋+女の絶対的水準も不明になる。

  

   例えば、東京0円追加、大阪0円追加、名古屋1000円追加だとしても、このモデルだと、東京大阪が-1000で有意となってしまう可能性がある。 つまり、回帰の意味が名古屋に対してどうかが基準になってくる。 名古屋よりやすいことは確認できるが、男女の係数と比較した時にどうなるか?

真の構造

  y=A1+A2+A3+V1+V2+CC とする。

     

(C)別の方法2 a1+a2+a3=0, b1+b2=0 と制限をいれる方法

->回帰係数のいみは、追加の値段の平均からのずれになる。これだと、上の例だと、-500,-500,500になり、東京大阪名古屋有意にはなる。東京大阪は真の構造では有意ではないが、やはり有意になるが、相対的に、名古屋が得というのはよくわかる。

(D)別の方法2、平均が0になるように、a1*ma1+a2*ma2+a3*ma3=0 b1*mb1+b2*mb2=0

ma1はカテゴリが東京のデータ個数、ma2はカテゴリが大阪のデータ数とすると、

Cは自動的に平均になる。 しかし、問題は、データ数に偏りがあるとき係数の正当性はいかに、

 データが多いほど係数の影響が小さくなる。

 

例えば、y-<y> を係数0で回帰するといいのかも、、2×3で値6個で方程式6個に見える、

しかし、さいごの2つ以外みちびけるから、4つになって不定になる。

なんていうか(D)と(C)は、東京大阪じゃなくて、東京に重みをおいた変数大阪に重みを置いた変数

つくるイメージ。男女の場合+1,-1 で、名古屋は計算でだす。

2013-08-28

 実験計画法と連立方程式

実験計画法と連立方程式

(1)実験計画法は以下の連立方程式をとくことに対応

y1=mu+a+b+ab+q[ab1]

y2=mu-a+b+ab+q[ab2]

y3=mu-a-b+ab+q[ab3]

未知数 mu[1つ],a[1つ], b[1つ], ab[1つ?],q[N-3つ] を解ければよい。

(2)効率のよい連立方程式のたてかたが実験計画法に対応。

   ほしい自由度以下(未知数の数)の実験数では情報をえられないのはこのため。

(3)自由度は着目変数方程式をとくための最低限必要な式(データ)の数(ノイズ=0を仮定したとき)。

たとえば、muだけなら、y1=mu でとける。 自由度1はn個のデータから1/n個づつ情報をもらっていることに相当。

つまり、平均をとるだけなら、n-1つの情報(ノイズになっているけど)をすててる。

すてる情報を以下に少なくバランスよくデータをとるかが実験計画法。

(データが多くなれば、多くなるほど、平均の精度、真の値との関係、が高くなることとの関係は? 

データを捨ててるわけではない。やはり,1個よりn個のがよい)

(4)平均の情報だけでは、N個のデータの情報は復号化できない。1個だけと同じ。

(5)しかし、正規分布を仮定の知見があれば、データ2個より、平均と分散の情報のほうがはるかにありがたい。

(6)ノイズが分布にしたがっているという知見があることが大切

(7)構造方程式仮定+正規分布仮定で、理論的には8個のデータから、

   未観測パタンも含め、2^7 パタンの情報をえることができる。

 A,B

AB: 平均でもAでもBでもない効果

ABC: 平均でもAでもBでもcでもABでもBCでもCAでも説明出来ないこうか

2013-02-21

単位根検定のお勉強

参考:

http://www.jaceksuda.com/teach/pse2011/class9_slides.pdf

単位根検定→定常性の検定→RW性の検定

70年代後半まで:

線形トレンドでlog(GDP)をフィッティング、確率部分はここからのゆらぎ

70年代問題

log(GDP)が線形増加がスローダウン→予見できなくなる。

多項回帰などより複雑なトレンドが考えられるようになる

→低周波数領域では有効

GNPのショックはランダムウォークに似ている。

→1982年 Nelson and Plosserが単位根検定

→単位根の存在を棄却できず。

実際には、単位根を検定する

ARMAを利用



とりあえず、AR(1)の例:

wikiより,

http://ja.wikipedia.org/wiki/%E8%87%AA%E5%B7%B1%E5%9B%9E%E5%B8%B0%E7%A7%BB%E5%8B%95%E5%B9%B3%E5%9D%87%E3%83%A2%E3%83%87%E3%83%AB

x(t+1)=φx(t)+c+ε

→t→∞で

x(t)=c/(1-φ)+Σφ^kε[t-k]

φ^n→0 だと定常へ

ARMAは力を過去から推定していることに相当。現在力のノイズ消し。

ラグ演算子

ε(t)=(1-Σφ_iL_i)X(t)=φx(t)

φ=1-Σφ_iL_i (自分-過去の影響の演算子)

AR:φX(t)=ε(t)

MA: X(t)=(1+Σ_{i=1 to q}θ_iL_i=θε_(t)

ここでθ=Σ_{i=1 to q}θ_iL_i←過去のショックと今のショックの合計(外力演算子)

したがってARMAは

φX(t)=θε(t)

→線形代数的に、

φ(z)=0 の単位根が1だと定常にならない。

→ということで単位根=1を検定

DF検定:場合わけがある

(1)単純なランダムウォークの場合

x(t+1)=φx(t)+ε

φを線形回帰で推定A

t_{φ}=(A-φ)/Se(A)を検定

もしφが0.9とか0.8とかなら、tはt分布

しかし、φ=1のときは解析的に単純にかけない分布に収束

→DF分布。DF分布はウィナー過程の積分的なもの。

→DF分布で検定。

(2)Constant+AR(1)

(y(t+1)-μ)=φ(y(t-1)-μ)+ε(t)

t^{μ}_{φ}というμに依存する検定量をつかう。

これはμ=1のときDF^{μ}分布に従う。

(3)ドリフト+定数+ランダムウォーク

y(t)-c-β(t)=φ(y(t-1)-c-β(t-1))+ε(t)

→t^{β}_{φ=1}統計量はDF^{β}分布に従う。

β=0のときは、DF^{μ}の分布のが強力な検定。

ADFtest(1984) DFtestの拡張版ARMAに対応

系列相関がある場合

Cov(u_i,u_j)≠0 ということ

Phillips-Perron Unit-Root Test:どんな系列相関に対応

検定量

Z(t)=(σ^2/λ^2)t_{ρ=0}-1/2*1(T・SE(ρ)/σ^2)

t_{ρ=0}=ρ/SE(ρ)

σ2^とλ^2の推定

σ^2=T^{-1}ΣE(u^2)

λ^2=ΣE(T^-1S_T^2)

S_T=Σu_t

u_tは誤差項、

Under the null hypothesis (帰無仮説 ρ=0 ←ん?、帰無仮説はρ=1?)では、tはADF分布に従う。

ほかのλの推定

λ^2=γ0+2Σ(1-j/(m+1))γ_j

γ_0=TΣu(t)^2

γ_j=1/TΣ_{t=j+1? to T}u_tu_{t-j}

http://faculty.washington.edu/ezivot/econ584/notes/unitroot.pdf

http://economia.unipv.it/pagp/pagine_personali/erossi/rossi_unit_roots_PhD.pdf

ssqrtl(λ^2)の計算

if (lshort)

l <- trunc(4 * (n/100)^0.25)

else l <- trunc(12 * (n/100)^0.25)

ssqrtl <- .C("R_pp_sum", as.vector(u, mode = "double"), as.integer(n),

as.integer(l), trm = as.double(ssqru), PACKAGE = "stats")

ssqrtl <- ssqrtl$trm

lshortはデフォルトがTらしい

.CはCの関数を実行。

R_pp_sum はおそらくこれ

http://fossies.org/dox/R-2.15.2/PPsum_8c_source.html

で、lが謎杉

実質やっていることは、nは回帰のデータ数、lは上で決まる定数

gamma2<-rep(0,l);for(j in 1:l){for(t in (j+1):n){gamma2[j]<-gamma2[j]+u[t-j]*u[t];}};gamma2<-gamma2/n

gamma2<-gamma0+2*sum([1-(1:l)/(1+l)]*gamma2)

gamma2

この部分 tableでいくつかのパタンで近似

tableは25,50,100,250,500,10^5の値が入ってる。

例えば、一つ目だと4.38,4.15,4.04,3.99,3.96

各パタンについて、それぞれnのときの値を計算.

approx(x,y,n,rule=2)

xとyの値を内挿して, nの値を求める、かつ, rule=2とする。

ルール2:外装状態のときの処理:最大値を使う (n=10^8 のときは3.96になる)

ルール1:計算しないNAにする

table <- cbind(c(4.38, 4.15, 4.04, 3.99, 3.98, 3.96), c(3.95,

3.8, 3.73, 3.69, 3.68, 3.66), c(3.6, 3.5, 3.45, 3.43,

3.42, 3.41), c(3.24, 3.18, 3.15, 3.13, 3.13, 3.12), c(1.14,

1.19, 1.22, 1.23, 1.24, 1.25), c(0.8, 0.87, 0.9, 0.92,

0.93, 0.94), c(0.5, 0.58, 0.62, 0.64, 0.65, 0.66), c(0.15,

0.24, 0.28, 0.31, 0.32, 0.33))

table <- -table

tablen <- dim(table)[2L]

tableT <- c(25, 50, 100, 250, 500, 1e+05)

tablep <- c(0.01, 0.025, 0.05, 0.1, 0.9, 0.95, 0.975, 0.99)

tableipl <- numeric(tablen)

for (i in (1L:tablen)) tableipl[i] <- approx(tableT, table[,

i], n, rule = 2)$y

*1:λ^2-σ^2)/(λ^2