2020-03-15

層で標準化した平均

ポエム

概要

周辺構造モデルで出てくる逆確率の最も簡単なパターンの計算練習

今回は特にcounter factualな話は出てこない。確率と期待値の機械的な計算に慣れる。

treatmentをA、共変量をLと読めばこれが何に通じるかがイメージできそう。 $I(\cdot)$ は指示関数。

基本的に $p(\cdot)$ は十分性質の良いものとする。(counter factual関係ないとは言ったけど、positivityが微妙に絡んでいる気はする。)

ゴール

多分次を導けること

$\displaystyle \tag{1.1} E \left( \frac{I_{A=a}Y}{p(A|L)} \right)=E _{L} \left\{ E_{Y|(A=a,L)} \left( Y| A=a,L \right) \right\}$

変形

積分についている ${S(\cdot)}$ はその台を示す。

(1.1)の左辺を積分に直す。

$\displaystyle \tag{1.2} E \left( \frac{I_{A=a}Y}{p(A|L)} \right)=\int _{L \in S(L)}dL \int _{A \in S(A)}dA \int _{Y \in S(Y)} \frac{I_{A=a}Y}{p(A|L)}p(Y,A,L) dY$

ここで、同時分布は次のように修正できることを思い出す。

$p(A|L)$ がipwとして扱われる量になっている。

$\displaystyle p(Y,A,L) = p(Y|A,L)p(A|L)p(L)$

すると、(1.2)の右辺は

$\displaystyle(1.2)=\int _{L \in S(L)}dL \int _{A \in S(A)}dA \int _{Y \in S(Y)} \frac{I_{A=a}Y}{p(A|L)}p(Y|A,L)p(A|L)p(L)dY$

で、p(A|L)が消えるから↓のようになって

$\displaystyle=\int _{L \in S(L)}dL \int _{A \in S(A)}dA \int _{Y \in S(Y)} I_{A=a}Yp(Y|A,L)p(L)dY$

p(L)とIはいったん考えないように追い出して

$\displaystyle=\int _{L \in S(L)}p(L)dL \int _{A \in S(A)} I_{A=a}dA \underbrace{\int _{Y \in S(Y)} Yp(Y|A,L)dY}_{E_{Y|(A,L)}(Y|A,L)}$

条件付き期待値 $E_{Y|(A,L)}(Y|A,L)$ が出てきた。

$E_{Y|(A,L)}$ の下付きの部分は、何の上で計算しているから自明じゃないから残しておく。

ということで、(1.2)はここまで計算できた

$\displaystyle \tag{1.3} E \left( \frac{I_{A=a}Y}{p(A|L)} \right)=\int _{L \in S(L)}p(L)dL \int _{A \in S(A)} I_{A=a}E_{Y|(A,L)}(Y|A,L)dA$

でAに関する積分だけど、 $\int_{A}f=\int_{S(A)}I_{A}f$ 的な計算を多分することになって、A=aの上でだけ残るから、次のように変形して良いはず。

$\displaystyle (1.3)=\int _{L \in S(L)} E_{Y|(A=a,L)}(Y|A=a,L)p(L)dL$

あとは普通に期待値として計算できるから

$\displaystyle =E_{L} \left\{ E_{Y|(A=a,L)}(Y|A=a,L) \right\}$

参考

Naimi, A. I., Cole, S. R., & Kennedy, E. H. (2017). An introduction to g methods. International journal of epidemiology, 46(2), 756–762. https://doi.org/10.1093/ije/dyw323 オンラインで読める

これのsupplementaryで周辺構造モデルのIPWを導いている。

第二部

で、結局何に使えるの？みたいなやつ
（ここからはまだうまく消化できていない）

supplementaryの　derivation of g-formulaの行間を埋める作業

time varyingなやつで、とりあえず2点のtreatmentがあるとする。exchangeabilityとかconsistencyは都合のいい感じに成立していると仮定する。

また、law of iterative expectationは勝手に用いる

ゴール

$E(Y^{a_{0},a_{1}})$ がなんかいい感じに変形できて、IPWの計算が役に立ちそうな雰囲気を出していることを確認する。

計算

まず、繰り返し期待値の性質から次のように変形できる

