2022-06-09

Win10からUbuntu経由でPythonを使う

Win10からUbuntu経由でPythonを使おうとしたら思いのほか手間取ったので作業メモを残しておく．

問題：WindowsPowerShellからWSLをインストールできない

原因：AMDのCPUが仮想環境を使える状態になっていない

解決法：BIOS→アドバンス設定→OC→CPUの機能→SVM mode→状態を「無効」から「有効」に切り替え．AMDのコアのデフォルト設定ではSVMが無効になっているようだ．インテルの場合も似たようなCPUの設定を探すと仮想環境対応にできる． CPUの状態が整ったかどうかは，タスクマネージャのCPUのプロパティで確認できる．仮想環境が有効になっていればWindowsPowerShellからWSLをインストールできる．コアの設定後にマイクロソフトストアから好きなディストリビューション（たとえばUbuntu 20.04LTS）を選択してインストールしてもよい． WSLからインストールする場合はディストリビューションまで指定すればうまくいく．

参考：

intint000.blog96.fc2.com

問題：WSLからgitにSSH接続できない

原因：SSH接続用の公開鍵をGithubで設定していない

解決法：鍵を作成して，Githubアカウントで設定する．まずWSLにて

cd ~/.ssh
ssh-keygen

で鍵を作成．clipコマンドだとうまくいかないので

cat ~/.ssh/id_rsa.pub

で中身を表示してコピーする．その後Win上のブラウザからgithubにログインし，マイプロフィールにて公開鍵を設定． WSLに戻ってgitでクローンすると，アクセスが許可される．

なお別のローカルクライアントからgitにSSH接続する場合は，鍵を追加で作成後，githubにログインしてssh keyを追加する．githubにログインせずにシェルから直接公開できるようだが，うまくいかなかった．

参考：

tonariho.com

問題：WSLからjupyter notebookが起動しない

解決法：エクスプローラでwslのファイルにアクセスし，jupyter_notebook_config.pyを直接書き換える

\\wsl$

でWSL上のファイルにWindowsから直にアクセスして

jupyter_notebook_config.py

を探す．メモ帳などで開けたらファイル冒頭に

c.NotebookApp.use_redirect_file = False

を書き込む．

その後WSLからjupyter notebookを入力．するとブラウザから自動でnotebookが開く．編集したファイルはWSL上に保存される．

なおディストリビューションが異なる別のローカルクライアントは，この方法ではうまくいかなかった． jupyter notebook listでローカルアドレスとトークンを表示して，直接ブラウザに打ち込んでアクセスした． configファイルの指定が違うのかもしれない．

参考：

cuxe.livedoor.blog

問題：Win10上のpythonでGinzaを使おうとしたらSudachiPyのインストールで失敗する．

解決策：未解決．Ubuntu（WSL）でCコンパイラを先にインストールするとうまくいった．

感想

自分だけが使うぶんには便利だけど，ちょっと授業では使えないかな．

2022-05-27

OSごとSSDを別のマザーボードに移植する作業の手順

OSがインストールされたM.2SSDドライブを，異なるマザーボードにCPUと一緒に移植するとOSはそのまま起動できるのでしょうか？またチップセットとCPUが異なる場合はどうでしょうか？

実験した結果を記録を残しておきます．

実験1　（CPUは同じだがチップセットのグレードだけ異なるマザボへの移植）

元の構成（メーカーPC）

マザーボード: （Dell製？）チップセットZ370
CPU: CORE-i7 9700
ストレージ: M.2-MVRe SSD500GB　
メモリ: DDR4-16GB
OS: win11

移植先

マザーボード: （ASRock製）チップセットB365

手順

あとでWindowsのPINコードを再設定する必要があるために，WindowsアカウントのログインIDとパスワードを準備しておきます．

またマイクロソフトからの確証コードを受け取るために携帯電話等も必要です．

パスワードを控えたら，もとのマザーボードからCPU，SSD，メモリを外して，B365マザーボードに載せ替えます *1．構成としてはマザボ，CPU，SSD，メモリのみで映像出力はCPU内蔵グラフィックスを使います．キーボードとマウスはUSB2.0接続の無線タイプをつなげました．

実験1の結果

電源スイッチを押すとWindowsは起動しました．しかしログインしようとすると，PINを再設定しろというメッセージが出て，Windowsのアカウント入力を要求されます．ここであらかじめ控えておいた情報が役立ちます（PCが一台しかない場合は，ここで困ると推測）． Windowsアカウントを入力するとSMSで確証コードが携帯電話に送られてきます．コードを入力してPINを再設定すると，Windowsのデスクトップが立ち上がります．結局必要なのはWindowsアカウントの確証のみでした．確証コードの受け取りに携帯電話かタブレットが必要です．確証後の動作は特に問題ありません．

実験2　（CPUもチップセットも異なるマザボへの移植）

元の構成

マザーボード: ASUS P7P55D-E EVO（チップセットP55）
CPU: CORE-i7 870（初代？）
ストレージ: SATA　SSD500GB　
メモリ: DDR3-8GB
OS: win10

移植先

マザーボード: （ASRock製）チップセットB360
CPU: CORE-i7 9700（第9世代）
メモリ: DDR4-16GB

実験2の結果

世代が異なるチップセット（初代から第9世代）のマザボにSSDを載せ替えてもOSは起動した（ちなみにP55は10年以上前に組んだマシンだが現役）
win10のためか，PINの再設定やWindowsアカウントの入力も必要なし．ただし起動後の動作は明らかに遅く，なにか不具合が生じている可能性がある．
数日後にログインすると，Windowsアカウントの確証を求められた．マザボのMACアドレス等で識別しているのかもしれない．

分かったこと

