『AICをぶっ壊す』とAICの復習の学習記録

【学習動機】

分析に用いているモデルの評価にAICが使えるのか不安になったので, AICについて学び直すことにした. 小西・北川『情報量規準』でAICの導出を追い直しつつ, Rで数値実験しているものを探していたら, ノンパラでよく見かける竹澤先生の書いた『AICをぶっ壊す』というものに出会った.
http://cse.naro.affrc.go.jp/takezawa/AIC%E3%82%92%E3%81%B6%E3%81%A3%E5%A3%8A%E3%81%99.pdf
このドキュメントを読んで, コードがもっと速くできると思ったので, コードを書き換えつつ, AICについて学び直す.

理論面はこちらで学ぶ. 記法や呼称もこちらに準拠する.

小西・北川『情報量規準 (シリーズ・予測と発見の科学)』

www.amazon.co.jp

統計学の入門を一通り終えた人には, この本をオススメする. とくに, 3章まで読むことをオススメ. というか他の章はまだ読めていないので, どんなものなのか分からない.

 

【学習内容】

AICの復習

最初に, カルバック・ライブラー情報量(以下, KL情報量, 相対エントロピーともいう. 符号を逆にするとボルツマンのエントロピーという.)からAIC(赤池情報量規準)の導出への流れを高速で確認する. ところどころ飛躍が激しいところもあるが, 詳細はテキストで確認されたい.

以下, 確率変数\(X\)の真の分布を\(G(X)\), あてはめるモデルの分布を\(F(X)\)とする. 統計学では, データというものは真の分布から生成されてくるものだという大前提を置く. この大前提は頻度論の枠組みだろうとベイズの枠組みだろうと変わらない. そして, その真の分布がどのようなものかは分からないが, 真の分布に近似したものをモデリングという営みによって作成していく. モデルをつくったとして, そのモデルがどれくらい良いものなのかを知りたいというのは自然な気持である. そのモデルがどれくらい真の分布に近いのかを測る尺度の1つがKL情報量である.

確率分布間の近さのようなもの(数学でよく出てくる距離の公理は満たしていない)として定義されるKL情報量は\[I(g;f)=E_G[log(\frac{g(X)}{f(X)})]=E_G[log(g(x))]-E_G[log(f(X))]\]と分解でき, 第1項は真の分布に依存し, モデルに依存しない. 真の分布\(G(X)\)の負のエントロピーともいう. 第2項の\(E_G[log(f(X))]\)は平均対数尤度と呼ばれる. これはモデルの分布\(F(Z)\)で定義される汎化損失の符号が逆になったものでもある.

KL情報量の符号を逆転させたボルツマンのエントロピーは, モデルから得られたサンプルの相対度数が真の分布と一致する確率の対数に比例する量として解釈できる.

また, KL情報量は\(G(X)\)\(F(X)\)でシミュレートしたとき, \(G(X)\)からの乱数のように見える確率の減少速度として捉えられる話やKL情報量が自然に出てくる話は黒木玄先生のアップしている資料や過去のツイートから学ぶことができる. まだしっかり理解できていないので, ここでは説明しない.
『Kullback-Leibler情報量とSanovの定理』
http://www.math.tohoku.ac.jp/~kuroki/LaTeX/20160616KullbackLeibler.pdf
(ちなみに, KL情報量をテイラー展開すると二次の項にフィッシャー情報量が出てくる.)

したがって, KL情報量を最小化すれば, 真の分布に近い分布がモデルから得られると考えられる. 上に書いたように第1項は真の分布にのみ依存するので, モデルの変化に対して不変である. このため, 第2項の平均対数尤度を最大化すれば(汎化損失を最小化すれば), 良いモデルが得られると考えることができる.

