論理の流刑地

地獄の底を、爆笑しながら闊歩する

RのプログラムをMacからWindowsに移行する時に使うコード

Mac R Windows

備忘 of 備忘 of 備忘

Motivation

ちょっと色々あってMacからWindowsに乗り換えることになり（8年振りくらい）、
自分で作った関数や色々な分析の記録も含めRのコードをそのまま新しい環境にうつした。
基本的にはOSに依存しないで動くものであるので、あまり考えるところはない。

しかし、置かれるディレクトリの場所がどうしても変わってしまい、source("")やload("")などのdirectoryの指定先が変わるので、そこに関しては一括で書き換えなくてはいけない。
具体的には、Macだと

source("/Users/Name1/xxx.R")

となっていたのを、Windowsでは

source("C:/Users/Name2/yyy/xxx.R")

みたいに該当のフォルダにかえなくてはならない。

この処理が単純そうに見えて意外と思いつかなかったのでメモしておく。

Implementation


How to read txt file and replace word in R? - Stack Overflow
↑このURLを参考にした。

ファイルのパスの一覧の取得

まず、Rのファイル一式のパスを取得する。
自分はRのコードファイルの拡張子を「.R」にしているので、それを正規表現で拾ってきている。

dir_tgt <- list.files( "C:/Users/Name1/Documents/R/", pattern = ".*\\.R" , recursive =T , include.dirs=F)
dir_list <- paste0( "C:/Users/Name1/Documents/R/" , dir_tgt)

recursive=Tで配下のディレクトリまで探知範囲にすることを、include.dirs=Fでフォルダのパスは取得しないことを指示している。
二行目では後々のためにフルパスにしている。

つぎに、中身を一括で書き換えるための関数をつくる

中身を書き換える関数を実装する

ファイルのパスを指定すれば、load()やsource()で使われている中身のディレクトリを書き換えてくれる関数をつくる

reToWin <- function(tdir){
old_txt <- readLines( tdir , encoding = "UTF-8")
new_txt <- gsub("/Users/Name2/Documents/" , "C:/Users/Name1/Documents/", x = old_txt)
writeLines(new_txt , con = tdir , useBytes=T)
} #function

readLines()のほうのencoding = "UTF-8"と、writeLines()のuseBytes=Tをいれないと自分の環境では文字化けした。
多くのMac→Windowsの移行でも一緒だと思う（たぶん）

作った関数で一括して変更する

for( i in seq_along( dir_list)){
reToWin(dir_list[i])
} #for


終わり！わかれば簡単。

あとはMac環境でsave()したものをWindowsでload()すると文字化けするみたいな問題もあるが、
それはその度にiconv()とかで対応していくしかないかなぁ....


Def Tech - Catch The Wave

積み上げ折れ線グラフをggplot2で描く

R

Motivation

ggplot2()は割と綺麗なグラフが描けるので重宝しているが、
いかんせんそこまで高頻度で使うわけでないので、忘れてしまう。

備忘として、積み上げ折れ線グラフの作りかたを書き記しておく。

参考URL

URL1
ill-identified.hatenablog.com
URL2
mukkujohn.hatenablog.com

上のブログだとgeom_line()とgeom_ribbon()を用いた方法が、下のブログではgeom_area()を用いた方法が紹介されている

やってみる

データ準備

データに関しては、参考URL1で用いられているSeatbeltsデータを使う

data(Seatbelts) 
df <- data.frame(Seatbelts, year=rep(1969:1984,each=12) ) 

df_byYear <-  dplyr::summarize( dplyr::group_by( df , year), drivers= sum(drivers) , front=sum(front) , rear = sum(rear) ) #年ごとに集計
df.melt <- reshape2::melt( dplyr::select( df_byYear, year , drivers , front, rear)  , id.vars = "year" ) #縦に溶かす

geom_area()で描画

見た目が思い通りであれば、別にどの関数を使用しても構わない派である。
ゆえに、参考URL2で紹介されているgeom_area()を用いた方が簡易な気がするため、こっちにする。

まず普通にやる。

tsumiage1 <- ggplot( data = df.melt) + geom_area( mapping= aes( x = year , y = value , fill = variable) )
plot(tsumiage1)