メーカー製のマザボが独自のバックプレートを装着している場合は，トルクスドライバ（T20）で外す．ただしCPUソケットを固定するには専用のステーかナットが必要．このステーはチップセットが異なっても共通なのでジャンクマザボから使い回すことが可能．
高さ164mmのクーラが入るならAS500(PLUS)が最善の選択．他のパーツとの干渉なし，冷える，静か，低価格で最高．OCしないのであればアサシンとかNH-D15は必要なさそう．
CPU用ATX12V8ピンコネクタに，4ピンだけ挿しても起動する．今回使用した電源ユニットで使える補助電源はATX12V4ピンのみ．TDP65WくらいのCPUなら問題なさそう．
PCケースのフロントパネルから電源スイッチだけ外してマザボ（Front_Panel）につなげば，バラックでOSを起動できる（この方法だとCPUクーラの取り付けとOS起動確認が簡単）．
M.2 Wifiコネクタにカードを挿せば，PCIスロットを節約できる．アンテナ線はカードにスナップで取り付け可能（ハンダ付けじゃない）．アンテナ線の先はケース内部に貼りつけでOK．っていうかM.2にWifi接続できるとか，おっさんは知らんし．

極私的結論

ストレージはM.2-MVRe SSD2TBをマザボに2枚挿せば十分．というか起動ドライブもSATA-SSDで十分．
もはやPCケースに5インチベイとHDDベイは不要．つまり165mmの空冷クーラー用ハイトさえ確保できればよい．
マザボの拡張スロットは2個で十分．wifiもM.2接続で十分．つまりマイクロATXで十分．
クーラはAS500かAK400で十分．簡易水冷は不要．　
静音化最後の壁は電源ユニット．ハイブリッドで無音化可能．
BIOSでCPUの電圧下げて静音化しても，パフォーマンスはそれほど低下しない．

おまけ　CPUの設定

RyzenのTDPを下げる方法の一つは，BIOSでEco modeを指定する． MSIの場合，delキーでBIOS設定．Setting→AMD overclocking→Eco mode．Ryzen9-3900xの場合は65Wに下がる．

クロックスピードの調整

BIOSをアドバンスモードに切り替える
CPU Ratio 43→33くらいに倍率を下げる
CPU Core voltage 1.4→1.0（Autoのままでもよい）

なお，DragonCenter（MSIのPCパーツ管理ソフト）でサイレントモードを選ぶとクロックスピードが定格の半分になるようです．

*1:マザボをケースに入れると作業しにくいので，ケースの外に置いたマザボにフロントパネルコネクタから電源スイッチのケーブルだけを挿しておくと便利です．高級なマザボはテスト用にスイッチが基盤に付属してますが，なくても電源スイッチだけ外付けできます．

2021-12-16

ggplot2で複数の折れ線をまとめてplot

すぐ忘れるのでメモ．

library("ggplot2")

y1 <- rnorm(10,5,1)
y2 <- rnorm(10,2,1)
y3 <- rnorm(10,3,1)
years <- 2001:2010

data <-data.frame(years,y1,y2,y3)
ggplot(data, aes(years)) +
  geom_line(aes(y = y1, colour = "y1")) +
  geom_line(aes(y = y2, colour = "y2"))+
  geom_line(aes(y = y3, colour = "y3"))

f:id:hamada7418:20211216173840p:plain

2021-12-14

確率変数の組の独立性について

確率変数の組の独立性

とあるエコノメのテキストに確率変数の組 $$(X_i,Y_i), i=1,2,...,n$$ が独立ならば $$X_j, j=1,2,...,n$$も独立である（大意），と書かれていた．

まあようするに確率変数のペア$(X_i,Y_i)$と$(X_j,Y_j)$が独立ならば，ペアの片方だけを取り出した $X_i$と$X_j$も独立である，ということなのだが，どうも釈然としない．

そもそも $$(X_i,Y_i), i=1,2,...,n$$ が独立であることの定義も書かれていない．

定義して確かめる

テキストに書いてないのだから，自分で定義してもよかろう，ということで1変数の$n$次元独立のアナロジーで，確率変数の順序対の独立を以下のように定義する．

定義（確率変数の組の独立性）

$$(X_i,Y_i), i=1,2,...,n$$ が独立であるとは順序対$(x_1,y_1), (x_2,y_2), ..., (x_n, y_n)$の同時分布について $$ f( (x_1,y_1), (x_2,y_2), ..., (x_n, y_n) )=f(x_1,y_1)f(x_2,y_2)\cdots f(x_n, y_n) $$ が成り立つことである．

なぜこのような定義が必要かと言えば，次のような理由である．上記の定義において$x$と$y$との間には依存関係があってもよい．よって$x$と$y$との間に推測したい関数関係を念頭に置いているが，データとして個人$i$の情報と個人$j$の情報とは無関連であるとき上記の定義が利用できる．

実際，説明変数が確率変数であるような回帰モデルを考える場合には，この仮定が必要である．

さて知りたいことは「 $(X_i,Y_i), i=1,2,...,n$が独立である$\Longrightarrow X_j, j=1,2,...,n$が独立である」が成立するかどうかである．

証明

仮定より $$ f( (x_1,y_1), (x_2,y_2), ..., (x_n, y_n) )=f(x_1,y_1)f(x_2,y_2)\cdots f(x_n, y_n) $$ が成立している．

左辺を$y_1, y_2, ..., y_n$で積分する． $$ \iint \cdots \int f( (x_1,y_1), (x_2,y_2), ..., (x_n, y_n) )dy_1 d y_2 \cdots d y_n=f( x_1,x_2, ..., x_n) $$ 次に右辺も$y_1, y_2, ..., y_n$で積分する． $$ \iint \cdots \int f(x_1,y_1)f(x_2,y_2)\cdots f(x_n, y_n)dy_1 d y_2 \cdots d y_n=f(x_1)f(x_2)\cdots f(x_n) $$

この二つは等しいから $$ f( x_1,x_2, ..., x_n)=f(x_1)f(x_2)\cdots f(x_n) $$ である．

このことはつまり， $$ X_j, j=1,2,...,n $$ が独立であることを意味する．

ふう．