$\tag{2.1} E(Y^{a_{0},a_{1}})= E(\underbrace{E(Y^{a_{0},a_{1}}|A_{0})}_{\because =Y^{a_{0},a_{1}}})$

次に、繰り返し期待値の部分が交換可能であることを利用し、 $A_{0}$ に具体的な割り当てを考えてもよい

$\tag{2.2} (2.1)=E(E(Y^{a_{0},a_{1}}|A_{0}))=E(\underbrace{E(Y^{a_{0},a_{1}}|A_{0}=a_{0})}_{\because exchangeability})$

さらに、繰り返し期待値の性質から、Yをさらに条件付き期待値にする

$\tag{2.3} (2.2)=E(E(Y^{a_{0},a_{1}}|A_{0}=a_{0}))=E(E( \underbrace{E(Y^{a_{0},a_{1}}|A_{0}=a_{0},X,A_{1})}_{\because =Y^{a_{0},a_{1}}}|A_{0}=a_{0}))$

もう一度交換可能であることを利用し、 $A_{1}$ に具体的な割り当てを考えてもよい

$\tag{2.4} (2.3)=E(E( E(Y^{a_{0},a_{1}}|A_{0}=a_{0},X,A_{1})|A_{0}=a_{0})) \\ =E(E( \underbrace{E(Y^{a_{0},a_{1}}|A_{0}=a_{0},X,A_{1}=a_{1})}_{\because exchangeability}|A_{0}=a_{0}))$

$A_{0}=a_{0}$ と $A_{1}=a_{0}$ となっているので、 $Y^{a_{0},a_{1}}$ にconsistencyを考えることができる

$\tag{2.5} (2.4)=E(E(E(Y^{a_{0},a_{1}}|A_{0}=a_{0},X,A_{1}=a_{1})|A_{0}=a_{0})) \\ =E(E(\underbrace{E(Y|A_{0}=a_{0},X,A_{1}=a_{1})}_{\because consistency}|A_{0}=a_{0}))$

まとめると、次のように整理できたことになる。

$\tag{2.6} E(Y^{a_{0},a_{1}})=E(E(E(Y|A_{0}=a_{0},X,A_{1}=a_{1})|A_{0}=a_{0}))$

残念ながら、ここの条件付き期待値がなんの分布の上で計算しているのかきちんと整理できていないのだけど、 $E(E(Y|A_{0}=a_{0},X,A_{1}=a_{1})|A_{0}=a_{0})$ が第一部で計算した層で標準化した量に対応するはず。

というので、第一部で求めたのはなんだか色々役に立ちそうだ、と言うところらしい。

第三部

第一部で計算した量に対して、適当な条件が仮定できるときに潜在反応の期待値に変形できることを確認する。

ゴール

consistency

conditional exchangeability $p(Y^{a},A|L)=p(Y^{a}|L)p(A|L)$

positivity

のもとで次の関係になること

$\displaystyle \tag{3.1} E \left( \frac{I_{A=a}Y}{p(A|L)} \right)=E \left\{ Y^{a} \right\}$

計算

まず、consistencyが成立するなら(3-1)式左辺は、潜在反応に置き換えられる。

$\displaystyle E \left( \frac{I_{A=a}Y}{p(A|L)} \right)=E \left( \frac{I_{A=a}Y^{a}}{p(A|L)} \right)$

ついでに積分に戻す。

$\displaystyle =\int _{L \in S(L)}dL \int _{A \in S(A)}dA \int _{Y^{a} \in S(Y^{a})} \frac{I_{A=a}Y^{a}}{p(A|L)}p(Y,A,L)dY^{a}$

同時分布と条件付き分布の関係 $p(Y^{a},A,L)=p(Y^{a},A|L)p(L)$ を用いて次のように変形できる*1。

$\displaystyle =\int _{L \in S(L)}p(L)dL \int _{A \in S(A)}dA \int _{Y^{a} \in S(Y^{a})} \frac{I_{A=a}Y^{a}}{p(A|L)}p(Y^{a},A|L)dY^{a}$

続けて、conditional exchangeabilityから $p(Y^{a},A|L)=p(Y^{a}|L)p(A|L)$ と変形できて

$\displaystyle =\int _{L \in S(L)}p(L)dL \int _{A \in S(A)}dA \int _{Y^{a} \in S(Y^{a})} \frac{I_{A=a}Y^{a}}{p(A|L)}p(Y^{a}|L)p(A|L)dY^{a}$