f:id:ronri_rukeichi:20181214144444p:plain

多少カスタマイズする

quartzFonts(HiraMaru=quartzFont(rep("HiraMaruProN-W4",4)))  #Font指定
tsumiage2 <- ggplot( data = df.melt) 
+ geom_area( mapping= aes( x = year , y = value , fill = variable) , color = "black", alpha = 0.7) 
+ scale_fill_brewer(palette="Oranges",name="死亡内訳" , breaks = c("drivers" , "front","rear") ,labels = c("運転者" , "前部座席", "後部座席" ) )
+theme_bw(base_family = "HiraMaru") +xlab("年度") + ylab("死亡者数")

f:id:ronri_rukeichi:20181214144559p:plain

カスタマイズ上参考になるURLはこちら


Enjoy!

カーネル密度推定を君へ

R 統計学

Motivation

数値ベクトルxが与えられているとき
xの値が与えられていない区間についても、なめらかな形で確率密度関数を求めたいときがある。

そんなときに多く用いられるのが、カーネル密度推定である。
カーネル関数をK(y)を用いた以下の式で、（実際の観測値の有無にかかわらず）任意のf(X = x)の値を求めることができる。

f(x) = \frac{1}{nh}\Sigma_{i = 1 }^n K(\frac{x - x_i}{h})

しかし、Rのstats::densityは実はこの式を直接実行することはできない。
（from,to,n などの引数をいじって無理やり推定点をずらしても、内部ではapproxを用いて近傍の2つの点の補間を行なっているので、上記の式が直接適用されているわけではない）
ksパッケージのkdeも同様で、内部でapprox()を使っているので直接上記の式を計算しているわけではない。
（結構な落とし穴だと思うがわかりやすく明言されていない）

なので、カーネル密度推定を用いた統計手法を用いる時は、自分で上記式の計算関数を実装するのがよい。

ガウス関数の場合

カーネル関数としてガウス関数（標準正規密度関数）を用いるとすると、
f(x) = \frac{1}{nh}\Sigma_{i = 1 }^n K(\frac{x - x_i}{h})
f(x) = \frac{1}{nh}\Sigma_{i = 1 }^n \frac{1}{\sqrt{2\pi}}exp({-\frac{(\frac{x-x_i}{h})^2}{2}})
f(x) = \frac{1}{nh}\Sigma_{i = 1 }^n h \times \frac{1}{\sqrt{2\pi h^2}}exp({-\frac{({x-x_i})^2}{2h^2}})
f(x) = \frac{1}{n}\Sigma_{i = 1 }^n \frac{1}{\sqrt{2\pi h^2}}exp({-\frac{({x-x_i})^2}{2h^2}})

ここで\Sigmaの中を見てみると、これは平均x_i、分散がh^2の正規分布の密度関数となっている。
だから、所定の方法でband幅hさえ求めてしまえば、あとはdnorm(x,x_i, h)の平均を求めればそれが密度となる。
各ケースにweightが付与されている場合は、平均が重み付き平均になるだけだ。

海外サイトのMLでも実装例が示されている（weightには対応してない）

実装

Rで実装する。weightにも対応するようにしている。
バンド幅は、bw.SJ()とbw.nrd0()関数を利用して求めている

density_pntEst <- function(x,  at , w=NULL ,bw = "sj-ste"){
# list-wise deletion
if(is.null(w)){
w <- rep(1 , length(x))
}
df <-  data.frame(x=x , w = w)
idx <- complete.cases(df)
x <- df[ complete.cases(df) , "x"]
w <- df[ complete.cases(df) , "w"] 
bw <- tolower(bw) #小文字化

# band width
if(bw %in% c("sj", "sj-ste")){
bwval <- bw.SJ( x , method = "ste")
}else if( bw == "sj-dpi"){
bwval <- bw.SJ( x , method = "dpi")
}else if(bw == "nrd0"){
bwval <- bw.nrd0(x)
}else{
stop("band width決定手法が適切ではない")
} #if

w <-  (w * length(x) )/ sum(w) # weightを平均1に調整

at <- matrix(at , ncol = 1) 
y <- apply(at , 1, function(a,x ,bw,w ){
return( mean(  dnorm( a, x , bw)  * w )) 
} , x= x , bw = bwval,w =w)

return( list(x = x , y = y , bw = bwval , index =idx  ) )

} #function