これで安心して $$(X_i,Y_i), i=1,2,...,n$$ が独立であるとき $$ E[u_i|X_1,X_2,...,X_n]=E[u_i|X_i] $$ っていう仮定が使える．

2021-12-02

傾向スコアによるバイアス補正についてのメモ

はじめに

介入の割り当て確率が共変量に依存して異なるとき，統制群と介入群の応答変数の平均の差によって介入効果を推定するとバイアスが生じます．もし，共変量の値毎に割り当て確率を決定する関数が分かっていると，推定量のバイアスを取り除くことができます．

ただし，確率を決定する関数が分かっている点がポイントで，テキトーに関連がありそうな共変量でロジスティック回帰をすればいいわけではありません．

以下にテキトーなロジスティック回帰だと失敗する例を示します．

バイアス割り当て実験

分析に使うパッケージの読み込み

library("MatchIt")
library("Matching")
library("cobalt")
library("WeightIt")
library("broom")

共変量xによって介入（トリートメント）の割り当て確率が異なる状況を考えます．

xは-5から5までの範囲に値をとる一様分布とします．潜在的結果y0は，この共変量と相関しています．

真の介入効果はtau=3とします．よって潜在的結果y1はy1=y0+3で定義できます．

n <- 1000;# sample size   
tau <- 3;# 真の介入効果    
x <- runif(n,-5,5)# 共変量
y0 <-5+2*x+ rnorm(n,mean=0,sd=1)#潜在的結果yo, xの影響を受ける
y1 <- y0+tau #潜在的結果y1

介入の割り当て確率をロジスティック関数で定義します．

p <- 1/(1+exp(-x))# make probabilities by logistic function
d <- rbinom(n,size=1,p)# treatment with bias

pは割り当て確率です．共変量xの値が大きいほど割り当て確率が高くなっています．

dは介入を表すダミー変数です． dを確率pのベルヌーイ分布で定義するため，xの値が大きい個体ほど介入を受けやすくなります．

以上のコードで割り当てにバイアスがある状況を表現できました．

観測データyを以下のように定義します*1．

y <- c()
  for(i in 1:n){# 観測値yの生成
    if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}}
data <- data.frame(y,d,p,x)

介入があった個体については潜在的結果y1を，介入がなかった個体についてはy0を応答変数として観測します．各個体は，y0もしくはy1のどちらかしか観測できないことに注意します（因果推論の根本問題）

共変量xと割り当て確率pの関係を図で確認しておきましょう．

以下，割り当て確率のことを「傾向スコア」と呼びます*2．

f:id:hamada7418:20211202225451p:plain — 共変量`x`と傾向スコア`p`の関係

共変量xと介入dの関係も図で確認しておきましょう．

f:id:hamada7418:20211202225737p:plain — 共変量xと介入dの関係

共変量xの値が小さいほど介入の割り当てがなく，値が高いほど介入が割り当てられています．

このようなランダムでない割り当ては，介入効果の推定にバイアスをもたらします．

計算して確かめてみましょう．

> summary(lm(data=data,y~d))

Call:
lm(formula = y ~ d, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.3599  -2.7284   0.0375   2.8791  11.4583 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8242     0.1699   4.852 1.42e-06 ***
d            11.8107     0.2384  49.549  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.768 on 998 degrees of freedom
Multiple R-squared:  0.711,  Adjusted R-squared:  0.7107 
F-statistic:  2455 on 1 and 998 DF,  p-value: < 2.2e-16

介入効果の推定値は11.8107でした．真の効果は3ですから，過大推定になっています．

回帰ではなく，条件付き期待値の差をとっても同じ結果を得ます．

mean(y[d==1])-mean(y[d==0])

傾向スコアマッチング

傾向スコアはすでに計算済みですから，それをつかってマッチングデータを作り，ATEを推定します．

> matched.data <- Match(Y = data$y, Tr=data$d, X =p,estimand="ATE")
> summary(matched.data)#Matched number of observations...............  500 

Estimate...  3.4254 
AI SE......  0.23556 
T-stat.....  14.541 
p.val......  < 2.22e-16 

Original number of observations..............  1000 
Original number of treated obs...............  508 
Matched number of observations...............  1000 
Matched number of observations  (unweighted).  1152

直感的に言えば，関数Matchは「傾向スコアが同じだけど介入の有無が異なる個体をペアにしたデータ」を作ります．

結果，バイアスが減り，真の効果に近いATE推定値を得ることができました．念のためcobaltパッケージで共変量の調整結果を可視化します

data2 <- matchit(d~x,data=data,method = "nearest",replace = T)
love.plot(data2, thresholds = 0.1)#SMD

f:id:hamada7418:20211202231936p:plain — 調整前後での共変量平均値の比較

傾向スコアによるマッチングで共変量のバランスがとれたことが分かります．

真の傾向スコアが手に入らない場合

上記の実験では，真の傾向スコアが分かっている条件下でマッチングデータを作りました．しかし現実の観察データを使って分析する場合，傾向スコアが事前に分かっていることなどありません．そもそも傾向スコアを決定する共変量の全てがデータとして入手できる状況など，極めて希でしょう．

そこで次に，必要な共変量が手に入らない場合にGLMで傾向スコアを推定してマッチングに利用した場合，どのくらい真の効果量の推定にバイアスが生じるのかを確かめます．まず分析用のデータの生成です．

n <- 3000;# sample size   
tau <- 3;# 真の介入効果    
x1 <- runif(n,-5,5)# 共変量
x2 <- runif(n,-5,5)# 共変量
y0 <-x1+x2#潜在的結果y0, x1,x2の影響を受ける
y1 <- y0+tau #潜在的結果y1
p <- 1/(1+exp(-x1-x2))# make probabilities by logistic function
# pは真の傾向スコア．共変量x1とx2の値に依存する．
d <- rbinom(n,size=1,p)# treatment with bias
y <- c()
for(i in 1:n){# 観測値yの生成
  if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}}