さて, 平均対数尤度を最大化したいわけだが, 平均対数尤度は普通未知である(つまり, 人工データでも作らない限り分からない)真の分布に依存しているので, この量を上手く推定したいというモチベーションが生まれる. 自然に思いつくものとして, 対数尤度を計算してサンプルサイズで割ったものは平均対数尤度の一致推定量となる. (最尤推定法は自然な推定法なのは, この事実があるためでもある.). つまり対数尤度は, \(nE_G[log(f(Z|\hat\theta)]\)の推定量となる.

しかし, この対数尤度による推定量はバイアスを持ち, \[b(G)=E_{G(x_n)}[log(f(X_n|\hat\theta(X_n)))-nE_{G(z)}[log(f(Z|\hat\theta(X_n))]]\] と表され, この大小はパラメータの次元に依存している. このバイアスは3つの項に分解され, 頑張って計算すると漸近的に\[b(G)=tr(I(\theta_0)J(\theta_0)^{-1})\]となる. ただし, \[I(\theta)=\int{g(x)\frac{\partial log(f(x|\theta))}{\partial \theta}\frac{\partial log(f(x|\theta))}{\partial \theta'}}dx\], \[J(\theta)=-\int{g(x)\frac{\partial^2 log(f(x|\theta))}{\partial \theta \partial \theta'}}dx\]である.

想定しているパラメトリックなモデルの族の中に\(g(x)\)が含まれている場合は, \(I(\theta_0)=J(\theta_0)\)となる.

このIとJをデータから推定する一致推定量\(\hat I\), \(\hat J\)として, バイアスの推定量\(\hat b = tr(\hat I \hat J^{-1})\)が得られる.

モデルの対数尤度をバイアスで補正した量が平均対数尤度を推定したものとなっている. これに-2を乗じたものが情報量規準TIC(竹内情報量規準)と呼ばれる.

f:id:moratoriamuo271:20190828011213p:plain

 

上に書いたように, 真のモデルがモデル族に含まれる場合, IとJが等しくなり, バイアスが\(p\), つまりモデルの自由パラメータ数となる. このバイアスで補正して, 平均対数尤度を推定したものを-2倍したものが, 情報量規準AIC(赤池情報量規準)である.

 

f:id:moratoriamuo271:20190828011249p:plain


-2倍しているのはカイ二乗統計量との関係でスケールを合わせているだけである. KL情報量に戻ってみると, 第二項の平均対数尤度が大きければよく, つまり平均対数尤度の符号を反転させたものを小さくすればよい. ゆえにAIC, またはTICが小さければ小さいほど, モデルは真の分布に近いことと解釈できる. TICに比べてAICは簡単に計算できて便利だ. ただ, 導出に使われる仮定があってその限界を持っていることや, AICもTICも確率変数なので不確実性な揺れを持っていることを忘れてはならない.

数値実験

AICの復習を終えたところで本題に入る.

AICをぶっ壊す』というドキュメントでは, 真の分布が確率\((cr1,1-cr1)\)\(N(m_1,sd_1^2)\)\(N(m_2,sd_2^2)\)の混合正規分布であるが, 想定するモデルが正規分布となっているケースで, 平均対数尤度と対数尤度のバイアスを数値計算し, それとAICのバイアス部分, つまり正規分布の自由パラメータの個数である2との違いを見ている. ドキュメントでは混合正規分布の他に一様分布と対数正規分布で計算しているが, ここでは混合正規分布のみ扱う.

以下のように関数を書きかえて定義する. 引数はサンプルサイズと真の分布のパラメータ5つをとる.

rewrite <- function(nd, cr1, m1,sd1,m2,sd2){
  set.seed(20190824)
  # generate data (nd : sample size, cr1 : miture proportions)
  # 確率(cr1,1-cr1)でN(m1,sd1^2)とN(m2,sd2^2)から生成する混合正規乱数生成
  gdata1 <- function(nd, cr1){
    n1 <- sum(runif(nd) < cr1)
    c(rnorm(n1,m1,sd1),rnorm(nd-n1,m2,sd2))
    #c(dqrng::dqrnorm(n1,m1,sd1), dqrng::dqrnorm(nd-n1,m2,sd2))   
    #{dqrng}は乱数生成高速化パッケージだが、この程度だと関数呼び出しあたりが悪いのかむしろ遅かった。
  }
  lf1h <- numeric(3000)
  lf1f <- numeric(3000)
  for(iii in 1:10){
    for(kkk in 1:3000){
      xxt <- gdata1(nd, cr1)
      mean1 <- mean(xxt)
      sdxxt <- sqrt(sum((xxt-mean1)^2)/nd)
      lf1h[kkk] <- nd * log(1/(sqrt(2*pi)*sd1)) - 0.5 * nd
      lf1 <- numeric(20)
      for(kk in 1:20){
        xxf <- gdata1(nd, cr1)
        lf1[kk] <- nd * log(1/(sqrt(2*pi)*sd1)) - 0.5 * sum((xxf-mean1)^2)/sdxxt^2
      }
      lf1f[kkk] <- mean(lf1)
    }
    print(mean(lf1h - lf1f))
  }
}

混合正規乱数の最適化についてはこのページを参考にした.
http://www.okadajp.org/RWiki/?R%E3%82%B3%E3%83%BC%E3%83%89%E3%81%AE%E6%9C%80%E9%81%A9%E5%8C%96%E4%BE%8B%EF%BC%9A%E6%B7%B7%E5%90%88%E6%AD%A3%E8%A6%8F%E4%B9%B1%E6%95%B0%E3%81%AE%E7%99%BA%E7%94%9F%E3%82%B3%E3%83%BC%E3%83%89

これにより, ドキュメントと同様の結果が出るコードが走る. 混合正規分布の割合パラメータを0, 0.25, 0.5, 0.75, 1のケースで確認する.

rewrite(nd=50, cr1=0.0, m1=2,sd1=3,m2=2,sd2=7)
## [1] 2.181348
## [1] 2.236288
## [1] 2.09223
## [1] 2.095036
## [1] 2.510326
## [1] 2.165027
## [1] 2.051087
## [1] 2.084752
## [1] 2.171209
## [1] 2.248626
rewrite(nd=50, cr1=0.25, m1=2,sd1=3,m2=2,sd2=7)
## [1] 2.530909
## [1] 2.58092
## [1] 2.373811
## [1] 2.473825
## [1] 2.734282
## [1] 2.486
## [1] 2.389523
## [1] 2.442735
## [1] 2.366403
## [1] 2.529645
rewrite(nd=50, cr1=0.5, m1=2,sd1=3,m2=2,sd2=7)
## [1] 3.094842
## [1] 3.027269
## [1] 2.839038
## [1] 2.985219
## [1] 2.986768
## [1] 2.98658
## [1] 2.938266
## [1] 2.913371
## [1] 2.812462
## [1] 3.133202
rewrite(nd=50, cr1=0.75, m1=2,sd1=3,m2=2,sd2=7)
## [1] 3.630774
## [1] 3.527391
## [1] 3.205531
## [1] 3.396918
## [1] 3.753438
## [1] 3.520168
## [1] 3.45365
## [1] 3.335673
## [1] 3.162852
## [1] 3.696705
rewrite(nd=50, cr1=1.0, m1=2,sd1=3,m2=2,sd2=7)
## [1] 2.181348
## [1] 2.236288
## [1] 2.09223
## [1] 2.095036
## [1] 2.510326
## [1] 2.165027
## [1] 2.051087
## [1] 2.084752
## [1] 2.171209
## [1] 2.248626

混合正規分布の割合パラメータが0と1のケースは普通の正規分布になっていて, 真の分布が想定しているモデル族に含まれるのでAICのバイアス部分の2で上手くいってそうなことが分かる. 一方でそれ以外のケースではAICのバイアス部分の2よりも大きくなっていて, AICのバイアス部分では過少推定しているであろうことが見て取れる.

せっかく書き換えをしたので速度比較をするため, ドキュメントのコードからコメントを抜いてそのまま書き写した関数も用意した.

allcopy <- function(){
  gdata1 <- function(iseed, nd, cr1){
    cr1 <- 0.5
    xxa <- NULL
    set.seed(iseed)
    for(ii in 1:nd){
      xx1 <- rnorm(n=1, mean=2, sd=3)
      xx2 <- rnorm(n=1, mean=2, sd=7)
      uni1 <- runif(n=1, min=0, max=1)
      xxa[ii] <- xx1
      if(uni1 > cr1) xxa[ii] <- xx2
    }
    return(xxa)
  }
  nd <- 50
  lf1h <- NULL
  lf1f <- NULL
  for(iii in 1:10){
    for(kkk in 1:3000){
      iseed <- kkk + iii*912
      xxt <- gdata1(iseed, nd, cr1)
      mean1 <- mean(xxt)
      sd1 <- sqrt(sum((xxt-mean1)^2)/nd)
      lf1h[kkk] <- nd * log(1/(sqrt(2*pi)*sd1)) - 0.5 * nd
      lf1 <- NULL
      for(kk in 1:20){
        iseed <- kkk + kk * 471 + iii * 412
        xxf <- gdata1(iseed, nd, cr1)
        lf1[kk] <- nd * log(1/(sqrt(2*pi)*sd1)) - 0.5 * sum((xxf-mean1)^2)/sd1/sd1
      }
      lf1f[kkk] <- mean(lf1)
    }
    print(mean(lf1h - lf1f))
  }
}

速度比較をしてみると

system.time(allcopy())
## [1] 2.865864
## [1] 2.775913
## [1] 2.744208
## [1] 2.960808
## [1] 2.986057
## [1] 3.040051
## [1] 2.98391
## [1] 2.904098
## [1] 2.917919
## [1] 2.878093
##    user  system elapsed 
##  218.67    0.01  218.95
system.time(rewrite(nd=50, cr1=0.5, m1=2,sd1=3,m2=2,sd2=7))
## [1] 3.094842
## [1] 3.027269
## [1] 2.839038
## [1] 2.985219
## [1] 2.986768
## [1] 2.98658
## [1] 2.938266
## [1] 2.913371
## [1] 2.812462
## [1] 3.133202
##    user  system elapsed 
##   11.53    0.02   11.54

だいぶ速くなったように思える. 並列化や関数呼び出しの省略を上手く使えば, もっと速くなりそうな気もするが, コードの分かりやすさ(と面倒くささ)から今回はこれでとどめておく.

さて, 速度は今回の話の本質ではない. AICがいい感じに平均対数尤度の-2倍を推定できていないことが数値実験から分かったというのが今回の本質である.

AICの定義するバイアス部分では平均対数尤度を対数尤度で推定したときのバイアスの補正が上手くいかない例が今回の数値実験例であった. 導出のところで確認すべきところだったのだが, 上手くいかないケースは, 小標本の場合, p/nが大きい場合, 真の分布が含まれていない場合, 大標本の場合とがある. 加えて, 非正則の場合もある.

庄野『モデル選択手法の水産資源解析への応用-情報量規準とステップワイズ検定の取り扱い-』
https://www.jstage.jst.go.jp/article/jjb/27/1/27_1_55/_pdf
は, 小標本の場合, 大標本の場合,真の分布が含まれていない場合の数値実験をして, 各情報量基準でモデル選択をして正解しているか否かを見ている. これを見てみると, AICはそこそこ頑張っているように見える.

バイアスの推定が上手くいっていない場合にどの程度上手くいかないかが分からないが, 多少仮定から外れていてもモデル選択は上手く行えるのかもしれない. ただ, 渡辺澄夫先生曰く, 複雑なモデルつまり,構造をもつモデル(混合正規分布、神経回路網、ベイズネットワーク、隠れマルコフモデル、 縮小ランク回帰、深層学習などの隠れた変数・階層構造・モジュール構造を持つモデル)ではAICなどではモデル評価が上手くできないらしい. 渡辺澄夫先生のHP『広く使える情報量規準(WAIC)』
http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/waic2011.html

今回の数値実験混合正規分布もそれにあたり, このような場合はWAICを使えば良いらしい. 今回の実験は, 平均が同じ2つの正規分布の混合分布であったが, これダメなタイプなのかいまいち分かっていない. 今度, 単峰と複峰で乱数を発生させてAICとWAICの挙動を実験してみたい. いつかそんな記事を書く.

 

 

ドキュメントのタイトルは『AICをぶっ壊す』という過激なものであったが, ぶっ壊すというよりはAICの数理的な構造を理解して, その妥当性と限界を知ることが重要なように思われる. これらを理解した先にWAICを理解して使えるようにしていけばいいような気がしてならない.

真の分布が正則で学習モデルに含まれいるなら, WAICはAICに, 正則で学習モデルに含まれていないならTICに, 確率変数として等価だそうだ。
http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/cross-val.html

【学習予定】

WAICは存在は知っているがまだ全然理解できていない. 数値実験したり, 数式を読み解きながら, ゆっくりでもいいから理解していきたい.

渡辺『ベイズ統計の理論と方法』

https://www.amazon.co.jp/dp/4339024627/ref=cm_sw_r_tw_dp_U_x_dPuzDbACQ08KB

 

うなぎのマスカラ

 関東地方で梅雨が明けた日くらいに、『天気の子』を観た。観終えた直後の個人的評価は曇り気味だった。ところが、見えないものが見えている人の書いた映画の感想のような何かの記事を読んで、自分にも、見えなかったものが見えてきたおかげで、面白い作品だったと結論付けられた。この記事である。

cr.hatenablog.com


見える人には見えるし、見えない人には見えない何かで、作品の奥行が変わってしまうものだと再度認識した。

 美しい女性が世界を大きく動かすというのは、昔からある話だ(美しいは、時代や場所、文化、その他諸々によって変化する。ここで書かれる美しいは私主観の形容詞ではなく、私主観の客観的な形容詞である。また、因果関係、事実関係についても、以下あやしいところがあるかもしれない。)。歴史の教科書でよく出てくるのは、百年戦争の戦況を覆したジャンヌ・ダルクだろうか。間接的かもしれないが、楊貴妃なんかもそうだ。傾城の美女なんて言葉もある。そんな話は前近代的だなんて言う人もいるかもしれないが、ここ数年で、2人の美人の死が日本を動かした。広告代理店の電通の労働環境で過労の末、自殺した美人な女性がいて、その後、社会全体の労働のルールが大きく変わった。軽井沢のスキー場へ向かっていたが墜落してしまった深夜バスに美人が乗り合わせていたため、大きく取り上げられることとなり、深夜バスの運転のルールが変わったし、料金も動いたように思う。政策形成の俎上に載るには、いくつかのパターンがあるが、その一つが、世間で騒がれるというものだ。美人は世間に大きなインパクトを呼び起こす。そのへんの中年男性じゃ同じ衝撃を起こせない。価値が平等なんてのは空想的な建前に過ぎないのだ。現代になろうと、政治は未だに「まつりごと」で、巫女の祈祷や白羽の矢の人身御供の時代(本当にあったかどうかは知らないが)となんら変わってやいないのだ。

 これが変わらぬ性質だとするならば、それを上手く利用してみたいと思うのも人間の性だろう。利用方法として一つ思いつくのは、バーチャルな美人を作り上げることだ。今の画像・映像の合成技術や音声技術をもってすれば、かなりリアルな、でも実際には存在しない美人を作り上げることは可能ではないかと思う。ただし、その存在がバーチャルなものであることは世間には知られないようにして、存在を世界に信じ込ませる。世界を騙しきったうえで、その存在を生贄に捧げれば、世界を変えることができるかもしれない。この方法ならば、感情や労力を無視すれば、現実では誰も死なずに世界を変えることができるのだ 。

 こんなくだらない話を考えながら、戦争捕虜が収容所で架空の美少女を作って、まるで本当に存在するかのように扱っていたことで、精神を健康に保って捕虜生活を耐えきったなんていう、昔読んだコピペ話を思い出した。実話なのか気になって少しググってみたところ、どうやら実話ではなさそうだ。

togetter.com


実話かどうかはともかく、人は(他の生き物については分からないとこだが)、ないものをあるものとして捉えられる。見えないけど見えるものが見えてくる。もちろん見えるものは人によりけりだ。六畳間をろくじょうかんと見える人だっているように。

大正義と美少女

 

 大学という場所を離れてからというもの、大学というものを相対化して見るようになった。違う大学に通った人や大学に通わなかった人に会う機会が増えたからだろう。高校を出てから、浪人を含めると7年間、一体、何をしていたのかを振り返って、大学教育について少し書き綴ってみる。

 大学に入ると、初年度くらいに、物事を見たままそのまま簡単には信じることができないものだということを、教養の講義か何かで学ぶ。それは、目の錯覚であったり、認知に関するバイアスであったり、何かしらのパラドックスであったり、誰も理解できないように小難しく表現しているけど中身は何もなかったよっていう事件だったりする。某大学の教養英語読本が書いていることは大体そのようなことだったように思う。トトロ見たもんといくら感情的に激しく言っても仕方ないのだ。

 教養を学んだか学んでいないかくらいで、簡単に信じられなくなった物事をいかに疑っていくかを学んでいくことになる。過去の文献等の調べ方であるとか、論理の展開の仕方だとか、疑ったことの表現の仕方のお作法だとか、そういったことを徐々に学んでいく。それらに関して、あらゆる分野に精通することはできないので、専攻する分野の全体を知って、その分野の一部に深く知っていくという形にはなる。

 大学には一つしか通っていないので、普遍的に上で書いたことが言えるのかは知らないし、他学部や他学科についても同様なのかは知らないが(最近、少し気になっているのは、モデルという考え方を身につけるのはどの学部なのかということだ。モデルについては今後いつか何か書く。)、ある程度まともに学んでいたのなら、ある程度一般的に正しいのではないかと信じている。大学院やその先については、君のその眼で確かめてくれといった感じで、ここには書かない。

 まぁ、だからといって、大学に通ったからといって、突然すごい能力が目覚めるわけでもない。「10時間で学べる」なんて表題掲げた書籍が本屋の棚に並んでいるが、実際学んだのがそれくらいの内容で過ごした人もいるだろう。そもそも、最低限の基礎体力もなく大学に通ってしまう人さえ存在する(話は少し逸れるが、予備校業界で6年間アルバイトして、東大入試を見ていたが、東大入試はよくできていて、基礎体力をもって学習した基礎知識を持ったうえで、正しく読めて正しく書ければ、合格できるようになっていると思う。理系については知らないが。このような話はいつか詳しく書くかもしれない。)。そういうわけなので、大学に行ったから偉いとかそういうわけではないが、教育機関として十分な教育機会は与えられていたように思う。

 よく名前を見かける人の書いたことがけっして正しいとは限らないし、声が大きい人が言っていることはよく聞こえるだけなんてことだってあるし、見た目だけ整えてやつの中身が全くスカスカだなんてことだってある。正確な情報にたどり着けるだけのぐぐらびりてぃ(googlability)がないと、不適切で作為的な情報にぐーぐらびてぃ(googravity)によって引っ張られてしまうことだってある。情報が非対称な労働市場も、有能な労働者が半数を下回れば、崩壊するレモン市場の出来上がりであり、陪審定理も、正答率50%を下回ればどうしようもない判別器の出来上がりである。教育機会の価値に気付けていさえすれば。あぁ、教育は重要だ。あぁ、教育は重要だ。

Bradley-Terryモデルの基本の学習記録

【学習動機】

とあるスポーツをやっていて, 過去に試合に出るメンバーを選ばないといけない役割を背負ったことがあった. その際に, 各メンバーの強さを数値化できればいいのにと思ったことがあった. 他にも連携による強さが数値化できたらいいなと思ったことがあった. そんな経験があった中で, Bradley-Terryモデルという、勝敗データから強さを推定するモデルの存在を知った. そんなわけで, Rの{BradleyTerry2}パッケージを用いて, 作成したシミュレーションデータでBradley-Terryモデルの基本を学習する.

【学習内容】

 

speakerdeck.com

 

 以下では, この発表スライドの分析に使用したコードを記載していく.

ライブラリ読み込み.

library(tidyverse)
library(ggplot2)
library(BradleyTerry2)
library(qvcalc)    #擬似標準誤差計算
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
## -- Attaching packages ------------------------------------------------------------------- tidyverse 1.2.1 --
## √ ggplot2 3.1.1     √ purrr   0.3.2
## √ tibble  2.1.3     √ dplyr   0.8.1
## √ tidyr   0.8.3     √ stringr 1.4.0
## √ readr   1.3.1     √ forcats 0.4.0
## -- Conflicts ---------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

自作の対戦データの読み込みとその表示.

dat <- read_csv("winlose1.csv")
## Parsed with column specification:
## cols(
##   player1 = col_character(),
##   player2 = col_character(),
##   win1 = col_double(),
##   win2 = col_double()
## )
dat
## # A tibble: 15 x 4
##    player1 player2  win1  win2
##    <chr>   <chr>   <dbl> <dbl>
##  1 どら    のび        7    13
##  2 どら    しず        5    12
##  3 どら    たけ        9     6
##  4 どら    すね        8     4
##  5 どら    でき        1     5
##  6 のび    しず        3    17
##  7 のび    たけ        4    14
##  8 のび    すね        6     9
##  9 のび    でき        1     7
## 10 しず    たけ       10     3
## 11 しず    すね        9     2
## 12 しず    でき        6     1
## 13 たけ    すね       17     2
## 14 たけ    でき        2     1
## 15 すね    でき        1     4

player1とplayer2の対戦成績でwin1がplayer1の勝利数, win2がplayer2の勝利数を記録している.

データ構造の確認.

str(dat)
## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 15 obs. of  4 variables:
##  $ player1: chr  "どら" "どら" "どら" "どら" ...
##  $ player2: chr  "のび" "しず" "たけ" "すね" ...
##  $ win1   : num  7 5 9 8 1 3 4 6 1 10 ...
##  $ win2   : num  13 12 6 4 5 17 14 9 7 3 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   player1 = col_character(),
##   ..   player2 = col_character(),
##   ..   win1 = col_double(),
##   ..   win2 = col_double()
##   .. )

{BradleyTerry2}パッケージのBTm関数は比較対象(ここではplayer)因子型しか受け入れてくれないので因子型に変換しておく. playerのところが同じレベルの因子型になっていないといけないので, そこも注意.

playerset <- c(dat$player1,dat$player2) %>% unique()
dat$player1 <- dat$player1 %>% factor(levels = playerset)
dat$player2 <- dat$player2 %>% factor(levels = playerset)

まずは, 勝率を確認する.

dat4count <- dat %>% mutate(nfights = win1 + win2)
datlong <- rbind.data.frame(dat4count %>% select(player1,win1,nfights) %>% 
                              rename(win = win1, player = player1),
                            dat4count %>% select(player2,win2,nfights) %>%
                              rename(win = win2, player = player2))

datlong %>% group_by(player) %>% 
  summarise(totalwin = sum(win), totalfights = sum(nfights)) %>% 
  mutate(winrate = totalwin/totalfights, rank = rank(-winrate)) %>% arrange(rank)
## # A tibble: 6 x 5
##   player totalwin totalfights winrate  rank
##   <fct>     <dbl>       <dbl>   <dbl> <dbl>
## 1 しず         54          68   0.794     1
## 2 でき         18          29   0.621     2
## 3 たけ         42          68   0.618     3
## 4 どら         30          70   0.429     4
## 5 のび         27          81   0.333     5
## 6 すね         18          62   0.290     6

すねを基準点としてモデルをあてはめる. すねを0として戦闘力パラメータを推定する.

res <- BTm(cbind(win1,win2), player1, player2, ~ player, id = "player", data = dat, refcat = "すね")
res
## Bradley Terry model fit by glm.fit 
## 
## Call:  BTm(outcome = cbind(win1, win2), player1 = player1, player2 = player2, 
##     formula = ~player, id = "player", refcat = "すね", data = dat)
## 
## Coefficients:
## playerどら  playerのび  playerしず  playerたけ  playerでき  
##     0.5662      0.2364      1.9715      1.1319      1.3129  
## 
## Degrees of Freedom: 15 Total (i.e. Null);  10 Residual
## Null Deviance:       60.62 
## Residual Deviance: 17.27     AIC: 65.41

パラメータの区間を擬似標準誤差を使って求める.

res.qv <- qvcalc(BTabilities(res))
plot(res.qv, levelNames = playerset)

2番目より3番目の方が強いとは言い切るのはとても難しそう.

参考にしたものはスライドに載せた.

【学習予定】

Cattelan(2012) “Models for Paired Comparison Data: A Review with Emphasis on Dependent Data”に拡張のモデルが紹介されているので, 興味があるものを追っていけたらいいなと思っているが, 優先度はあまり高くないので, いつになるかは分からない. やる気が出たらそのうち.

Sports Analyst Meetup #3の後日譚々

Sports Analyst Meetup #3@ハロー貸会議室半蔵門(協賛:ダイナミックプラス株式会社)にてLTを行ってきました。

そのスライドがこちらになります。

speakerdeck.com

 

コードや分析の詳細記事は今週中に書きたいところです。

【追加 : 2019/07/02】書きました。

moratoriamuo.hatenablog.com

今回はBradley-Terryモデルの基本的なところを学習して、それを紹介するというものでした。本当はもう少し拡張したものまで学んで、ダブルスの試合結果からダブルスとしての戦闘力やペア間の相性などがデータから見出せたら面白いかなと思っていたのですが、そこまでは学習が間に合わずというところでした。最近、離散選択モデルの勉強をしていたおかげで、基本モデルの理解はスムーズに行きました。心理学のサーストンモデルや計量経済学の離散選択モデルを学んだことのある方はなじみがある数理構造じゃんってなったように思います。離散選択モデルのデータって自前では用意しにくいところもあるし、スタートとしてBradley-Terryモデルからこのような数理モデルを学んでいくのは、ありな気がしました。あと、レーティングとの繋がりがあるみたいなので、そのへんもこれから知れたら面白いかなというところです。

 

今回のspoanaはLT大会ということで18本ものLTがあって面白かったです。

spoana.connpass.com

以下でロングトークとLTに対する振り返りをメモ的に書いていきます。当日深夜に自分の思考絡めながらの走り書き小学生の感想的なメモなので、とても読みにくいとは思われます。発表内容と思考内容の区別がつきにくかったり、言葉の選択のゆれが酷かったりもします。ご了承ください。間違った理解をして書いている可能性もあるので、もしそういうところがあればご指摘いただけると嬉しいです。

 

・ダイナミックプライシングの話

自分が経済学部出身なので興味はある話題だった。ダイナミックプライシングの話は航空券やホテルの価格設定で使われているという話も聞いたことがあって気になっていた。やっていることは供給側が需要関数に合わせて利潤を最大化するように動学的に価格設定する戦略決定なのだろうと予測しているのだけれど、理論と実際の実装が気になる。ダイナミックプライスの需要側と供給側の戦略の均衡が気になる。需要側としてはどういう戦略が良い戦略なのか。属性にも依存する話だろうが。均衡に達するまでに、我々、需要側はどういう買い方が合理的なのか分からず右往左往する気がするのだけれど、そういうところも気になる。一需要者として、どういう対応すべきか、どのタイミングで買うべきか、安く買うことがいいのだろうけど、どのタイミングなのか分からないし、考えることによる効用の低下もなんて考えたりして。賢い人はみんなダイナミックな話やっとるって再認識。

・バドミントンの話

バドミントン!!動画から自力でデータ化しててすごかった。データでいかに攻めに持っていけるかっていうことを確認していたが、まさにその通り!経験から来るものと一致する!と頷きながら聞いていました。これからもバドミントンネタが増えるといいなと思います。

 

・フリスビーの話

Kinect使ってフォーム修正して初級者を中級者にする話。感覚で上達してきたところをデータ取って明確にして言語化、可視化してもらえると初級者はとてもありがたいよなぁと思うし、こういうのすごくいいと思います。

 

・オリンピックの自国開催でパワーアップの話

メダルでスコア化して正規化してカルマンフィルタで単純な予測。メダルが増える原因が気になった。応援で力が出るのか、環境に適合しているからか、育成にお金が使えるようになるからか、などと考えながら聞いてました。

 

MTGの話

カードゲームも不確実性の中でいかに情報をもぎ取って有利に確率ゲームを展開するかって話だなと共感。自分のLT直前でメモがとれず。

 

・レーティングの話

レーティングの話の入門ができてよかった。イロレーティングとBradley-Terryモデルはつながっているぽいので、そのへんが気になった。『レイティング・ランキングの数理』という本を買ったので、ちょっとずつ読み進めたいところ。レーティングっていろんな指標から一次元に乗せるってとこだと思うのですが、これって面白いけど難しい問題だなって。とくに自分が使わなきゃいけない話ではないので、ゆっくりになるとは思いますが、もう少し話を追っていけたらなというところです。

 

・自分のテニスのスコア分析

セットの間で、メモしてデータ化。40-0でリードしているときは試したいことやってみたり気が緩んだりでミスりやすいみたいなのが条件付きで平均とることで見えてきて、あぁいいなってなってました。そういうところは自分の中でも実際にあって、少しレーティングの話に戻りますが、得点差で評価するタイプだとどう対応するのだろうかって気になったりしました。感覚としては、個人競技であれば、そこでミスしても気を引き締めていけば勝敗としては問題なさそうだけど、団体競技だと連携が崩れて勝敗が転がるということが起こりそうっていうものがあって(たとえばサッカーの2-0があぶないみたいなのはよく聞きます。データとしてはどうなのと今思いました。)、そのへんをモデルで分析出来たらいいなってアイデア模索中です。懇親会で話してて、『ベイビーステップ』はいいぞってなりました。

 

ラグビー

すごいタックルに10ポインツ!!DEA(包絡分析法)って名前となんとなくやっていることは知っているけど、なんとなく以上を知らないのでもう少し知っておかなきゃだなと。病院とか図書館とかの効率性評価に使われますよね。

 

・野球の戦評の自動生成の話

Win Probability Added(WPA)。試合展開打撃状況投手敗因の項目で自動生成。仕事の諸々文書作成も、こんな感じで自動でやっていけないかなと思いながら聞いてました。自動生成した戦評を数パターン用意しておいて、1パターンを発表するっていうやり方のライブ感がいいなって思いました。他のパターンもどこかで見れたら面白いなって思いました。

 

・野球の盗塁に関する心理バイアス

評価手法と発表がとても面白かった。ビハインドではリスクとりたくなるっていうのはすごく分かる(競馬で負ける人並みの感想)。後でもう一度スライドを見返してみたいと思った。リスクに関しては前回、自分が発表したLTで扱ったのですが、自分の勝負における関心事として大きなもので、リスクコントロールについて、まだ言えそうなことがあるので、ぼんやりとある思考を形にできたらいいなと思っています。

 

・野球の満塁の話

スライドじっくり見ながら聞いてて全然メモがない。DeNA、結構満塁策とること多いよなと思っていたけど、やっぱり多いし、今日の広島戦の国吉に痺れました。野球データは豊富にあるので、そのへんをいじれるようになったら楽しいだろうなと思いつつ手を出せずにいる。"Analyzing Baseball Data with R, Second Edition "が気になるところ。

 

軟式野球のバッドの話

1点とれば勝てる世界の軟式野球への革命。びよんどまっくす。金属に比べてバッドの買い替え頻度が高まり、学生には金銭的に辛いという。ビジネスな何か。素材で競技が変化するってのはバドミントンのラケットでもあったはず。ネタになるかも。

 

・サッカーのトラッキングパスの話

走っている総距離はプロといい勝負。1本あたりのダッシュの距離が短くなっていくところは違うんだとか。サッカーが5分で数回程度しかボールを蹴らないっていうような一言があったと思うんですが、これはこういう競技における本質で、大事なのは触れていないときの動きなんだ!!みたいなのがたぶんあって、そのへんを上手く示せたり、そういう動きを定量的に評価できたらいいよなぁと思う次第です。

 

・サッカーの筑波大学の蹴球部データ班の取り組みの話

分析体制半端ないってと感心。負担が重すぎる分析は軽くしていくみたいな取り組みはたしかになってなりました。

 

・sportech sanityの紹介の話

スポーツデータ分析の知見はいろんな分野に散らばっているのでレポジトリを作ってまとめていこうというもの。統計学の勉強してても、経済学だったり理工系だったり心理学系だったり医学系だったりプログラミングのところにあったりと散らばっていて大変だなと昔から思っていて共感。

 

・サッカーのアナリストの話

気付かぬところでドラえもん送りバントしてた。勘と経験だけでなく、勘と経験とデータだ!って話と不正解を選ばせないっていう話は、うんうんとなりながら聞いてました。後者に関しては、相手が強いと悪手は厳しく咎められるということがある。対戦ものに関しては本当にそれを実感する。「ミスを咎める」みたいな言い方は将棋の言葉だと思うんだけど、広く知り渡るといいなと思っています。

 

・サッカーのPKの先攻有利の話

順番を「ABBA方式」にすることで先攻有利が変わるかどうかを確率シミュレーションで確認。こういう単純なところから確認していくスタイル、個人的に好きです。

 

・サッカーのサラ―効果の話

リヴァプールでのサラーというエジプト人選手の活躍で、イギリス人のイスラム教への感情が変わっているのかどうかを因果推論している論文の紹介。この論文は気になるので要チェックや。サラー以外でもそういうことありそうだし、他への適用例が気になる。あと、クリロナがレアルから抜けてどうこうみたいなのを前読んで、クリロナのピアエフェクトとか見れたら面白そうだなって。

 

登壇する機会を与えてくださったSports Analyst Meetup運営の皆さま、参加者や関係者の皆さま、ありがとうございました。人生初安打報告への温かい拍手もうれしかったです。今回も楽しく参加できたので、またネタを探して発表出来たらいいなと思っています。その際はよろしくお願いいたします。

瞬間、心、重ねて

 人生は瞬間瞬間における選択の連なりである。選択をするまでに与えられる猶予の時間も使用可能なエネルギーもけっして多くない。制限された時間の中で、持っているエネルギーの中で、選択を続けていく。そんなわけで、余計な時間やエネルギーを消費しないように、いつもやっている動作が無意識に出てしまうものだ。それらの行為は時間とエネルギーに対する制約のもとで最適化されているものなのだ。周りの人の真似を選択するのも、情報を集めて最善の手を計算する余裕がない場合にはある程度有効になるため起こりうると考えられる。制約に縛られながら人は選択を繰り返していく。
 普段からやっている動作は、過ごしてきた環境の中で見よう見まねで無意識に覚えてきたものであったり、どこかで意識して実行していく中で習慣化して無意識でもできるようになったものであったりする。そういった動作は、ぎこちなさがなく、スムーズに行われる。演技で一時的に真似することはできるかもしれないが、それなりのコストを代償にしない限りどこかできっとぼろが出る。演技の中に不自然さが垣間見える。自然かどうかは見てれば分かる。
 人はそんな自然な一挙手一投足に魅了されたり、幻滅させられたりする。ぽろっと出た、本当に些細な言葉で、テンションが上がったり、モチベーションを損なったりする。受け身目線の言葉を並べてしまったが、きっと私も誰かに対してそのような感情の変動を起こさせてしまっているのであろう。誰がどう思っているのかなんてことの正解不正解の判断なんて基本的に不可能なことなのだが。
 このことは、表現と印象の問題として一纏めにできる(表現と印象という言葉が対義語であることを知ったのは高校生の頃だったが、かなり衝撃的だったのを覚えている。ExpressionとImpression なわけだ。)。日常会話も芸術作品も何でもかんでもキャッチボールの投げ手と受け手は真逆の立ち位置にあって、互いにそれぞれの意図するところは分からない。表現者が何を考えて表現を発信したのかに関わらず印象者は好きに何かを考えて受信する。そんなものである。この文章もそのような中で書かれている。
 表現をする側としてはできるだけ正しく意図したものを伝えていきたいし、印象を受ける側としてはできるだけ正しく意図されたものを受け取っていきたい。表現するときは、先に書いたように制約がつきまとうとはいえ、エネルギーを割いて推敲して表現を選択していきたい。印象を受けるときは、表現の意図の拡大解釈や全く異なる解釈をしてしまうこともあるだろうし、それもまた感情を動かしていくものなのでが、第一としてはその限界を頭に置きながらも正しく理解することを心掛けていきたいものだ。
 「白は、完成度というものに対する人間の意識に影響を与え続けた。」という一文を思い出す。原研哉の『白』だ。そして、高校の美術の教師が、白いキャンバスに私が気分で重ねた油絵具を一目見て、何を重ねたのか言い当てたことも思い出した。頭の中はそんな感じでこの文章を書き終える。

{mlogit}パッケージ(混合ロジットモデル)の学習記録

  

【学習動機】

前回の続きで少し複雑なロジットモデルをRで推定するために, {mlogit}パッケージの使い方を覚えたい.

moratoriamuo.hatenablog.com

調べていたら, Train先生が{mlogit}の使い方のエクササイズを書いているのを発見した.
Kenneth Train’s exercises using the mlogit package for R
http://www.idg.pl/mirrors/CRAN/web/packages/mlogit/vignettes/Exercises.pdf
これに沿って学習していく. 今回は, 混合ロジットモデル(ミックスドロジットモデル(mixed logit model))の節を扱う.
ドキュメントの全てを扱うわけではない. また変数名やコードなど自分の理解しやすい形に変えているところもあるので注意されたい.

Exercise 3: Mixed logit model

ググって探索していたら, 二週間前にこれがアップロードされていたことに気付いた. 英語が読めるならこっちを読んだ方が絶対にいい気がする.

加えてこれも発見した. 理論面はTrainのテキストに沿っていて, 補完的に読めてとてもよい.
{mlogit}パッケージ作った人の解説
Yves Croissant “Estimation of multinomial logit models in R : The mlogit Packages”
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.303.7401&rep=rep1&type=pdf

【学習内容】

3. 混合ロジットモデル

ライブラリ読み込み.

library(mlogit)
library(tidyverse)

データの確認

選択実験データ(離散選択実験データ, 選好意識データ(Stated Preference Data), コンジョイント分析のデータ)と呼ばれるもの(だと思う). それぞれ異なる属性を持つ仮想の電力会社を4つ設定し, 実験参加者の個人がどの電力会社を選ぶかというデータ.
361人が実験に参加. 各個人は12個の想定パターンに回答. 12個全てに回答していない参加者もいる. そのため, 総回答数は4308となっている (全員が12個全てに回答していれば, 361*12=4332となる)

data("Electricity", package = "mlogit")
head(Electricity)
##   choice id pf1 pf2 pf3 pf4 cl1 cl2 cl3 cl4 loc1 loc2 loc3 loc4 wk1 wk2
## 1      4  1   7   9   0   0   5   1   0   5    0    1    0    0   1   0
## 2      3  1   7   9   0   0   0   5   1   5    0    0    1    0   1   1
## 3      4  1   9   7   0   0   5   1   0   0    0    0    0    1   0   1
## 4      4  1   0   9   7   0   1   1   0   5    0    0    1    0   1   0
## 5      1  1   0   9   0   7   0   1   0   5    1    0    0    0   0   1
## 6      4  1   0   9   0   7   0   0   1   5    0    0    1    0   0   0
##   wk3 wk4 tod1 tod2 tod3 tod4 seas1 seas2 seas3 seas4
## 1   0   1    0    0    0    1     0     0     1     0
## 2   0   0    0    0    1    0     0     0     0     1
## 3   1   0    0    0    1    0     0     0     0     1
## 4   0   1    0    0    0    1     1     0     0     0
## 5   0   1    1    0    0    0     0     0     1     0
## 6   0   1    0    0    1    0     1     0     0     0

変数
choice:個人の選択(1,2,3,4) id:個人のid
pfi:price fixed. 選択肢iの固定価格[cent/kWh] 価格は供給元や実験によって異なる、kWhあたりの表示セントでの固定価格. 価格はサプライヤーや実験によって異なる.
cli:contract length. 契約期間[年]. 他のサプライヤーに切り替えた場合, 違約金を支払うことになる. いつでも契約を解除できる契約も結ぶことができ, その場合は0と記録される. (携帯の2年縛りと同じようなものか.)
loci:local. サプライヤーが地元企業かどうか.
wki:well-known. サプライヤーが有名な企業家どうか.
todi:time of date. 午前8時から午後8時まで11cent/kWh, 午後8時から午前8時まで5cent/kWhというレートになるプラン.
seasi:夏季は10cent/kWh, 冬季は8cent/kWh, 春と秋は6cent/kWhというレートになるプラン. 価格は電気のみで, 地域の公益事業者による送配電は含まれていない.

wide形式から変換

Electr <- mlogit.data(Electricity, id = "id", choice = "choice", 
                      varying = 3:26, shape = "wide", sep = "")
head(Electr)
##     choice id alt pf cl loc wk tod seas chid
## 1.1  FALSE  1   1  7  5   0  1   0    0    1
## 1.2  FALSE  1   2  9  1   1  0   0    0    1
## 1.3  FALSE  1   3  0  0   0  0   0    1    1
## 1.4   TRUE  1   4  0  5   0  1   1    0    1
## 2.1  FALSE  1   1  7  0   0  1   0    0    2
## 2.2  FALSE  1   2  9  5   0  1   0    0    2

モデルのあてはめ

mlogit()で混合ロジットモデルするときの引数:

rpar : 変数名に分布を指定することで, その変数のパラメータに分布を当てはめる. “n”ならば正規分布, “l”ならば対数正規分布, “t”ならば切断正規分布, “u”ならば一様分布.

R : heterosc = TRUEの場合はガウス求積法に使用する関数の数. rparがNULLの場合は擬似乱数の数. (heteroscのデフォルトはNULL, Rのデフォルトは40)

print.level : printメッセージの詳細の指定. 0,1,2のいずれかをとる. 0ならば,最適化プロセスに関する情報は記載しない. 1ならば, 尤度の値,ステップおよび停止基準を記載する. 2ならば, パラメータと勾配ベクトルも記載する.

panel : rparがNULLではなく, かつデータが同じユニットの観測の繰り返しの場合に関係する. TRUEの場合、混合ロジットモデルはパネル手法を使用して推定される. (panelのデフォルトはNULL)

halton : rparがNULLでない場合のみ関係する引数. 擬似乱数の代わりにHalton sequence(ハルトン数列)が使用される. Halton数列は, 素数を選んで,(0,1)区間を分割することを繰り返して 分布範囲を均等に網羅し,かつ相互に相関などを持たない数列を準乱数として用いる方法.
使用する素数と削除する要素の数のリストを引数にとる. NAを指定したら, デフォルト値が使われれる.

Halton数列については東京大学工学部都市工学科都市生活学/ネットワーク行動学研究室の混合ロジットモデルに関する説明資料に依っている. http://bin.t.u-tokyo.ac.jp/kaken/pdf/mxl.pdf
詳しくはTrainのテキストのChap.9に説明がある. (https://eml.berkeley.edu/books/choice2.html) ここの研究室は東大の社基の羽藤研あたりが集まってるところで, 離散選択モデルに関する資料がたくさんアップロードされている. 土木・交通分野の個人の行動予測, 交通需要推計の研究が盛んみたい.
この研究分野の研究事例としては, 神奈川県のシーサイドラインが開通される前に需要予測を行っていたものなどがある.
http://library.jsce.or.jp/jsce/open/00039/1988/11-0455.pdf

using 100 draws, halton sequences and taking into account the panel data structure

混合モデル1

切片なしで, モデルの係数パラメータ6つに対して正規分布を当てはめるモデル. Halton数列を使用し, パネルデータ構造を考慮するケース.

推定には多少時間がかかる.

Elec.mxl <- mlogit(choice ~ pf + cl + loc + wk + tod + seas | 0, data = Electr,
                   rpar = c(pf="n", cl="n",loc="n",wk="n",tod="n",seas="n"),
                   R = 100, halton = NA, print.level = 0, panel = TRUE)

summary(Elec.mxl)
## 
## Call:
## mlogit(formula = choice ~ pf + cl + loc + wk + tod + seas | 0, 
##     data = Electr, rpar = c(pf = "n", cl = "n", loc = "n", wk = "n", 
##         tod = "n", seas = "n"), R = 100, halton = NA, panel = TRUE, 
##     print.level = 0)
## 
## Frequencies of alternatives:
##       1       2       3       4 
## 0.22702 0.26393 0.23816 0.27089 
## 
## bfgs method
## 21 iterations, 0h:0m:53s 
## g'(-H)^-1g = 3.3E-07 
## gradient close to zero 
## 
## Coefficients :
##          Estimate Std. Error z-value  Pr(>|z|)    
## pf      -0.973384   0.034324 -28.359 < 2.2e-16 ***
## cl      -0.205557   0.013323 -15.428 < 2.2e-16 ***
## loc      2.075733   0.080430  25.808 < 2.2e-16 ***
## wk       1.475650   0.065168  22.644 < 2.2e-16 ***
## tod     -9.052542   0.287219 -31.518 < 2.2e-16 ***
## seas    -9.103772   0.289043 -31.496 < 2.2e-16 ***
## sd.pf    0.219945   0.010840  20.291 < 2.2e-16 ***
## sd.cl    0.378304   0.018489  20.461 < 2.2e-16 ***
## sd.loc   1.482980   0.081305  18.240 < 2.2e-16 ***
## sd.wk    1.000061   0.074182  13.481 < 2.2e-16 ***
## sd.tod   2.289489   0.110731  20.676 < 2.2e-16 ***
## sd.seas  1.180883   0.109007  10.833 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -3952.5
## 
## random coefficients
##      Min.     1st Qu.     Median       Mean     3rd Qu. Max.
## pf   -Inf  -1.1217350 -0.9733844 -0.9733844 -0.82503376  Inf
## cl   -Inf  -0.4607190 -0.2055565 -0.2055565  0.04960589  Inf
## loc  -Inf   1.0754783  2.0757333  2.0757333  3.07598832  Inf
## wk   -Inf   0.8011189  1.4756497  1.4756497  2.15018054  Inf
## tod  -Inf -10.5967791 -9.0525423 -9.0525423 -7.50830550  Inf
## seas -Inf  -9.9002649 -9.1037717 -9.1037717 -8.30727842  Inf

平均的な係数を持っている客が契約期間を1年伸ばす(減らす)のに追加的に欲しい(支払う)額を調べる.

coef(Elec.mxl)["cl"]/coef(Elec.mxl)["pf"]
##        cl 
## 0.2111772

契約期間の係数は平均-0.2, 標準偏差0.37の正規分布に従っていると推定されており, その場合, 係数が負となっている割合を計算すると

pnorm(-coef(Elec.mxl)["cl"]/coef(Elec.mxl)["sd.cl"])
##      cl 
## 0.70656

となり, 7割が契約期間が長くなることを嫌うと解釈できる.

同様に固定価格の係数についても見る.

pnorm(-coef(Elec.mxl)["pf"]/coef(Elec.mxl)["sd.pf"])
##        pf 
## 0.9999952

ほぼ100%, 係数が負となって整合的.

混合モデル3

よく知られた企業の方が好む人もいるだろうが, 好まない人も一定数いると考えがある場合, well-knownの係数の分布は正規分布ではなく一様分布であるというモデルを当てはめる

Elec.mxl3 <- update(Elec.mxl, rpar = c(cl="n", loc="n", wk="u", tod="n", seas="n"))
summary(Elec.mxl3)
## 
## Call:
## mlogit(formula = choice ~ pf + cl + loc + wk + tod + seas | 0, 
##     data = Electr, rpar = c(cl = "n", loc = "n", wk = "u", tod = "n", 
##         seas = "n"), R = 100, halton = NA, panel = TRUE, print.level = 0)
## 
## Frequencies of alternatives:
##       1       2       3       4 
## 0.22702 0.26393 0.23816 0.27089 
## 
## bfgs method
## 22 iterations, 0h:0m:55s 
## g'(-H)^-1g = 3.71E-07 
## gradient close to zero 
## 
## Coefficients :
##          Estimate Std. Error z-value  Pr(>|z|)    
## pf      -0.882232   0.032818 -26.883 < 2.2e-16 ***
## cl      -0.217127   0.013776 -15.761 < 2.2e-16 ***
## loc      2.099319   0.081056  25.900 < 2.2e-16 ***
## wk       1.509410   0.065496  23.046 < 2.2e-16 ***
## tod     -8.607002   0.282983 -30.415 < 2.2e-16 ***
## seas    -8.602408   0.280671 -30.649 < 2.2e-16 ***
## sd.cl    0.381070   0.018150  20.996 < 2.2e-16 ***
## sd.loc   1.593850   0.087802  18.153 < 2.2e-16 ***
## sd.wk    1.786377   0.125764  14.204 < 2.2e-16 ***
## sd.tod   2.719078   0.119357  22.781 < 2.2e-16 ***
## sd.seas  1.945377   0.103686  18.762 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -3956.7
## 
## random coefficients
##            Min.     1st Qu.     Median       Mean     3rd Qu.     Max.
## cl         -Inf  -0.4741552 -0.2171273 -0.2171273  0.03990058      Inf
## loc        -Inf   1.0242834  2.0993191  2.0993191  3.17435473      Inf
## wk   -0.2769665   0.6162218  1.5094101  1.5094101  2.40259840 3.295787
## tod        -Inf -10.4409924 -8.6070022 -8.6070022 -6.77301190      Inf
## seas       -Inf  -9.9145449 -8.6024084 -8.6024084 -7.29027193      Inf

well-knownの係数の分布情報を確認する.

rpar(Elec.mxl3,"wk")
## uniform distribution with parameters 1.509 (center) and 1.786 (span)

plotで分布の形状も確認できる. 係数が負になるところは色が塗られて表題にその割合が書かれるようだ. 先のサマリーで確認できるが, -0.27から3.29の一様分布.

plot(rpar(Elec.mxl3,"wk"))

【学習予定】

プロビットモデルもある(https://cran.r-project.org/web/packages/mlogit/vignettes/e4mprobit.html)が, ロジットモデルに関してとりあえず{mlogit}の使い方が分かったので, いったん離散選択モデルの扱いの学習はここで止める.   プロビットの学習やネスティッドロジットや混合ロジットのシミュレーションデータを作成と推定をしてみたりしたいところだが, 他に時間と労力を割かなければなのでまたいつか.
今回やっていて数値計算や最適化周りの知識の足りなさを感じたのでそのへんの補強も必要そう. このへんはTrainのテキストの後半パートが役立ちそう.  

経済学の産業組織の分野のBLPアプローチは混合ロジットモデルを, 内生性を考慮して何かしているものっぽい. 産業組織は分からないので正しいか分からないがここにメモしておく.