汚いコードなのは重々承知だが、何百回も呼び出すわけじゃないからご愛嬌。

www.youtube.com

車輪の再発明としての二項ロジスティック回帰 feat. 指数分布族

R 統計学

Motivation

ちょっと複雑なモデルを開発・実装しなくてはならないので、単体テスト的な意味で車輪の再発明をする。
サービス提供の質や速度の観点からすると車輪の再発明は望ましくないが、学習者は車輪の再発明をすることを恐れてはいけない、ってどっかの偉人が言っていた気がする（言っていないかもしれない）。

指数分布族に属する分布は対数尤度関数の偏微分（＝スコア関数）が、解析的に求められることができるし、それを利用してGLM（といえる統計モデル）は高速に最適化できるよって話。
新たなモデルを思いついた時、最尤関数だけを求めるのは比較的簡単なのだが、実用性（計算スピード）を考えると勾配も求められないと厳しいのである。

参考資料

  1. 大阪大学の舟木剛教授の授業資料（応用システム工学）
  2. 高井・星野・野間（2016）『欠測データの統計科学』
  3. WEBに落ちてたプリンストン大学の講義資料 by Barbara Engelhardt（機械学習とかやってる人っぽい）

欠測データの統計科学――医学と社会科学への応用 (調査観察データ解析の実際 第1巻)

欠測データの統計科学――医学と社会科学への応用 (調査観察データ解析の実際 第1巻)

指数分布族とは

f( y| \theta ) = exp(a(y)^T b(\theta) + c(\theta) + d(y))

という形でその密度関数を表現可能な確率分布は、指数分布族（exponential family）である。
ここでa(y)は、確率変数のみで構成された関数のベクトルであり、b(\theta)は、パラメータのみで構成された関数のベクトル、c(\theta)はパラメータのみで構成された関数、d(y)は確率変数のみで構成された関数である。

正規分布をはじめ、統計学で頻出する有名な確率分布はその多くが指数分布族である。

指数分布族の便利な性質

指数分布族にはいくつか便利な性質があるが、その最たるものは対数尤度関数の微分*1を解析的に求められることである。

対数尤度関数をl(\theta) , その（偏）微分をl'(\theta)と表すこととすると、

l'(\theta) = a(y)^T b'(\theta) + c'(\theta)

と書ける。
よく知ってる形の密度関数をわざわざまどろっこしい形に書き直すことの利点はここにある。
比較的用意に勾配の関数が特定できるので、（最尤推定等の）最適化計算を準ニュートン法などで高速に行える。

他にもいくつか便利な性質はあるが割愛（詳しくは舟木先生のスライドにて）。

指数分布族の性質とGLM

さて、我々が慣れ親しんでいるgeneralized linear model（GLM, 一般化線形モデル）は、

  • 従属変数が従う確率分布 f(y | \theta )
  • 独立変数の影響をとらえる線形予測子 \mu = \beta X
  • 線形予測子と確率分布上のパラメータとの関係を表すリンク関数 \eta(\mu) = \theta

といった三つの部分からなる。
我々が知りたいのは、\beta_jの最尤推定値なので、用いたい勾配というのは \frac{\partial l}{\partial \beta_j}である。

ここで、
\frac{\partial l}{\partial \beta_j} = \frac{\partial l}{\partial \theta} \frac{\partial \theta}{\partial \mu} \frac{\partial \mu}{\partial \beta_j}

となるが、前述の指数分布族の性質を利用すれば第一項は解析的に容易に求められる。
また第三項は言うまでもなく \frac{\partial}{\partial \beta_j}\Sigma_{i=1}^{K} \beta_i x_i = x_j
第二項は（逆）リンク関数を偏微分すればよい。

かくして、GLMの勾配は容易に求められる。

指数分布としてのベルヌーイ分布

道具立ては整ったので、実装に向けて具体的に二項logistic回帰モデルについてかんがえていく。
上で書いたGLMの三要素に沿って書くと、

  • 従属変数はベルヌーイ分布にしたがう
  • 線形予測子とベルヌーイ分布のパラメータは、logit変換によってリンクされる

というのがlogit modelである。