$Y^{a}$ とAのそれぞれの積分に直すことができて

$\displaystyle =\int _{L \in S(L)}p(L)dL \underbrace{\int _{A \in S(A)} \frac{I_{A=a}}{p(A|L)} p(A|L) dA}_{=E (\frac{I_{A=a}}{p(a|L)} |L )} \underbrace{ \int _{Y^{a} \in S(Y^{a})} Y^{a} p(Y^{a}|L) dY^{a} }_{=E(Y^{a}|L)}$

で、難しそうな方の条件付き期待値が1になるから*2

$\displaystyle =\int _{L \in S(L)} \int _{Y^{a} \in S(Y^{a})} Y^{a} p(Y^{a}|L)p(L) dY^{a}dL$

あとは同時分布に直したら、これが $E[Y^{a}]$ になる。

$\displaystyle =\int _{L \in S(L)} \int _{Y^{a} \in S(Y^{a})} Y^{a} p(Y^{a},L) dY^{a}dL =E[Y^{a}]$

ということで(3.1)式が証明できた。

追記

スティルチェスでもなんでもない記法で展開しているから、 Aがdiscreteだったらそもそも記法的にアウトなのである。(そこだけΣで書き直せば解決するのだけど。)

で、 $A=a \in \mathbb{R}$ の状況下(一点)では $\int _{A \in S(A)} \frac{I_{A=a}}{p(A|L)} p(A|L) dA =0$ となる。

もしaが適当な長さをもつ区間であったなら、 $\int _{A \in a=(x,y)} \frac{I_{A=a}}{p(A|L)} p(A|L) dA=\int _{A \in a=(x,y)} I_{A=a} =y-x$ となり、1にならない。

となると、上記の式変形はAがdiscreteでないと成立しないのだろうか。困った。

*1:ついでにp(L)を左側に追いやっておく

*2:a={0,1}しかなく、しかも0か1かの1点なら確かに1なんだけど、 $A \subset \mathbb{R}$ で適当な部分集合をaにとると成立しない気がする

2020-03-14

潜在反応を交換する練習

ポエム

概要

よく忘れる内容の復習。どうしても覚えられない。

無視できそうな感じの独立性まで

条件付き分布が存在するかとかラドン=ニコディム微分がとかそう言うのではなく、ただ式変形の意味での復習

なお、causal inference bookとの整合的には次に注意

大文字と小文字は区別は面倒なので使い分けない（分けないと困る場合のみ使い分けようと思うがこの記事ではその必要はなかった）
割付はAであるがこの記事ではzを用いている

復習1：条件付き独立の定義

条件付き独立 $\displaystyle x \perp y |z$ はこう書けること

$\displaystyle p(x,y|z)=p(x|z)p(y|z)$

左辺を同時分布にして戻すとこうなる

$\displaystyle p(x,y,z)=p(x|z)p(y|z)p(z)$

復習2：交換可能性の定義

交換可能(exchangeable) $\displaystyle (y^{1},y^{0}) \perp z|x$ はこう書けること

$\tag{1} \displaystyle p(y^{1},y^{0},z|x)=p(y^{1},y^{0}|x)p(z|x)$

勝手に()で括られて一瞬「?」となるけど、xで条件づけると $\displaystyle (y^{1},y^{0})$ の同時分布とzが独立になっている、というのがポイント。

で、(1)の左辺を同時分布に戻すとこうなる。(p(x)を両辺にかけてある)

$\tag{2} \displaystyle p(y^{1},y^{0},z,x)=p(y^{1},y^{0}|x)p(z|x)p(x)$

復習3：交換可能性から期待値を計算する

同時分布については、 交換可能性には関係なく 、単純に次のように変形できる

$\tag{3} \displaystyle p(y^{1},y^{0},z,x)=\frac{p(y^{1},y^{0},z,x)}{p(y^{1},y^{0},x)}\frac{p(y^{1},y^{0},x)}{p(x)}p(x)=p(z|y^{1},y^{0},x)p(y^{1},y^{0}|x)p(x)$

(2)と(3)は同時分布 $p(y^{1},y^{0},z,x)$ を別の書き方で表したものに過ぎない。