真の傾向スコアは共変量x1とx2によって定まると仮定しています．次に，必要な共変量が入手できない場合とできた場合を比較します．

out.glm <- glm(d~x1,family = "binomial")#傾向スコアを過小定式化モデルで推定
out.glm2 <- glm(d~x1+x2,family = "binomial")#フルモデルで推定

p.predict <- out.glm$fitted.values
p.predict2 <- out.glm2$fitted.values

p.predictはx2が足りない場合の傾向スコアの推定値で，p.predict2はx1とx2が正しくGLMに含まれる場合の推定値です．結果を比較してみます．

まず共変量が足りない場合

matched.data <- Match(Y = y, Tr=d, X =p.predict,estimand="ATE")
summary(matched.data)

Estimate...  7.2681 
AI SE......  0.10906 
T-stat.....  66.645 
p.val......  < 2.22e-16 

Original number of observations..............  3000 
Original number of treated obs...............  1509 
Matched number of observations...............  3000 
Matched number of observations  (unweighted).  7892

真の効果3に対して推定値は7.2681ですから過大推定しています（GLMのMcFadden R² は 0.2209495でした）．

このとき共変量のバランスをチェックすると，共変量x1に関しては調整後にバランスはとれています．

dat <- data.frame(d,x1,x2,y)
match.out <- matchit(data=dat,d~x1,method = "nearest",replace = T)
love.plot(match.out, thresholds = 0.1)

f:id:hamada7418:20211215125227p:plain — 調整後の共変量

つまり「投入した共変量のバランスがとれたとしても，効果量の推定には失敗する可能性がある」ので注意が必要です．

次に必要な共変量が入手できる場合のマッチング推定の結果です．

 matched.data2 <- Match(Y = y, Tr=d, X =p.predict2,estimand="ATE")
 summary(matched.data2)

Estimate...  3.4448 
AI SE......  0.1594 
T-stat.....  21.611 
p.val......  < 2.22e-16 

Original number of observations..............  3000 
Original number of treated obs...............  1509 
Matched number of observations...............  3000 
Matched number of observations  (unweighted).  4921

バイアスが補正されました（GLMのMcFadden R² は 0.5866154でした）．

実際の観測データを用いた分析では，傾向スコアの正しい定義に必要な共変量は不明です．データに含まれる共変量のバランスがとれたとしても，処置への割り当て確率を計算するモデルに欠落変数があると，処置効果の推定に失敗します．つまり処置に影響しそうな共変量を分析者が（勝手に）考えて傾向スコアを予測しても，うまくいく保証はありません*3．

補足

必要な共変量が揃っている場合の，誤差項の分散の影響も確かめてみました．複数回マッチングさせるために関数化します．

f1 <- function(s){
  n <- 1000;# sample size    
  tau <- 3;# 真の介入効果  
  x <- runif(n,-5,5)# 共変量
  y0 <-x+ rnorm(n,mean=0,sd=3)#潜在的結果yo, xの影響を受ける
  y1 <- y0+tau #潜在的結果y1
  p <- 1/(1+exp(-x-rnorm(n,mean=0,sd=s)))# make probabilities by logistic function
  # pは真の傾向スコア．共変量xの値が大きいほど割り当て確率が高い
  d <- rbinom(n,size=1,p)# treatment with bias
  y <- c()
  for(i in 1:n){# 観測値yの生成
    if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}}
  out.glm <- glm(d~x,family = binomial)
  Lm <- out.glm$deviance/(-2)
  L0 <- out.glm$null.deviance/(-2)
  p.predict <- pout.glm$fitted.values
  #予測した傾向スコアを使いMatchでマッチさせる
  matched.data <- Match(Y = y, Tr=d, X =p.predict,estimand="ATE")
 c(matched.data$est[1], 1-Lm[1]/L0[1])# McFadden R^2
}

7行目p <- 1/(1+exp(-x-rnorm(n,mean=0,sd=s)))の部分で，共変量xから傾向スコアを計算する過程で誤差を導入します．誤差項の標準偏差sが大きくなると，共変量xによって予測できる傾向スコアの精度が低下します．

f:id:hamada7418:20211212001754p:plain — 共変量`x`と傾向スコア`p`の関係．誤差項の標準偏差`s=3`

トリートメントの割り当て確率は，この傾向スコアによって決まるので，割り当てdを共変量xに回帰すると傾向スコアの予測値を得ます．

最後にMatch関数に予測した傾向スコアを与えてマッチングデータを作成します．計算の手続きをまとめ関数f1の引数は，誤差の標準偏差sのみです．アウトプットは「マッチングデータから推定したATE」と「傾向スコアを予測したロジスティック回帰モデルのMcFadden 疑似決定係数R²」です．

計算例を確認してみましょう．

 Mac <- c()
 Est <- c()
 for(i in 1:50){
   Est <- append(Est,f1(s=3)[1])
   Mac <- append(Mac,f1(s=3)[2])
 }
 
> mean(Mac)
[1] 0.2395028
> mean(Est)
[1] 2.998437

McFadden疑似決定係数が0.24程度でも，真の効果をうまく推定できています．

 Mac <- c()
 Est <- c()
 for(i in 1:50){
   Est <- append(Est,f1(s=10)[1])
   Mac <- append(Mac,f1(s=10)[2])
 }
 
> mean(Mac)
[1] 0.03893244
> mean(Est)
[1] 3.002857

McFadden疑似決定係数が0.039程度でも，真の効果をうまく推定できています．この計算例から，疑似決定係数の大きさはバイアス補正にとって本質的ではないことが分かります．たとえ疑似決定係数が高くても，必要な共変量が足りなければ推定には失敗します．

逆確率による重み付け

傾向スコアを使ったIPWによる調整も試しておきます．ウェイトの計算にはWeightItパッケージを使いました． weightit関数がdataを引数にとるのでdata.frameを定義しておきます．

データセットの定義