そこでまず、ベルヌーイ分布を正規分布族として表してみると、
f(y) = \mu^y (1-\mu)^{1-y} = exp{(y log{\mu} + (1-y)log{(1-\mu)} )= exp(y log{\mu}-ylog(1-\mu) + log(1-\mu))}
となるので
a(y) = \begin{pmatrix} y & y \end{pmatrix}^T
b(\theta) = \begin{pmatrix} log{\mu} & -log(1-\mu) \end{pmatrix}^T
b'(\theta) = \begin{pmatrix} \frac{1}{\mu} & \frac{1}{1-\mu} \end{pmatrix}^T
c(\theta) = log(1-\mu)
c'(\theta) = -\frac{1}{1-\mu}

したがって、上の公式を利用して、
l'(\theta) = a(y)b'(\theta) + c'(\theta) = y\frac{1}{\mu(1-\mu)} -\frac{1}{1-\mu}
となる。ここから、
\frac{\partial l}{\partial \beta_j} = \frac{\partial l}{\partial \theta} \frac{\partial \theta}{\partial \mu} \frac{\partial \mu}{\partial \beta_j}= (y\frac{1}{\mu(1-\mu)} -\frac{1}{1-\mu})\mu(1-\mu)x_j
= (y_i - \mu_i)x_{ij}
が導かれる。これを全ケースに関して足し合わせたものを勾配とすればよい。

Rによる実装

以上勾配を求めることができたので、これで車輪の再発明をしてみる。

まず、勾配を返す関数を組む（第一引数がパラメータ、第二引数が従属変数、第三引数が独立変数）

#############------ 勾配関数 ------#############
logreg_grad <- function(prms , y , x, decreasing=F){

#- - 線形予測子の列をつくる- - #
x_mat <- cbind(1,  as.matrix(as.data.frame(x))) #切片項目のため
y <- as.numeric(y)

# N × 1の線形予測子列をつくる
eta <-  as.numeric(matrix(prms , nrow=1 )) %*% t(x_mat)
mu <- 1 /  ( 1+ exp(-1*eta) ) #逆logit 変換

#--- 勾配を計算して返す- - #
y_mu <-  matrix( rep( (y-mu) , ncol(x_mat) ) ,ncol=ncol(x_mat) ) #

grad_mat <- y_mu * x_mat
grad_vec <- apply( grad_mat , 2, sum) 

# y-mu: N × 1 , x_mat = N * K 
return(grad_vec)
} #func

つぎに、尤度関数を組む

#############------ 尤度関数 ------#############
logreg_llik <- function(prms , y , x, decreasing=F){

#- - 線形予測子の列をつくる- - #
x_mat <- cbind(1,  as.matrix(as.data.frame(x))) #切片項目のため
y <- as.numeric(y)

# N × 1の線形予測子列をつくる
eta <-  as.numeric(matrix(prms , nrow=1 )) %*% t(x_mat)
mu <- 1 /  ( 1+ exp(-1*eta) ) #逆logit 変換

# 対数尤度のデータ列をつくる
llik <-  sum( y * log(mu) + (1-y)*log(1-mu))
if(decreasing){
llik <- llik * -1 
} #if

return(llik)
} #logreg_llik(func)

これらを使ってlogistic regression（最尤推定）する関数を組む

############ 上の尤度関数、勾配関数をつかって、logistic regressionする関数 ############
logreg2 <- function(y_var ,x_var, tdf){

y2 <- tdf[[y_var]]
x2 <- tdf[, x_var]
tdf2 <-  as.data.frame(na.omit( cbind( y2, x2)) )
y <- tdf2[,1]
x <- tdf2[ , (2:ncol(tdf2))]

#- - -初期値の設定- - -#
prm0 <- rep(0 , ncol(x)+1 ) 

opt <- optim(par = prm0 , fn = logreg_llik ,gr=logreg_grad , y=y , x = x, control=list(fnscale = -1, maxit=1000,abstol=1e-10),method="BFGS")
return(opt)
} #function logreg2 . . .

#- - 実行- - -#

res1 <- logreg("y" , c("x1","x2"),test_df) #勾配を使わないver

#  計算速度
# ユーザ   システム       経過  
#     0.110      0.004      0.113 