つまり「(2)の右辺＝(3)の右辺」になる。でもよく見ると同じ項が多いので

$\require{cancel} \displaystyle p(z|y^{1},y^{0},x)\cancel{p(y^{1},y^{0}|x)}\cancel{p(x)}=\cancel{p(y^{1},y^{0}|x)}p(z|x)\cancel{p(x)}$

結局こうなる

$\tag{4} p(z|y^{1},y^{0},x)=p(z|x)$

でも、この左辺は次のように変形できる

$$ \begin{align} \displaystyle p(z|y^{1},y^{0},x) & =\frac{p(y^{1},y^{0},z,x)}{p(y^{1},y^{0},x)}\frac{p(z,x)}{p(z,x)}\frac{p(x)}{p(x)}\\ \displaystyle & =\frac{p(y^{1},y^{0},z,x)}{p(z,x)}\frac{p(x)}{p(y^{1},y^{0},x)}\frac{p(z,x)}{p(x)}\\ \displaystyle & =p(y^{1},y^{0}|z,x)\frac{1}{p(y^{1},y^{0}|x)}p(z|x) \end{align} $$

最初の＝は条件付き確率を戻したのと、1になるp()を追加している。次の＝は順番を入れ替えただけ、最後は条件付き確率に戻した。

これは(4)の左辺をただ変形したもので、p(z|x)に等しい。
結果として、次が成立する。

$\tag{5} p(y^{1},y^{0}|z,x)=p(y^{1},y^{0}|x)$

なので、たとえば $y^{0}$ を周辺化してから $y^{1}$ の期待値を計算すれば、次のようになる。

$E(y^{1}|z,x)=E(y^{1}|x)$

$y^{0}$ についても同様 $E(y^{0}|z,x)=E(y^{0}|x)$

これがmean的なexchangeable。

2020-03-08

A structual approach to selection bias

どくしょかんそうぶん

というのを読んだ。 2004年の論文らしい。

自分用のメモ

associationの指標がcausalな指標とならない場合の一つに「共通の結果でconditioningしてしまった、というのが挙げられる。例えば

case-controlで人選ミスってしまう
RCTでもなんかいろいろあってloss to follow upして（欠測データになる）しまう
（被験者をボランティアとかで募ってしまって）自己選択バイアス

で、まぁ選択バイアスみたいなのがあって、適当にconfounderで調整しようとしてもselectionに対策できていないと結局それはcausalな指標のestimatorとしてはbiasedだよねっていう。

どうやって調整するのかで、stratificationでダメな時があるから、IPWとかg-estimationとか使ったほうがいいよねー。

あと、このpaperのメインの話ではないけど、選択バイアスとか交絡とか言葉の意味が人によってまちまちだから、統一的に使えるようになった方がわかりやすくていいよね。

みたいなことが書いてあった。

感想

causal inference bookで読んだような話が書いてあった。（同じ著者なんだからそりゃそうだよね。）

なので、例の本の8章、14章を読み直すのが良さそうである。

あの本、time varyingのあたりから難しくてボコボコにされた。 22章のtarget trialについては読んですらいない。

きちんと読み直さないと。

2020-02-15

パーセンタイル的なあれを数値的に計算する2

ポエム

前回、回帰もできるし罰則もつけられるよねー

ってことで試した。

tau適当にいじってもそれらしい結果が出てきた。

全く検証してないけど。

多分、定数項には正則化項をつけないのが正しい気がする。

n<-50
a<-3
b<-2.5

x <- runif(n)
y <- a+ b*x + rnorm(n,mean=0,sd=b*x)

dmat <- cbind(rep(1,n),x)

#初期値をOLS解にする
lres<-lm(y~x,data=data.frame(x=x,y=y))
beta<-lres$coefficients

#%点
tau<- 0.9
#正則化項　15くらいまでにしないと発散して解が出ない
lambda<-0.01

#学習率とか計算の回数
learning_rate <- 0.1
itermax<-1000

#1000くらいじゃ収束しない時がある
for(i in 1:itermax){
  #ざんさ
  resd<- y - dmat %*% beta
  #一回微分
  temp<- tau-(resd<0)
  beta_new <- colSums(temp[,1] * (-dmat))/n + lambda * beta
  #目的関数(いらないけど)
  #objF<- sum((tau-(resd<0))*resd)/n
  
  #更新
  beta <- beta - learning_rate*(beta_new)
}