n <- 3000;# sample size   
tau <- 3;# 真の介入効果    
x1 <- runif(n,-5,5)# 共変量
x2 <- runif(n,-5,5)# 共変量
y0 <-x1+x2#潜在的結果y0, x1,x2の影響を受ける
y1 <- y0+tau #潜在的結果y1
p <- 1/(1+exp(-x1-x2))# make probabilities by logistic function
d <- rbinom(n,size=1,p)# treatment with bias
y <- c()
for(i in 1:n){# 観測値yの生成
  if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}}
dat <- data.frame(d,x1,x2,y)

過小定式化モデル（共変量x2が足りないモデル）と適切なモデルによって傾向スコアを計算し，それぞれ重みに変換して回帰に使用します．

## 重み付きデータでの効果の推定 傾向スコアを過小定式化モデルで推定
weighting <- weightit(data=dat, d ~ x1,method = "ps",estimand = "ATE")
IPW_result <- lm(y ~ d,weights = weighting$weights) %>% tidy()
IPW_result


## 重み付きデータでの効果の推定 フルモデルで推定
weighting2 <- weightit(data=dat, d ~ x1+x2,method = "ps",estimand = "ATE")
IPW_result2 <- lm(y ~ d,weights = weighting2$weights) %>% tidy()
IPW_result2

結果の比較です．

> IPW_result
# A tibble: 2 x 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    -2.16    0.0679     -31.9 4.56e-192
2 d               7.45    0.0959      77.7 0.       


> IPW_result2
# A tibble: 2 x 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)   -0.269    0.0923     -2.92 3.58e-  3
2 d              4.02     0.134      30.1  6.63e-174

共変量が足りないモデルではうまくバイアスが補正できませんでした．共変量がそろったモデルではバイアスが（相対的にはうまく）補正できました．ただし傾向スコアが0に近い値がデータセットに多く含まれるせいで，推定結果にはまだバイアスが残っています．

　参考文献

安井翔太，2020『効果検証入門』技術評論社
ローゼンバウム，2017=2021『統計的因果推論入門』共立出版
星野崇宏，2009『調査観察データの統計科学　因果推論・選択バイアス・データ融合』岩波書店

*1:データフレームはなくても構いませんが，ここでは他のコードで利用するために定義しています

*2:傾向スコアとは「共変量Xの値がxであるすべての個体iに関する処置群への割り当て確率の平均値」です

*3:もちろん，偶然成功する場合もあります

2021-11-05

数理モデル解説本の内容紹介

はじめに

2018年12月に『その問題，数理モデルが解決します』をベレ出版さんから刊行しました．

当blogではこの本の内容紹介をしていなかったので，以下にざっと解説します*1．

その問題、数理モデルが解決します

作者:浜田宏
ベレ出版

Amazon

本書は数学，統計学，経済学（行動経済学，ゲーム理論），社会学（数理社会学）から，筆者が特におもしろいと思う《人や社会に関する数理モデルやアルゴリズム》を紹介しています．2人のキャラクターが対話しながら「ある問題」を解決するための「数理モデル」や「アルゴリズム」を考えます．基本的には1章で1話が完結します．

続編『その問題，やっぱり数理モデルが解決します』と世界観は共通ですが，内容は独立しているので，どちらから読み始めていただいても問題ありません．概ね高校・大学初年度レベルの数学を使っています．

特徴

単にこんなモデルがあるよ，と紹介するだけでなく，その数学的な原理をフォローできるように証明や計算を省略せずに書きました．この本で紹介するモデルやアルゴリズムは，「一度はどこかで聞いたことがあるもの」であり，目新しさはないかもしれません．本書は「こういう命題は知ってるけど，どうして成立するんだろう」とか「こういうアルゴリズムは知ってるけど，どうして有効なんだろう」という疑問に答えます．

たとえば，5章の「秘書問題」の解決アルゴリズムが「全体の36.8%を観察してから，次にベストを超えた対象を選ぶ」であることは，よく知られています．しかしどうして36.8%（より正確に書くと1/e）を見送る必要があるのか，その《しくみ》を説明することは簡単ではありません．本書は，強力なアルゴリムが成立する《しくみ》をなるべく丁寧に解説します．

対象読者

本書の読者として以下のような方を想定しています．

数学って何の役に立つの？と考えている中高生
「数理モデル」*2に興味のある大学生
統計学を使うけど，仕組みがよく分からないという研究者
データ分析を勉強しようかなと考えている社会人

各章のタイトル・難易度・トピック

章タイトル	難易度	トピック
序	☆	モデルとはなにか
第1 章　隠された事実を知る方法	☆☆	回答のランダム化、集合、確率変数、期待値、分散
第2 章　卒業までに彼氏ができる確率	☆☆	ベルヌーイ分布、確率変数の合成、2 項分布
第3 章　内定をもらう方法	☆☆☆	2 項分布の期待値、確率変数の和の期待値、インプリケーション、ベータ分布、ベータ2 項分布
第4 章　先延ばしをしない方法	☆☆	先延ばしのメカニズム、時間割引、準双曲型割引、先延ばしの防止
第5 章　理想の部屋を探す方法	☆☆☆	秘書問題、グーゴル・ゲーム、最適停止問題、アルゴリズム、全体の36.8% を見送る理由、生活満足度と通勤時間
第6 章　アルバイトの配属方法	☆	選好、安定マッチング、DA アルゴリズム、パレート効率性
第7 章　売り上げをのばす方法	☆☆☆	ランダム化比較試験、条件付き期待値、潜在的結果、不偏推定量、統計的検定
第8 章　その差は偶然でないと言えるのか？	☆☆☆☆☆	検定の仕組み、棄却域、対立仮説、正規分布の性質、サンプルサイズの設計
第9 章　ネットレビューは信頼できるのか？	☆☆☆☆	陪審定理、チェビシェフの不等式、大数の弱法則、陪審定理の一般化
第10 章　なぜ0 円が好きなのか？	☆☆☆	ゼロ価格アノマリー、効用関数と導関数、プロスペクト理論、価値関数、ゼロ価格効果の一般化
第11 章　取引相手の真意を知る方法	☆☆☆	ゲーム理論と支配戦略、第2 価格封印入札、メカニズムデザイン、ナッシュ均衡
第12 章　お金持ちになる方法	☆☆☆	ギャンブルでお金持ちになる方法、所得分布のカタチ、累積効果、対数正規分布の生成