res2 <- logreg2("y" , c("x1","x2"),test_df) #勾配を使うver

#  計算速度
#ユーザ   システム       経過  
#     0.049      0.004      0.052

勾配を用いたほうが2倍ほど早くなっている。
ちなみに収束するまでのiterationはlogreg（勾配なし）が182回、logreg2（勾配あり）が32回で、やはり勾配があったほうが早く収束する。

Enjoy!

*1:スコア関数（score function）とも呼ばれる

「泣き虫しょったんの奇跡」というか渋川清彦の奇跡（監督：豊田利晃）

映画 将棋


『泣き虫しょったんの奇跡』本予告

客がいない

出先からの帰りに見た＠池袋HUMAX。
公開初日だというのに、めっちゃ空いてて心配になるほどだった。いくら平日とはいえ。
やはり年齢層は50代以上の年配の男性が多い、という感じだ。
将棋ブーム（なのか藤井ブームなのかはさておき）による需要を見込んでの公開だとは思うが、なかなか実際のところは難しいのだろう。
せちがれぇ。

全体的な感想：長いハイライトムービー

予告を見るとまず抱くのは、キャストが豪華である、といった印象である。
妻夫木も早乙女太一も藤原竜也もRADWIMPSのアイツも出るのかよっていう。
永山絢斗に板尾創路までいるじゃん、予算大丈夫かよっていう。
なんかみんな曲者っぽい役だし、もうワクワクがとまらねぇぜ！っていう予告である。

んで、実際に見た後に抱くのは、
「あ、予告でだいたい全部終わってた」感というかダイジェストムービー感である。
いろんな人物が入れ替わり立ち変わり現れては、2-3分で退場していくっていう。
早乙女太一も妻夫木聡も藤原竜也も、全然彼らがやる必要性も蓋然性もなかった。

早乙女太一がやっているあの変な対局姿勢の奨励会員とか、絶対なんかやってくれるって期待しちゃうじゃん。でもすぐいなくなるもん、そんなんできひんやん普通。

ある人物の半生を2hで描くっていう制約、
しかも①将棋と出会った幼少期、②奨励会時代、③退会後のアマ時代、④編入試験
の4つの時期をそれなりに描かなくてはいけない、という制約があるから大変だったのはわかるんだけど、
90秒の予告を、すべて均等の濃度で2時間に引き伸ばした、という感があった。
もう少し掘り下げて書くところとそうでないところ、濃淡をつけてもよかったのでは。

ラストの対局シーンでの盛り上げ方が、これまでに出会った人の言葉を走馬灯のように流していく、という方法にしか求められなかった点もこのダイジェストムービー感を増幅している。

映画批評をする枠組みも見識もないし技術論も語れないので、
じゃあどこをどうすればよかったのよ！って言われると口をつぐむしかないのが歯痒いところではあるんだが。

将棋を題材にする難しさ：「聖の青春」との比較

豊田監督は元奨励会員であることから、「将棋の面白さを伝える」という意味で、
これまでの将棋映画の監督とは異なる期待を寄せられていた（いる）と思う。

映像化という観点だけでなく普及という点からいっても、
将棋の魅力を将棋をよく知らない人に伝えるのはむずかしい。

ある程度将棋を知っているひとにとってみれば、
たとえば去年の竜王戦第4局渡辺-羽生戦の△6六飛や今年の竜王戦予選5組決勝石田-藤井の△7七同飛成といった手は、
その一手の意味がわかってくるとともに、背筋をゾワっとさせるような興奮や人間の叡智というものを感じさせるような力をもっている。


しかし、将棋をわからない人から見ればそれは、ほかの変哲もない手と同様に、
パチっという駒音とともにひとつ駒が移動したという現象が目に映るにすぎないのである。
（そういった意味で近年の将棋界が、将棋めしや棋士のキャラクター性を全面に押し出すことで、指し手を楽しむのとは違う方向性でのファン＝「観る将」の開拓を進めているのは興行的に正しい）