#あってるかどうか見る
plot(x,y,xlim=c(0,1),ylim=c(0,10))
par(new=T)
plot(x,dmat %*% beta,xlim=c(0,1),ylim=c(0,10))

2020-02-15

パーセンタイル的なあれを数値的に計算する

ポエム

なんとなく、パーセンタイルって数値的に計算できないのかなと思った。

二回微分計算しようと思ったけど、一回微分で符号関数出てくるし、もう一回微分したらデルタ関数になりそうだし、そしたら逆行列なくなりそうな気がしてやめた。

線形計画法で解いた方がいいらしいけど知らん。

#適当にデータを作る
yvec <- runif(100)
#何%点か
tau<- 0.4
#計算されたパーセンタイル
mu<- 0

learning_rate <- 0.1
# 長さ
n <- length(yvec)
for(i in 1:50){
  #残差
  resd<-yvec-mu
  #一回微分
  #わざわざ符号関数を使って描くとこうなるけど、別に下側の式で良い。
  #mu_new <- sum(abs(tau-(resd<0))*sign(resd))/n
  mu_new <- sum(tau-(resd<0))/n
  #目的関数(いらないけど)
  #objF<- sum((tau-(resd<0))*resd)/n
  #更新
  mu <- mu + learning_rate*(mu_new)
}

これで、resdの部分を線形推定量なんかにしたら、回帰にできるね。さらに言うと、とても性質の良い正則化項なら簡単に追加できる。

SNCDなる方法で、elastic netな方法もできるらしい。読んでみたところ、coordinate descentで、アルゴリズムもしっかり解説してくれてたので、頑張れば作れそうだった。

2020-01-03

統計検定各級の対策事例(1級、準1級、2級、3級)

テストたいさく

TL;DR

基本的に過去問に過学習すれば受かる。

級	自分の合格状況	指導経験	お気持ち
1級	合格	なし	実は数学力いらない
準1級	合格	あり	浅く広く適当で受かる
2級	合格	あり	この辺から対策が必要になる
3級	未受験	あり	社会人はノー勉では取れない
4級	未受験	なし	社会人なら取れる

注意点

全て個人の試験対策や、幾らか他人を指導した範囲での内容です。

高校生〜大学生にはあまり役に立たない内容で、目安として社会人（20〜40才程度）向け。

さらに、数学を真面目に学んでこなかった人向け。数学にこだわる人はこの記事でなく鍋谷(1978)とか読んでください。

1級

状態

自身が合格(数理・応用)

指導経験は無し

対策

1. 検索で出てくる大学の数理統計の講義資料

「数理統計学」なんかで検索すると、大学の講義レジュメがたくさん出てくる。
こういうのを使えば、大体勉強できてしまう。

個人的な指針として、数学科レベルの、測度論に踏み込む資料は必要ない点を強調したい。

過去問を見ればわかるが、 $\varepsilon - \delta$ を用いて証明する問題は出ない。もちろん位相も測度も関数解析も不要である*1。

ぶっちゃけた話、ヤコビアンを含む変数変換と部分積分と $\sum$ の変形さえできれば、高校数学に毛が生えたレベルで十分戦える。

この辺をなんとなく抑えてから資料を選ぶと良い。

2. 公式の過去問

これは当然なので解説しない。

3. RSSのGraduate Diplomaの過去問

RSS(王立統計協会)は、以前に認定試験を実施していた。Higher CertificateとGraduate Diplomaであり、後者の方が難しい。

この試験は過去の実施問題が公開されており、graduate diplomaのmodule 1とかmodule 2あたりが妙に統計数理の対策に使える。
1級過去問解いた後に見るとなんとなく既視感を覚える。 higher certificateも少し触れておくと計算練習になる。

個人的な感想

以下に挙げる書籍は、1級の参考書籍として挙げられることが多い。
しかし、個人の見聞の範囲からするとやりすぎな内容も含む *2 。ここまでやらなくても合格できる 。

久保川「現代数理統計学の基礎」(共立出版)
稲垣「数理統計学」(裳書房)
竹村「現代数理統計学」(創文社)