以下，各章の概要です．

各章の内容

第1 章　隠された事実を知る方法

規則やモラルに反する行動の調査では，回答者が正直に答えないことが知られています．そこで「禁煙の場所と知りながら煙草を吸ったことがある」などの逸脱行動を把握するために《回答のランダム化》という手法を使って解決を試みます（この手法のもとで母比率を推定するためには，実は《大数の法則》という命題が必要です．その解説は9章で詳しくかきました）・

第2 章　卒業までに彼氏ができる確率

100人と出会ったとき，1人以上から好かれる確率はどのくらいだろうか？この確率は2項分布をつかったモデルで計算できます． 2項分布は統計学でよく使う基礎的な部分の1つであり，より単純なベルヌーイ分布の合成から《つくる》ことができます．その過程を詳しく解説しました．

第3 章　内定をもらう方法

第2章の内容を展開します．モデルからインプリケーションを引き出したり、ベータ分布を使ったモデルの拡張について考えます．

第4 章　先延ばしをしない方法

先延ばし行動に関する心理学や行動経済学の知見を紹介しつつ，時間割引モデルについて解説します．先延ばし防止に役立つ知識もまとめています．

第5 章　理想の部屋を探す方法

いわゆる「秘書問題」モデルを解説しています．ファーガソンのレビュー論文内の証明を，できるだけ丁寧に再構成しました．他の章と比べて少しだけ難易度が高くなっていますが，1/eの謎が解けるはずです．コンピュータで計算するための簡単なコードも記載しています．こちらを参照

第6 章　アルバイトの配属方法

ゲールとシャプレーが提唱したマッチングアルゴリズムを紹介します．アルゴリズム自体は簡単なのですが，なぜこれが安定マッチングを生むのか，その理屈（証明）を詳しく書きました．

第7 章　売り上げをのばす方法

ランダム化比較試験をルービンの因果推定モデルをベースにして紹介しました．平均処置効果の推定では条件付き期待値の計算がとても重要なので，その基礎部分を詳しく書きました．

第8 章　その差は偶然でないと言えるのか？

検定の仕組みを解説しました．サンプルサイズを計算する式を見たことのある人は多いと思いますが，あの式がどこから出てくるのか疑問に思うことでしょう．私自身もどういう理屈なのか気になったので，いろいろ調べて理解した内容をここで解説しています．

第9 章　ネットレビューは信頼できるのか？

陪審定理について解説しています．人文社会系では有名な定理ですが、証明は書かれていないことが多いように思いました．第1章で使った大数の法則を利用して、陪審定理の証明を説明します．また少しだけ定理の一般化についても考えています．

第10 章　なぜ0 円が好きなのか？

《100円の商品の10円の値引き》と《10円の商品の無料化》はどちらも値引き額は10円ですが，後者の方が人を引きつけることが知られています．どうして同じ値引き額なのに，われわれは無料を選ぶのでしょうか？行動経済学の実験を紹介しつつ，効用関数やプロスペクト理論の価値関数について学びます．

第11 章　取引相手の真意を知る方法

第2 価格封印入札というオークションのルールがあります．これは「入札額が一番高い人が勝者だが，支払い額は2番目に高い入札額となる」というルールです．このルールを使うと，理論上は，入札者は評価額を正直に表明するインセンティブを持つことが知られています．なぜそうなるのか？という理屈をゲーム理論のモデルを使って説明します．

第12 章　お金持ちになる方法

所得の分布はパレート分布や対数正規分布という確率分布で近似できることが知られています．しかし世の中の人々は，所得分布の形など意識することなく，生活していることでしょう．なぜ所得分布が特定の確率分布で近似できるのか，そのしくみを簡単な数理モデルで説明します．このモデルの原型は私が大学院生のときに作りました．

メッセージ

私は大学院に入るまでは，「自分は文系だから数学は人生に無用だ」と考えていました．

しかし師の影響で数理モデルの勉強をはじめたことがきっかけとなり，モデルの楽しさに目覚めました．これ以上の娯楽は，なかなかない，と今では思います．

本書をきっかけに，数理モデル（数学）っておもしろいかも，と思っていただけたら望外の喜びです．

*1:続編（『その問題，やっぱり数理モデルが解決します』2020年9月刊行）との違いはなんですか？という質問を読者から時々いただきます．この紹介を参考にしてください．

*2:数理モデルは，数学を使って現象を表現・説明・予測する方法の総称であり，観察した現象から本質を見抜く洞察やデータを分析するための指針（理論）を与えてくれます．現象が成立するメカニズムを数理モデルで特定できれば，現象の理解や未来の予測に役立ちます．

2021-10-21

自由エネルギーの計算例2

実験2

前回の続きです．

データを生成する分布として次のような複合分布を考えます．

$$ X \sim \mathrm{Bernoulli}(q) $$ $$ q \sim \mathrm{Beta}(a, b) $$ これはベルヌーイ分布のパラメータ$q$の事前分布がベータ分布$\mathrm{Beta}(a, b)$にしたがう複合分布です．

パラメータ$q$を消去すればこの複合分布はベータ2項分布 $$ X \sim \mathrm{Betabinomial}(1,a,b) $$ と一致します．

このベータ2項分布からデータを生成し，「正しい」統計モデルで分析します．

つまり $$ X \sim \mathrm{Bernoulli}(q) $$

$$ q \sim \mathrm{Beta}(a, b) $$ というベイズ統計モデルを仮定して自由エネルギーを計算したらどうなるかを実験して確かめます．　

まずデータを生成するためにRのextraDistrパッケージを使い，ベータ2項分布にしたがう確率変数の実現値を100個つくります．