だから映画をつくるにあたっては、どうにかして将棋を指している棋士の対局姿をうまく撮ることで、なんとか魅力を伝えていくしかない。
ただ、その点においても、ダイジェスト・ムービーたる「しょったんの奇跡」は弱いところがあった。
対局シーンが映し出されるシーンは数多くあるものの、ひとつの対局＝真剣勝負をじっくり見せるというシーンが皆無に等しかったのが残念であった。
結婚式で流れる新郎新婦の紹介ムービーではないが、どんなに長くても一局は2-3分で終わってしまっていた。


その点は映画としては「聖の青春」のほうがうまくやれていて、
村山-羽生戦における両棋士（というか松ケンと東出昌大）の
鬼気迫る姿をじっくり見せることによって、見る側が「これはただのボードゲームじゃない」と感じられるような真剣味を感じさせることができていた。

聖の青春 [DVD]

聖の青春 [DVD]

この棋譜はNHK杯のものであるものの、羽生と村山が和服で長時間対局するという設定は架空のものであり、この脚色自体については長年の将棋ファンからすれば賛否あるものだったと思う。


この変更点に関して、私は映画としてはアリだが、原作も素晴らしい作品なので否定派の気持ちもわからんでもない、という立場であったが、
映画「しょったん」を見た後に改めて考えて見ると、これはやはり森義隆監督（かあるいは脚本の向井康介）の名采配だったのだ、という印象をもつ。

よかったところ

上述のように映画全体に関しては若干一本調子な印象をもったが、もちろんひとつひとつのシーンや演技、セリフには素晴らしい要素もたくさんあった。

まず、やはり心動かされるセリフがいくつもあったのは良い点だと思う。小林薫扮するアマ強豪の人（名前忘れた）のセリフで、

「人生は負けたら終わりだけど将棋は負けても次の対局がある。敗北も楽しめないと将棋を愛せない」（うろ覚え）

みたいな言葉があったんだけど、なんか沁み入るものがあった。
小林さんの演技がよかったのもあるんだけど、いや人生も同じだよなぁと。敗北も愛せるようにならないとたぶんいい味出せないよなぁと。黒星を積み重ねながら歳をとっていくおっさんとしてはそう思わされた。
あとキャッチコピーにもなっているけど、「負けっぱなしじゃ、終われない」っていうのはシンプルだけど力のある言葉で、やはりあのタイミングで出てくると感動してしまう。

プロ試験にのぞむ瀬川に旧友がかける

しょったんの弱点は、勝つことに慣れてないことだよ
勝つことの喜びを恐れるなよ

という言葉も重みがある。
「勝つことの喜びを恐れない」っていうのは本当に重要な言葉で、もうちょいやれるのに自分にブレーキをかけてしまうことが多くなる齢の自分には刺さるものがあった。

そして何より白眉は、豊川七段（がモデルとなった棋士）を演じた渋川清彦さんのすばらしい演技である。
何を言っているがわからないと思うが本人よりもマンモス感あった。なんというか濃縮還元100%マンモスである。
時折出るダジャレ（「こんばんワイン〜♪」）だけでなく、カラオケではしゃいでる姿や、後輩を暖かく見守る姿、そのすべてが我々の「豊川孝弘」像の理念型みたいなものをうまく表現できていて、個人的にはこの映画のMVPであった。


ってことでおつかれマンモス！


www.nicovideo.jp

Stataで出力したregression tableをcsv形式で出力 by estoutパッケージ

Stata 統計学

備忘用。年を重ねるごとに物忘れが激しくなってしまうな。

Motivation

Stataで回帰分析を行なった結果をExcelにうつすときのやりかた（これがbestかはわからん）。

Texで出力するのがスタンダードでスマートなのはわかっているが、
実際の仕事だと（例えば出版社に回帰分析の元表を提出するときなど）Excelベースの管理が一般的という事情もある。

あとMacユーザだと、Excelでregression tableを作っておけば、
その部分を選択コピーしてからプレビューを開いて、command + Nで比較的綺麗に画像形式（pdfファイル）に変換できるというのがコスパがよくて重宝する。

そうなると（Texの習得に費やす労力も考えると）やっぱりExcelにベタ貼りできるような形で整形・出力してくれるようにしたくなるのだ。

"estout"パッケージ。