ちなみに私は次の書籍を使用した。上の3冊より難易度が低いため、無能の私にはこちらの方が良かった*3。

鈴木「数理統計学 -基礎から学ぶデータ解析-」(内田老鶴圃)

もちろんこれも丸々一冊マスターする必要はない。

準1級

状態

自身が合格

部分的に指導経験あり

対策

WEBで簡単そうな資料を探して手当たり次第読んだ。検索方法は、クラスター分析、主成分分析、カーネル法、など関連しそうな用語で出てきた物だけ。

数理的な理解は一切ない。

個人の感想

2級と一緒に受験したら合格した。
体感の目安として2級の問題を解いて安定して8割以上解けるようになるあたりから挑戦権がある、くらい。

なんとなく勉強して受かっただけなので、あまり偉そうなことは言えない。

2級

状態

自身が合格

指導経験あり

内容

ここに以前まとめた。

紙のテストなら6割ちょっと(23/35)得点できれば受かる。

ben-key.hatenablog.com

3級

状態

受けてない

指導経験のみあり

教えていた時の感想

3級のためだけの対策をやったことはない。

2級の指導が多かったため、その部分集合で基本的に受かる。気合の入った試験対策が必要かという観点で言うと、そんなに必死に対策するレベルではない。

なお内容が改定され、ちょっと試験範囲が増えるので、今後難易度は少し上がるはず。

社会人相手に試験対策を指導していた際、指導を始める前に内容確認として4級と3級の内容を解かせていたことがある。

結果として、 何も勉強していない社会人にいきなり解かせると、4級はほぼ間違いなく合格点取れるが、3級はほぼ間違い無く合格点取れない という感じだった。

ただし指導した社員は「ほぼ文系＆しかしそれなりにいい大学を出ている」人が多かったので、本当に一般人相手に解かせるともう少しレベルは落ちると思う。

*1:もちろん知っていてマイナスになることはない

*2:もちろん統計検定の対策専用書籍じゃないから当然である

*3:久保川の本が発売されるよりも前であり、そもそも選択肢に存在しない時期だったのもある。今だったら久保川がいいのだろうか？

2019-12-14

2019年の半径数メートル

にっきちょう

1月

causal inference bookのpart1を読んだ
確率微分方程式の定期勉強会に参加した
技術書店で当選したため、何か書くことになった

3月

機械学習のための特徴量エンジニアリング　を読んだ
ガウス過程と機械学習を途中まで読んだ

4月

技術書店で本を出した
会社が嫌になり、転職の面接を受けた
集合と位相そのまま使える回答の書き方を読んだ

5月

椅子、机、本棚を買った。
causal inference bookのpart2を読んだ
転職先から内定が出たが、悩んでいた

6月

転職を決めた

7月

ルベーグ積分の定期勉強会を完走した
確率微分方程式の定期勉強会に脱落した
入門機械学習による異常検知を読んだ
とても大切な人が亡くなった

8月

複素関数入門を途中まで読んだ
引っ越した

9月

転職した。仕事はpower pointとwordばかりで虚無
関数解析の定期勉強会に参加した

10月

予測にいかす統計モデリングの基本　を読んだ
点過程の時系列解析　を読んだ

11月

家族の幸せの経済学を読んだ
政策評価のための因果効果の見付け方を読んだ
入門統計的因果推論(pearl のprimerのやつ)　を読んだ
kaggle で勝つデータ分析の技術を読んだ

12月

社会科学のためのベイズ統計モデリングを読んだ
イシューからはじめよを読んだ
急に具合が悪くなるを買ってきたので読んだ
twitterで有名な人の因果推論の勉強会的なイベントに参加した

総括

良い、悪いに関わらず、環境の変化する年であった。

反省

全く笑えないイベントが発生したこともあり、今後のあり方について考える時間を取っておくべきだった。実際には、無意味に残業したり、別の本読んだり、現実逃避してしまっていた。多分現在もそう。

近いうちにやりたいこと

causal inference bookのpart3を読む。*1
または既に読んだ2や1の復習

time-varyingなtreatmentとか自分には関係ないよね、って思っていたけど、一応読んでおくべき？みたいな感覚

失敗の本質を買ってきたので(第二部だけでも)読む

*1:2020年1月に結局読んだ