> library(extraDistr)
> rbbinom(100, 1, alpha = 2, beta = 8)
  [1] 1 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0
 [43] 0 0 1 0 1 0 0 0 0 0 1 1 0 0 0 0 0 1 1 0 0 0 1 0 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0
 [85] 1 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0

このデータを使って自由エネルギーを計算します．

> Fn(x,a=2,b=8)
[1] 49.84582
> Fn(x,a=20,b=80)
[1] 48.97739
> Fn(x,a=200,b=800)
[1] 48.69709

念のため，Mathematicaでも同種のデータを生成します．

In[1]:= (* ベータ2項分布からの実現値生成 *)
x = RandomVariate[ BetaBinomialDistribution[2, 8, 1], 100]

Out[1]= {1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, \
1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, \
0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, \
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1}

In[2]:= Mean[x]

Out[2]= 9/50

自由エネルギーを定義して計算してみましょう．

In[4]:= FEnergy[x_, a_, b_] :=
  N[-1*Log[Beta[Total[x] + a, Length[x] - Total[x] + b]/Beta[a, b]]];

In[5]:= 
FEnergy[x, 2, 8]
FEnergy[x, 20, 80]
FEnergy[x, 200, 800]
FEnergy[x, 2000, 8000]
FEnergy[x, 20000, 80000]

Out[5]= 48.3529
Out[6]= 47.5336
Out[7]= 47.3005
Out[8]= 47.271
Out[9]= 47.268

さすがはMathematica．パラメータが20000，80000でも余裕で計算してくれます（Rだと200，800が限界でそれ以上は計算してくれない．そういうとこやぞ）．

どちらも同じような結果になりました．さてデータを生成する分布として$\mathrm{BetaBinomial}(n=1,a=2,b=8)$を仮定していますが，データから計算した平均の最尤推定値は9/50です．

次に事前分布の平均値が9/50になるようなモデルの自由エネルギーを計算してみます．

N[FEnergy[x, 180000, 820000], 10 ]
N[FEnergy[x, 1800000, 8200000], 10 ]
N[FEnergy[x, 18000000, 82000000], 10 ]

47.13939868
47.13935368
47.13934918

真のモデルよりも小さい？　自由エネルギーは確率変数なので，データの出方によってばらつきがあります．以上の計算をふまえ，データを複数回生成することで自由エネルギーの平均を比較します．

out <- c()#真の事前分布beta(2,8)の自由エネルギー
for(i in 1:3000){
  x <- rbbinom(1000, 1, alpha = 2, beta = 8)
  out <- append(out,Fn(x,a=2,b=8))
}

out2 <- c()# 事前分布beta(20,80)の自由エネルギー
for(i in 1:3000){
  x <- rbbinom(1000, 1, alpha = 2, beta = 8)
  out2 <- append(out2,Fn(x,a=20,b=80))
}

out3 <- c()# MLの自由エネルギー
for(i in 1:3000){
  x <- rbbinom(1000, 1, alpha = 2, beta = 8)
  out3 <- append(out3,-1*log(mean(x)^(sum(x))*(1-mean(x))^(length(x)-sum(x))))
}

結果を比較すると......

最尤法モデル＜ベイズモデル beta(20,80) ＜ベイズモデル beta(2,8)　

となりました．真の分布と同じ統計モデル（確率モデル＋事前分布）の平均自由エネルギーより，最尤推定モデルに近いベイズモデルや最尤推定モデルの平均自由エネルギーのほうが小さいことが分かりました*1．

統計モデルが真の分布に一致している場合よりも，事前分布が最尤推定値に質点を持つディラックデルタ分布に近い方が，この条件下では自由エネルギーの平均値が小さくなっています．このことから自由エネルギーによるモデル選択は，データを生成した分布よりも，観察したデータそのものに尤度の意味で近い統計モデルを選んでいるらしいと考えられます．なぜこうなってしまうのでしょうか？引き続き考えてみます．

実験3

さきほどの実験では，$\mathrm{BetaBinomial}(n=1,a=2,b=8)$からデータを生成しましたが，そもそもこのデータでは複合分布を周辺化した結果がベルヌーイ分布と一致するので例としては特殊すぎた可能性があります．

違う条件で計算してみるとどうなるのでしょうか？自由エネルギーによるモデル選択は最大対数尤度による選択と一致するのでしょうか？

実験3の設定は次の通りです．はじめにデータを $$ X \sim \mathrm{Poisson}(\lambda) $$

$$ \lambda \sim \mathrm{Gamma}(a, b) $$ という結合分布から生成します．よく知られているように，この仮定の下で$\lambda$を周辺化すると，$X$は負の2項分布に従うことが知られています．したがってデータを $$ X \sim \mathrm{NegativeBinomial}(a, 1/(b+1)) $$ から生成したと考えても同じことです．以下では$a=1, b=3$を仮定してデータを生成します（$n=100$のデータを作りました）．

In[3]:= x = RandomVariate[NegativeBinomialDistribution[1, 0.25], 100]
Out[3]= {5, 0, 2, 1, 3, 5, 2, 9, 2, 0, 8, 2, 13, 2, 2, 0, 6, 0, 3, 6, \
4, 0, 1, 4, 1, 19, 2, 8, 0, 1, 0, 1, 5, 4, 3, 0, 5, 0, 9, 0, 3, 0, \
17, 1, 10, 5, 2, 8, 1, 0, 2, 2, 4, 0, 0, 9, 2, 0, 0, 3, 3, 2, 5, 0, \
0, 2, 0, 2, 1, 1, 4, 0, 7, 0, 1, 21, 5, 6, 0, 0, 0, 7, 4, 0, 0, 0, 0, \
0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0}

このデータをベイズモデル $$ X \sim \mathrm{Poisson}(\lambda) $$