ふだんあまりStataを使わないこともあって調べ方が悪いのかもしれないが、 標準で実装されている機能であまり御誂え向きなものはなかった。 今の所、Excelに貼り付けられそうな形での出力になじんだパッケージとしてestoutを利用するのが、 それなりに簡単で、実用にたえうる手段であるように認識している。
※outreg2というのもあるらしいけどよくわからない。

以下、参考URL

install

Stataのコンソールで以下を実行。 

ssc install estout, replace

eststo関数によりtableを保存

たとえば普通の重回帰分析を2本推定して、それを並べた表を出力したい時。
まず、以下のように、回帰分析の結果を保存する。

*Model1: 従属変数yをx1,x2で回帰 
quietly reg y x1 x2 

*mdl1という名前をつけて保存 
eststo mdl1 

*Model2 : 従属変数yをx1,x2, x3で回帰 
quietly reg y x1 x2 x3

*mdl2という名前で保存する
eststo mdl2

esttab関数により、csvファイルに出力。

上の手続きで、回帰分析の推定結果を保存したのちに、
esttab関数のusingオプションを使うことで、csvファイルとして出力できる。

以下は、上で保存した結果を

  • デスクトップにcsv形式で「tbl1.csv」という名前で保存。
  • 統計量としては調整済R自乗、対数尤度、AIC, BICを出力
  • 標準誤差は丸括弧（parentheses）でなく角括弧（brackets）にする。
  • Model名を「Model1」「Model2」と指定
  • 回帰係数、標準誤差ともに、小数点以下4桁まで表示する。

といった形式で出力するときのプログラムである。

esttab mdl1 mdl2 using "/Users/Ronri_Rukeichi/Desktop/tbl1.csv" , /*
*/replace b(4) se(4)  ar2  scalar(ll) aic bic brackets mtitles("Model1"  "Model2")

以下補足

  • scalarsの中に指定できる統計量については、このページのExample3を参照。例えばパラメータ数も表示したい時はscalars(ll rank)と指定する。
  • 詳しい構文はesttabのhelpを参照。


場合によってはspost9アドインが求められるが、そのときは以下のプログラムを実行する

net install spost9_ado

対応分析（correspondence analysis, CA）では一体何を「分解」しているのか

R 統計学

Introduction

小野滋さんの「読書日記」が久しぶりに更新を再開し、ひとり小躍りするGW中盤戦...

前々回の記事で、SVDとはなんぞや、について解説した。
ronri-rukeichi.hatenablog.com


この記事のなかで、対応分析も多重対応分析も結局はSVDだよーってちらっとふれたんだけど、とりあえず今回は対応分析がSVDであることを示さんとす。

目の前の仕事が煮詰まっていることからの現実逃避もかねてね...

[参考資料]

  1. Clausen（藤本訳）「対応分析入門」
  2. Greenacre, 1984 , Theory and Applications of Correspondence Analysis

対応分析入門 原理から応用まで 解説◆Rで検算しながら理解する

対応分析入門 原理から応用まで 解説◆Rで検算しながら理解する

前提知識

とりあえず前々回記事からの抜粋で復習すると、

行列AのSVDが得られているとき、降順に並んだj \ ( j <= k ) 個の特異値と、それに対応する左右の特異ベクトルからなる
A_{j^*} = \Sigma_{i = 1}^j \lambda_j u_i v_i は、
すべてのランクjの行列X_{j^*}のなかで、距離|| \boldsymbol{A - X_{j^*}}||_F を最小化するような行列となる。


すなわち、この意味で A_{j^*} は、Aの「ランクjの最小二乗近似」であるといえる。

フロベニウスノルムを最小化するような近似行列を得るような特異値& 左右特異ベクトルのセットを得るための方法として、SVDはあった。
だから、上の説明でいうところの、Aに分解の対象としたい、データ全体の適切な散らばりの表現の行列を指定しさえすれば、色々な応用がきく。

実際に、Greenacre（1984, pp.348-9）においては、Aに様々な行列を選ぶことで、さまざまな手法がSVDとしてとらえられることを一覧化している。
（ただし、Aだけでなく、距離計算のときの重み付けを指定する正定値対称行列\Omega , \Phiの指定も考える必要がある。これに関しては今まで説明していない一般化SVDという枠組みが関係してくるが、説明はおいおい...）

SVDとしての対応分析

いきなり答えから載せちゃうと、対応分析は、
行列 D_r^{-1/2} ( P - r c^T ) D_c^{-1/2}を特異値分解している。

  • D_rは、対角成分に各行の周辺相対度数をとる対角行列、D_cは、対角成分に各列の周辺相対度数をとる対角行列、
  • rは各行の周辺相対度数をとるベクトルで、cは各列の周辺相対度数をとるベクトルである。
  • Pは「同時確率行列」と呼ばれるものだが、要するに相対度数の行列で、クロス表を総度数で割ったものである。

アルファベットがゴツゴツとした感じで並んでて感じ悪い行列だが、
この意味は（初等的なカテゴリカル・データの扱い方を学んだ経験のある者からすれば少なくとも）比較的明確である。

まず、 D_r^{-1/2}, \ D_c^{-1/2}はそれぞれ対角行列であるが、対角行列を左からかける場合対角成分の値を各行にかけたものが、
対角行列を右からかける場合対角成分の値を各列にかけたものが、得られる。
よって、i行j列目の分母は、完全独立が成立している場合の相対期待度数の平方根となる。

つぎに、 P - r c^Tだが、これは相対度数から完全独立が成立している場合の相対期待度数を引いたものである。これが分子となる。

よって、行列 D_r^{-1/2} ( P - r c^T ) D_c^{-1/2}の各要素というのは
（相対度数- 完全独立が成立している場合の相対期待度数）／（完全独立が成立している場合の相対期待度数の平方根）
となっている。

さて、ここで D_r^{-1/2} ( P - r c^T ) D_c^{-1/2}にSVDを適用するということは、
フロベニウスノルム || D_r^{-1/2} ( P - r c^T ) D_c^{-1/2} || への最小近似を提供するような分解を得るということである。

フロベニウスノルムの定義より、 || D_r^{-1/2} ( P - r c^T ) D_c^{-1/2} ||というのは、
すべての行列の要素に関して、その要素の自乗を足し合わせたものであるので、
要するに、

（相対度数- 完全独立が成立している場合の相対期待度数）の二乗／（完全独立が成立している場合の相対期待度数） の総和

となる。

これはどこかで見たことがないだろうか....
そう、これはカイ自乗値の定義式そのものである

つまり、

対応分析というのは、分析対象となる二重クロス表から算出されるカイ自乗値への最小近似を実現するような特異値と特異ベクトルを、SVDによって得ているだけの処理である

といえるのである。

SVD後のこまかい処理の話:対応分析から得られるスコアはどうやって得られるか

そしてこっからは対応分析に関するSVD「後」の話になる。
A = D_r^{-1/2} ( P - r c^T ) D_c^{-1/2}とし、このSVDから A = USV^Tという分解が得られる。

ここから、 座標（行・列）のスコアに関して、重み付きの自乗和が固有値と一致するような形で特異ベクトルを変換したい。
すなわち、行スコアをF_r 、列スコアをF_cとすれば、
(D_r^{1/2} F_r)^T D_r^{1/2}F_r= F_r^T D_r F_r = S^2
(D_c^{1/2} F_c)^T D_c^{1/2}F_c= F_c^T D_c F_c = S^2
が成り立つようにしたい。

SVDの性質より、 U^T U = V^T V = I が成り立っているので、
F_r = D_r^{-1/2} U S , F_c = D_c^{-1/2} V S という変換を施せば、

F_r^T D_r F_r = S^T U^T D_r^{-1/2} D_r D_r^{-1/2} U S = S^T U^T U S = S^2
F_c^T D_c F_c = S^T V^T D_c^{-1/2} D_c D_c^{-1/2} V S = S^T V^T V S = S^2
という形で上制約を満たすスコアが得られる。

これが対応分析で得られる行／列のスコア（座標）である。

Conclusion

ってことで、対応分析という手法は、
フロベニウスノルムがカイ自乗値となるような行列に対してSVDを適用したのちに、
スコアの重み付き自乗和が固有値（=特異値の自乗）と一致するように左右の特異ベクトルを変換するだけの簡単なお仕事だとわかった。

ちなみにMCAもほとんど同じようなロジックであるが、ちょっと疲れたので気が向いたらまたまとめよう


Def Tech - Catch The Wave

そんじゃ、Enjoy!