$$ \lambda \sim \mathrm{Gamma}(a, b) $$ で学習して自由エネルギーを計算します．するとこのベイズモデルの自由エネルギーは $$ F_n=-\log\Gamma\left(\sum x_i+a\right) +\log \Gamma(a)+\sum \log x_i ! -a\log b +\left(\sum x_i +a\right)\log (n+b) $$ です（浜田・石田・清水2019: 111）．

自由エネルギーの引数はデータ$x_1, x_2, ... x_n$とパラメータ$a, b$です．このうちデータと片方のパラメータの真値$a=1$が与えられたとして，$b$の変化による自由エネルギーの変化を確認してみましょう．

f:id:hamada7418:20211019163804p:plain — 自由エネルギーの変化

データを生成したパラメータは$a=1, b=3$ですから，$b=3$のとき最も小さくなることを予想していましたが，全く違う箇所で最小化しているようです．

次にデータを生成した負の2項分布ではなく，ポアソン分布を確率モデルとして使った場合の自由エネルギーを計算します．事前分布としてデータの平均の最尤推定値に質点をもつようなディラックデルタ分布を仮定します．つまりポアソン分布で最尤推定した場合の対数尤度の符号反転を自由エネルギーとして定義します．

ポアソン分布のパラメータ$\lambda$の変化にともなう自由エネルギーの変化は次の通りです．

f:id:hamada7418:20211019163636p:plain — 自由エネルギー（モデルにポアソン分布を仮定したとき）

このデータの平均の最尤推定値は2.84なので，自由エネルギーの定義上，$\lambda=2.84$のとき最小化します．

サンプルサイズ1000のデータセットを3000回発生させ，自由エネルギーを計算した結果を比較します．もはやRでは実行できないのでMathematicaに計算してもらいます．

Module[{xvec, out, out2},
 out = {samplesize = 1000, itr = 3000};
 Do[
  xvec = RandomVariate[NegativeBinomialDistribution[1, 0.25], 
    samplesize];
  out = Append[out, Fn[xvec, 1, 3]], {itr}
  ];
 out2 = {};
 Do[
  xvec = RandomVariate[NegativeBinomialDistribution[1, 0.25], 
    samplesize];
  out2 = Append[out2, Fnpois[xvec, Mean[xvec]]], {itr}
  ];
 Print[{Mean[out], Mean[out2]}];
 Histogram[{out, out2}]
 ]

f:id:hamada7418:20211021150820p:plain — 自由エネルギーの分布の比較．オレンジは負の2項分布，青はポアソン分布

自由エネルギーの平均値は，

ポアソンモデル　2899.84　<　データを生成した分布（負の2項分布）　2910.24

でした．この条件でも平均自由エネルギーはデータを生成したモデルを選択しません．定義上，自由エネルギーの平均値は真のモデルと周辺尤度の交差エントロピーなので，周辺尤度を真の分布と一致するように作れば，自由エネルギーの平均も最小化するはずです．

式で書くと $$ \begin{align} E_{q(X^{n})}[F_n]&=-\int q(x^{n})\log p(x^{n}) dx^{n} \\ &=H[q(X^{n})]+D[q(X^{n})||p(X^{n})] \end{align} $$ なので，自由エネルギーの平均は真の分布のエントロピーと真の分布と周辺尤度のKL距離に分解できます．このうちKL距離は真の分布$q(X^{n})$と周辺尤度$p(X^{n})$ が一致しているのでゼロです．

しかし上記の条件では，データを発生させる真の分布とは異なる確率モデルのもとで，自由エネルギーがより小さな値をとっています．

なぜこうなるのか正直よく分かりません．原因が判明したら追記します．

*1:細かいことを言うと，ここでいう最尤推定モデルの自由エネルギーとは，事前分布が最尤推定値で質点を持つディラックデルタ分布であるようなベイズモデルの自由エネルギーと定義しておきます．これは実質的には最大対数尤度の－1倍と同じです

問題：WindowsPowerShellからWSLをインストールできない

問題：WSLからgitにSSH接続できない

問題：WSLからjupyter notebookが起動しない

問題：Win10上のpythonでGinzaを使おうとしたらSudachiPyのインストールで失敗する．

感想

実験1 （CPUは同じだがチップセットのグレードだけ異なるマザボへの移植）

元の構成（メーカーPC）

移植先

手順

実験1の結果

実験2 （CPUもチップセットも異なるマザボへの移植）

元の構成

移植先

実験2の結果

分かったこと

極私的結論

おまけ CPUの設定

確率変数の組の独立性

定義して確かめる

定義（確率変数の組の独立性）

証明

はじめに

バイアス割り当て実験

傾向スコアマッチング

真の傾向スコアが手に入らない場合

補足

逆確率による重み付け

参考文献

はじめに

特徴

対象読者

各章のタイトル・難易度・トピック

各章の内容

第1 章 隠された事実を知る方法

第2 章 卒業までに彼氏ができる確率

第3 章 内定をもらう方法

第4 章 先延ばしをしない方法

第5 章 理想の部屋を探す方法

第6 章 アルバイトの配属方法

第7 章 売り上げをのばす方法

第8 章 その差は偶然でないと言えるのか？

第9 章 ネットレビューは信頼できるのか？

第10 章 なぜ0 円が好きなのか？

第11 章 取引相手の真意を知る方法

第12 章 お金持ちになる方法

メッセージ

実験2

実験3

実験1　（CPUは同じだがチップセットのグレードだけ異なるマザボへの移植）

実験2　（CPUもチップセットも異なるマザボへの移植）

おまけ　CPUの設定

　参考文献

第1 章　隠された事実を知る方法

第2 章　卒業までに彼氏ができる確率

第3 章　内定をもらう方法

第4 章　先延ばしをしない方法

第5 章　理想の部屋を探す方法

第6 章　アルバイトの配属方法

第7 章　売り上げをのばす方法

第8 章　その差は偶然でないと言えるのか？

第9 章　ネットレビューは信頼できるのか？

第10 章　なぜ0 円が好きなのか？

第11 章　取引相手の真意を知る方法

第12 章　お金持ちになる方法