チラシの裏の落書き日記

統計とか,研究とか,日常についての備忘録的なもの。

統計学に主義はあってもいいけど,人を説得するときに使うものじゃないという気がする

以下の記事を読んで個人的に思ったことをまとめておく。
hidekatsu-izuno.hatenablog.com

 個人的な思いは,タイトルの通りだけどそれ以外にも色々思ったことがある。この文章はなにかの戦いに決着を付けるためのものでもないし,現状の交通整理をするためのものではない。また,私は渡辺研究室とはなんら関係ないし,この文章は渡辺先生を含め,関係者の総意を表すものではなく,独学で統計の理論を学んでいる個人の限定的な解釈に基づくものとして理解してください。内容についての責任はあくまでも書いた私にあります。ガーッと書いたので不正確な点もあるかもしれないので,その点は先にお詫びいたします(怖い議論は怖いので,優しくしてください)。

 「渡辺ベイズ解釈」というものは,便宜的に名前をつけたのかもしれないが,適切には思えない。渡辺先生が作ったのは,特異学習理論であって,特定のベイズ統計的な方法ではない。渡辺先生の定義によると,ベイズ推測は,統計モデルと事前分布の組から予測分布を作って,それが「真の分布(データが従う分布)」に近いだろうと考える方法である。そして,予測分布の真の分布による平均を汎化誤差を考えたり,自由エネルギーを考えて,としてその確率的な挙動を解明する方法論として,特異学習理論があるのだと思う。そして,特異学習理論のすごいところは,真の分布がわからかなくても,特異なモデルであっても,(かなり少ない仮定のもと)手元にあるデータのみから,汎化誤差や自由エネルギーの不偏推定量を計算できる,という点にある。
 そして,渡辺先生は,ベイズ統計学の構成要素である,真の分布・確率モデル・事前分布の3つについて,それらに依存しない数学的な構造を明らかにしたのであって,なにか新しい結果の解釈を見方を提供したわけではないと思う。真の分布を想定することの是非は議論があるかもしれないが,ポイントはそれが存在するというよりは,あると想定すると予測や汎化などの重要な概念が定義できる,ということだと思われる。ここで,確率モデル,事前分布は人間が設定するものであるし,正しい推測方法は存在しないということは繰り返し強調されている。また,真の分布を想定しているから特異学習理論の範囲が限定されているという意見もあるかもしれないが,この理論自体かなり広範囲で使える設定になっており,すごく強い制約には思われない。どの方法も多かれ少なかれ制約はあるだろうが,個人的には,特異学習理論の成果は現実でも広く使えるようなものに思われる。どうしてもこの仮定がいやなら,汎化誤差などについて,数学的な保証はない別の理論を使っても,もちろんよい。私は数理的に保証された安心できる方法を使いたいと思っているが,選択は分析者の自由だ。
 やや脱線してしまったが,確率モデルと事前分布が想定されていることが,ベイズ統計学の最小限の要求であると思われるし,それ以上の「解釈」の部分は,少なくとも渡辺先生の理論の中には存在しないと考えられる。また,データはなにかの分布から発生していると考えるのであれば,それは確率変数とみなすことに,特に問題があるわけではない。データを確率変数と考えなければ,汎化誤差のデータの出方についての平均なども考えられられないと思う。これは統計量の理論的な性質を明らかにするための設定としてのみならず,実際に手元にあるデータ(確率変数の実現値)を使って計算したものは,ゆらぎがあるものの一つの現れ方に過ぎない,という点を思い出させてくれる意味で意味があると思う。 そして,ベイズ統計なのにデータを定数として扱っていないとしても,汎化誤差の推定値を計算する場合には,データは確率変数の実現値として(いわゆる「定数」ではないが固定されたものとして),扱うわけで,実際のデータ解析上大きな違いがあるわけではない。
 また,予測の観点やデータの出自の構造を考えるというのは,統計に限らず,科学が目指すもののひとつなのではないかと思う。その意味で,渡辺先生のベイズの捉え方を,(一貫したものの見方としての)主義,としたいと考える事もできるのかもしれないが,普通の科学がどこでも目指していることなので,「主義」といって区別をつけることなのだろうか。さらに,他の推測法,たとえば最尤推測であっても,(モデルの)パラメタの推定値をプラグインして,予測を行うし(そのためにパラメタを推定する),やはり渡辺先生の独自のものではないと思う。
 以上のことから,渡辺先生の特有のベイズ統計学の解釈があるわけではなく,ベイズ統計学のうち数学的に定式化できる最小限度のパーツ(統計モデル・事前分布)を使って,予測を行ったときの数学法則を解明したのが重要な点であると思う。これは,ベイズ統計学の中の成果ではなく,数学を使ったベイズ統計学の評価についての成果といってよいと思われるため,やはり「渡辺ベイズ解釈」という言い方には疑問がある。

 また,主観的・客観的,という用語について,「主観確率ベイズ統計学を特徴づけている」というのが,よくあるベイズ主義と言われるようなものになっているが,必ずしもそうではないと思われる。先に述べたように,ベイズ統計学の論理構造の最小限の部分だけを見たときにには,主観確率を含む必要がないからである。私は,主観確率が間違っているとか使うべきでないと言いたいわけではない。
確率をどう解釈するかは分野や個人の自由だと思う。そこに「正しさ」はないと思う。しかし,「ベイズ統計学を採用したから確率の解釈を主観的にしてよい」,というのは違うと思う。こうした考え方は標準的ではないかもしれないが,統計的な数式上の推論の論理構造と,現実での解釈や意味を切り分けることは,むしろ,様々な統計的な方法の共通点や違いを理解する上で意味があると思う。また,確率の解釈を使い分けることによって,新しい発見ができるのであればそうすればいいし,とくに一つになにかの確率の解釈に固執しなくてもよいと私は思っている。しかし,こうしたが統計的手法を特徴づけていると考えることに,意義はあまり感じない。
 統計学が扱う範囲として,こうした確率の解釈も含むべき,と考えるのも一つの立場だと思う。しかし,統計的な方法が様々な分野で普及しており,各分野の関心に合わせて手法が開発されている状況で,統計がカバーすべき範囲というのは,様々にあって然るべきなのだろうと思う(機械学習ではデータは勝手に溜まっていくようなものだが,社会調査などではデータの取得方法から設計する必要がある,など)。その意味で,統計学の数学的な構造の部分を考えるというのは,どの分野でも共通してできる最小限の部分のように感じられる。数学的な成果としていえることはそれはそれとして,実際の分野で方法の前提や結果の具体的な解釈を考えていけばよいのではないかと思う。これは,むしろ数学でできることとできないことをきちんと区別した慎重な態度であろう。科学的なものとはなにかは私にはわからないが,現実とは関係のない数式としての統計のなかだけでなにかできるのではなく,現実のなかでのデータの取得方法や実験の方法などが関わってくるものだと思う。その意味で,統計結果の解釈は現実のデータをどうするのかという点において生じるものでしかないと思われる。モデルはモデルで,想像上の産物でしかない。しかし,そのモデルが現象を予測するとかよく機能したり,厳密なテストにパスして有用そうだと認識されるなど,実際のデータで確認されたときにはじめて意味や解釈がとなってくるように思われる。
 
 また,もとの記事では,「統計学とは哲学的概念を数学を使ってモデル化したものと捉えた方がよい」と述べて,「渡辺ベイズの主張は、~中略~(極端に言えば)ツールに過ぎない数学だけで何かが語れると思うのは相当に極端な態度に思える」といっているが,先に述べたように,特異学習理論により,ベイズ統計学の数学的な特徴を明らかにしたことで,現実によくある特異モデルでも安心して使える基盤を提供しているわけで,数学だけでなにかの現実を語っているわけではないと思う。むしろ,モデル(あるいは仮説)を現実のデータでテストするための非常に一般的な方法を提供している(余談だが,検定も情報量基準もモデル選択という意味で同じような枠組みとして考えられる)。予測や汎化というものはもしかしたら一般的な考え方ではないかもしれないし,「哲学的な概念」ではないかもしれないが,私には,渡辺先生の成果は極めてベイズ統計学がなぜ現実で機能するのかの保証を与えているように見えるし,かなり広い現実的な状況で一般に使える統計的なツールだと思う。統計学の一部の考え方が,「哲学的な概念の数学的なモデル」といえることもあるのかもしれないが,少なくとも渡辺先生の成果は,そういったたぐいのものには思われない。

 さらに,「「合理的な(逆向きの)推定 or 意思決定」という概念があり、その主要な実装としてベイズ統計学が存在しているのだ」と言っているが,ベイズ統計学の一部にこういう考え方があるのは否定しない。しかし,意思決定以外にも事後分布でモデルを平均することで予測分布を作って予測を行うという営みもあり,こうしたことは哲学的な概念が先に来ているというわけではないと思われる(予測が大事,という点を哲学的と呼ぶなら哲学的なものが先に来ているのかもしれないが,わざわざそういう必要はないように思われる)。さらに,ベイズ統計学が機能する理由は「主観にもとづいている」かどうかに関わらず,数学の定理として証明されている事柄である。意思決定と予測は目的が違うので,それぞれに合った確率の解釈をすればよいが,いいたいこととしては,「ベイズ統計学だから確率は主観的だ」ということは必ずしもそうではないということだ(渡辺先生成果がマイナーであり,無視できるものとみなしたい人もいるかもしれないが,WAICの成果などは順調に引用数が増えていて,かなり認知されていると考えて良いのではないかと思う)。もちろん,学習機(統計モデル=ただの同時分布)が,どう判断するのかという点を強調して,主観確率的な解釈をしたければそうすればいいかもしれないが,最終的な結論は対して変わらないし,本質的なものには思われない。(脱線:こう考えると,最尤推定したパラメタを使った学習機の,”主観的な確率”というものを考えられそうだし,その意味でベイズじゃなくても主観的な確率になるのかな。)
 また,もとの記事の最後に,「頻度論に基づく統計学は「合理的な(逆向きの)推定 or 意思決定」のために作られていない」と述べられているが,これにも同意できない。頻度論的な,統計的仮説検定は,特定のデータが与えられたときに仮説の採択・可否を判断しているし,信頼区間に推定の正確さも評価しているし,まさに推論・意思決定のための道具であるように見える。また,先に述べたように,検定も仮説的なモデルをアセスメントするためのツールとして考えることもできる。この意味でベイズ統計学だろうと頻度論的な統計学だろうと同じことをやっている(もちろん,これは一部の活動が重なっているという事実を指摘しているだけで,両方の統計学の目的がすべてが同じというわけでない)。こうしたこともあり,統計の主義を分ける必要はやはり,特にないと考えられる。もちろん,現場の状況によって何を使えるのか,使えないのかを選択する必要があるが,それは統計的な主義が決めることではないと考えられる。

 とくにまとまった文章ではないが,「渡辺ベイズ」として言及されていることについて,ディフェンドする形になった。いずれにせよ言いたいこととしては,
1.統計が扱う範囲はかならずしもコンセンサスがあるわけではないので,モデルや数学的な部分と,主義の部分は分けて考えて,とくに主義的なものは現実場面に即して考えたほうがよいのでは,
2.主義はあるかもしれないが,それは例えばベイズ統計を選択したから一緒にくっついてくるものというわけではない,
3.主義があるから,統計的なものの良さが決まるわけではない,
4.「渡辺ベイズ」という用語は適当には思われない,
5.個人の自由でどの主義に依拠するかは好きにすればよいが,人にまでそれを強要しないでほしい,
ということだ。統計を使う人は,使いたい理論や主義を使って色々と解釈を頑張ればいいんじゃないかな。

Slice sampling(スライスサンプリング)の簡単な例をRでやってみた。

テラモナギさんの記事と,PRMLの11.4節を参考にしてみた。
https://teramonagi.hatenablog.com/entry/20130203/1359888495#fn-b0191759
本当はもとの論文を見るべきだが,大体あっていそうなのでOKなような気がする。


(1変数の確率密度からの)スライスサンプリングのアルゴリズムまとめ

  1. サンプリングのための変数(x)の初期化(x_oldとする。)
  2. 目標とする確率密度(規格化していないものでOK)をfとしたとき,一様分布(min=0, max=f(x_old) )からスライスする位置uをサンプリングする。
  3. {x: u < f(x)}となる,領域から一様にx_newをサンプリングして,x_newをサンプルとして保存し,x_old = x_newとする。
  4. 2~3を繰り返す。

3を実行するときに,lower_xとupper_xを適当に決める。このときに,f(lower_x), f(upper_x) < uになるように,lower_x, upper_xを引き伸ばす必要がある(steppin out)。
さらに,その後,一様分布(min=lower_x, max=upper_x )からx_newをサンプリングし,f(x_new) < uとなってしまっていたら,lower_x, upper_xをx_newに置き換えて(shrinkage),再度x_newを一様分布からサンプリングする。
f(x_new) >= uであれば,x_newを新しいMCMCサンプルとして保存する。
x_newはスライスの範囲,つまり,{x: u < f(x)}に入っていなければならないので,この条件に入らないx_newはサンプルとして採用できない。
この辺のややめんどい手続きが必要になるが,これは,{x: u = f(x)}となるxの点を数値的に見つけないといけないから。
しかし,メトロポリスヘイスティングアルゴリズムのような,採択・棄却サンプリングのように,提案分布の調整がいらないので,とても楽。
しかも,確率密度関数カーネルだけわかっていればよいので,汎用性が高い。
多変量の場合も,サンプリングしたいパラメタ以外のパラメタなどをすべて条件づけたものを利用すればよい。
(確率密度の部分については,対数をとったものを考えるというのも手。数値がアンダーフローになるのを防ぐため。)

#
# Slice サンプリング
#

#
# サンプリングを目指したい混合分布
# 
f <- function(x, 
              mixing=c(0.7, 0.3),
              mu=c(-1, 4),
              sigma=c(0.5,2)){
  y <- 0
  for(i in seq_along(mixing)){
    y <- y + mixing[i]*exp(-(x-mu[i])^2/(2*(sigma[i]^2))) 
  }
  y 
}

#
# 目標の分布
# 
curve(f, from=-3, to = 10)

# 
set.seed(111)

is_f_x_in_s <- function(u, f_x) u <= f_x

slice_one_step <- function(x_old, width = 1.0){
  
  
  #x_old <- rnorm(1)
  
  f_x_old <- f(x_old)
  u <- runif(1, min = 0, max=f_x_old)
  #f_u <- f(u)
  
  
  lower_x <- x_old - width*runif(1)
  upper_x <- x_old + width*runif(1)
  
  f_lower_x <- f(lower_x)
  f_upper_x <- f(upper_x)
  
  
  #
  # Stepping out
  # 
  while(is_f_x_in_s(u, f_lower_x)){
    lower_x <- lower_x - width 
    f_lower_x <- f(lower_x)
  }
  while(is_f_x_in_s(u, f_upper_x)){
    upper_x <- upper_x + width 
    f_upper_x <- f(upper_x)
  }
  
  #
  # Shrinkage
  # 
  repeat{
    x_new <- runif(1, min=lower_x, max=upper_x)
    f_x_new <- f(x_new)
    if(is_f_x_in_s(u, f_x_new)){break}
    
    if(x_new < x_old){
      lower_x <- x_new
    } else {
      upper_x <- x_new
    }
  }
  
  x_new
}


n_iter = 1000
slice_sampling <- function(n_iter = 1000){
  x_samples <- rep(0, n_iter)
  x_samples[1] <- x_old <- rnorm(1)
  
  for(i in 2:n_iter){
    cat("i = ",i,"\n")
    x_new <- slice_one_step(x_old, width = 0.5)
    x_samples[i] <- x_new
    x_old <- x_new
  }
  x_samples
}



res_x <- slice_sampling(n_iter = 10000)

hist(res_x,freq = F,ylim=c(0,0.4),
     breaks = seq(min(res_x)-0.5, max(res_x)+0.5,by=0.1),
     col=rgb(0,0,1,alpha=0.7),
     ylab="")
f_const <- integrate(f, lower=-10, upper=10)

#
# 数値的に規格化した密度
# 
f2 <- function(x)f(x)/f_const$value
curve(f2, from = min(res_x)-0.5 ,
      to= max(res_x)+0.5,
      col=rgb(1,0,0,alpha=0.8),lwd=2, 
      xlim=c(min(res_x)-0.5, max(res_x)+0.5),
      ylim=c(0,0.4),
      ylab="",
      add=T)
f:id:kazzstat:20200727070222p:plain
サンプリングの結果とサンプリングをしたい密度

Rでforeachをつかって楽にシミュレーションを分割してみる雛形

この雛形を使うと並列計算ができる。個人的には返り値などは基本的に全部書きだしてあとで分析する形にしておくと安心。
シミュレーションを関数にしてしまうことと,foreach関数の指定のときに必要なパッケージや環境の変数(自作の関数とか)を全部持っていてもらうように指定するのがポイント。

#
# シミュレーションの雛形
# 
rm(list=ls())
library("MASS")
sim_times = 1:100

simulation_fun <- function(sim_times){
  setwd("hogedir")
  for(i in c(500, 5000)){
    for(pop_cor in c(0, 0.5, 0.8)){
      #結果保存用フォルダ  
      dir.create(paste0("sim_",i,": pop_cor = ",pop_cor)) 
      
      for(n in sim_times){
        set.seed(i*10 + 100*pop_cor)
        cat(format(Sys.time(), "%Y %b%a %d %X"),  ": i = ",i,": pop_cor = ",pop_cor,": n = ",n,"\n")

        Sig <- matrix(pop_cor, 10, 10 )
        diag(Sig) <- 1
        mu = rep(0,K)
        
        x = MASS::mvrnorm(I,mu = mu, Sigma = Sig)
        
        time1 = as.POSIXlt(Sys.time())
        res <- list(mean=colMeans(x),cov = cov(x))
        time2 = as.POSIXlt(Sys.time())
        use.time = difftime(time2, time1, units="secs")

        #use.time
        save(list = c("res", "x","time1","time2","use.time"),file = paste0("sim_",i,"_",pop_cor,"/sim_",i,"_",pop_cor,"_",n,".Rdata")) 
      } # sim times 
    } # pop cor
  } # sample size
  return()
}


library("foreach")
library(doParallel)

package_list <- c("MASS")

#
# 指定したいコアの数
# 
n_div <- 4
# 
# いくつからいくつのシミュレーションを分割するか指定
# 
start_ind <- seq(from=1, to = max(sim_times),by = max(sim_times)/n_div)
end_ind <- seq(from=max(sim_times), to = 1,by = -max(sim_times)/n_div)
end_ind <- end_ind[order(end_ind )]
divide_indicator <- mapply(function(x,y)x:y, start_ind, end_ind,SIMPLIFY = F)
divide_indicator

cores = n_div 
cluster = makeCluster(cores)
registerDoParallel(cluster)

#
# 環境を全部持っていってもらう
# ls(envir=parent.frame())
# 
foreach (sim_times = divide_indicator, .combine=c,.packages = package_list,.export = ls(envir=parent.frame())) %dopar% {
  simulation_fun(sim_times)
}
stopCluster(cluster)

>||

とある論文その2を読んで考えたこと

また読んでいても気になる論文があったので,備忘録的にツッコミをメモしておく。
http://team1mile.com/sjpr61-1/shimizu.pdf


第2章のベイズ統計学の流れの整理区分は後の議論を明確にする上でもわかりやすいように思う。

さて,第2章2節から本題。

統計モデリングにおいてベイズ推定が注目されるようになったのは,Gelfand and Smith(1990)によってベイズ推定のためのアルゴリズムとしてMCMC が利用できることが示されたことが一因である。

とあるが,まず,ベイズ推定の定義が明示されていない。
おそらくこの文章でのベイズ推定の定義は,”パラメタの事後分布を求めること”ということだろう。
しかし,例によって,渡辺先生のHP(http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/Bayestheory.html
)から引用を行うと,

ベイズ推測とは,
「予測分布 p(x|D) は,きっと,真の密度関数 q(x) に近いであろう」
と推測する方法です.

とのことで,事後分布を求めることそれ自体は,ベイズ推定ではないと思います。
ただ,”(パラメタ)推定”と”推測”ということばの使い分けが微妙なので,難しい。

これまで統計モデルは最尤法による推定が一般的であったが,変量効果を含む一般化線~中略~広くベイズ推定の有用性が理解されるようになった。

とあるが,何を解くのか?最適化問題を解くという意味?
とても気になるのが,この論文は一貫して”MCMCベイズ推定”という認識で議論がなされているようだが,その理解は個人的には正しくないと思う。MCMCはあくまでも,事後分布の確率的な近似方法,より広くは期待値計算の汎用的な方法であり”アルゴリズム”であって”推定方法”ではない。他の力技の積分計算を使っても事後分布を実現することはできる。つまり,最尤推定値を得るために,(対数)尤度関数の最大値を求めるために最急降下法という”最適化計算のアルゴリズム”を使うことがあるが,最急降下法のことは”最尤推定法”とは呼ばないのと同じである。モデルと推測方法と計算アルゴリズムは区別しなければならない。結局,この論文で言うところの”ベイズ統計学”というのは,”パラメタ推定を行う方法の一つ”として認識されている(次の第3章に明確に,”ベイズ統計学はパラメータ推定の枠組みのひとつであって,統計モデリングの本質ではない。”といっているので,私の論文の読み取りは正しそうだ。)渡辺先生のいうところのベイズ統計学とはだいぶ違いそうな雰囲気がある。

心理学では多くの場合~中略~しかし,確率モデルとして現象を記述する場合,データが仮定した確率分布から生成されるという想定さえ成り立てばよい。~中略~確率モデルの表現で重要なのはデータを得る具体的な手続きではなく,データがどういう確率分布に従って得られるか,という点だけである。

とあるが,いまいちよくわからない。母集団を仮定しているとかなんとかとコイン投げの話をしているが,コイン投げだって母集団を仮定していると思う。”データが仮定した確率分布から生成されるという想定さえ成り立てばよい”といっているが,これが成り立つかどうかがわかっていれば統計学なんていらないわけで,データ生成モデルはあくまで仮想的なものであるという認識がなさそうである。また,マルチレベルモデルではデータがどのように取られたかを考慮したモデリングがなされているし,データ取得の手続きを無視したらiidの仮定なども無視することになりそうだがそのあたりはどう考えているのか。

データ生成メカニズムを確率モデルで記述するというのは,データがどのように発生するかを,確率分布を使って数理的に表現する,という意味であることがわかる。

結局,確率モデルとはなんなのか明確にはわからない。

その点においては 2 つの平均値の差についての t 検定(以下,単に t 検定と呼ぶ)も統計モデリングとして~略

検定は検定であって,統計モデルではない。

3.2 心理学における伝統的なデータ分析手法との違い

データ分析手法と統計モデリングは直接比較できるものではないと思う。というか何を比較しているのかわからない。

次に参加者がどのような行動を行うかについての予測を導く目的で使われることが少なく,またその関心からも外れている。それは,伝統的な心理学の方法論がおもに行動の説明に重点が置かれており,予測が重視されていないからであると考えられる。

しかし,例えば実験計画法に則ってデータを収集し平均的な行動の傾向の違いを分析した場合,群Aと群Bの平均的な値の行動の傾向の違いは推定されているわけで,母集団においてこういう介入を行った場合にはこういう行動がこれくらい増える(減る)といったことは因果効果の推定から考えていることだと思うが違うのだろうか。また,”行動の説明”と”行動の予測”というという言葉の使い分けの意味はなんだろうか。行動が説明できたら予測もできると考えるのが自然ではないだろうか。また,結局細かい因果効果の説明を積み重ねていけば,ある種の行動の予測ができるようになってもよいと思うが,一度だけの実験結果と統計モデリング(とこの論文で読んでいるもの)を比較するのはミスリーディングのように思う。

ここからわかるように,確率モデルを仮定したうえでパラメータを推定することで,理論的に知りたい値の推論や予測が可能になる。

さきに示した2群の平均値差の話でも,群AかBかによって大雑把な値は予測できるのだからよいのでは?結局”伝統的な分析手法”といいながら,最初の方で除外した検定論を仮想の敵としたいつもの論法になっているように思われる。

4.ベイズ推定による統計モデリング

パラメタ推定方法とモデリングは別の話だから,ならべてお話をするのは変な気がする。あくまでどのようにパラメタを推定するのか,という話であればわかるが。

4.1 汎用的な推定アルゴリズムである MCMC

対数尤度の最大値を求める問題にだって汎用的なアルゴリズムはあるし,Mplusといったソフトや,Rのnlme関数などを使えばOKだと思う。MCMCだけが汎用的なパラメタ推定方法であるかのように述べるのは不自然。

最小二乗法や最尤法でもうまく解けないという問題がしばしば生じる

というのは,数値計算アルゴリズムの問題であってベイズ推定だから解消できる問題というわけではない。たしかにMCMCを使うことで計算量の問題が比較的現実的になったり,局所最適解などの問題に悩まされないですむことはあるが,MCMCだって極めて相関が高いパラメタがあった場合にはパラメタ空間を自由に動けなかったり,ポテンシャル障壁の壁に阻まれて双山の分布の片方からしかサンプリングができないという問題もある。また,mixtureモデルのラベルスイッチング問題や収束診断がどれほど正しいかもわからない。そのような問題があるにも関わらず,MCMCだけが最高という論調には到底賛同できない。

この長所により,数学的に高度な演算から標本分布の導出などをせずとも,パラメータの推定ができるようになった。

最尤推定法の文脈でも漸近論があるので,パラメタの推定量の分布が正規分布になると考えてよいこともある。たしかに,パラメタの推定量が漸近的に正規分布に従うというのは,近似的な方法であり,正確な標本分布が導出できているわけではないが実用的にはよく使われている。

9) なお,異分散を仮定した平均値の差の推定については,構造方程式モデルなどを利用すれば最尤法で解くことが可能である。

差の推定???全体的に文章が正確で無い気がして読んでいていちいち引っかかってしまう。

したがって,MCMCでは最尤法などに比べて,パラメータを求めるためのアルゴリズムを容易に,汎用的に自動化できるのである。

今どき自動微分もあるし,尤度関数を作ればパラメタを推定してくれるソフトだっていくらでもあると思う。マルチレベルモデルなど変量効果項を含むモデルにおいてMCMCを使うことが有効な理由は,アルゴリズムが自動化できるからではなく,数値的に厳しい積分計算をうまく回避できるからだと思う。

しかし,すでに述べたように,階層モデルを最尤法で解くのは簡単ではない。それは,頻度主義統計学がパラメータを固定値と想定することから,パラメータに確率分布を仮定することが頻度主義の理論的枠組みからは困難だからである。

まったく違う。マルチレベルモデルモデルではパラメタの分布を考えているし,そのパラメタの推定方法として最尤推定法を使うことも当たり前に行われている。しかしパラメタ推定が困難になるのは積分計算の大変さとデータが少ないことによって,変量効果項の分散の推定値が0になってしまうというような問題である。たしかにMCMCを使って事後分布の平均を推定値とすればこうした問題は避けられるし,繰り返しの計算のなかでランダムに生成された分散がきっちり0という値を取ることはありえない。かなり小さい正の値をとることはあるが,log0といった計算が発生することはない。これによって,一見して不適解は回避できる。また,ベイズ的にパラメタの事前分布を想定したとしてその事前分布のパラメタ(ハイパーパラメタ)は定数だと思う。そうだとすると,ベイズであってもパラメタが定数であることはある。

4.4 最尤推定ベイズ推定の違い

渡辺先生の最尤推測とベイズ推測の違いなどをきちんと読んだほうが良い気がする。

そのようなときには,最尤推定のように点推定を探索するアルゴリズムを用いる方法(たとえば EM アルゴリズム)では,初期値によって推定値が変わってしまうという,初期値依存性が生じ,一意に解を得るのが難しい。

理論的にはMCMCは無限にサンプルを続ければ正確な事後分布を再現できるが,現実的にはそうは行かない。その結果として,局所からしかサンプリングしていなかったりなんだりかんだりがあるし,初期値の選択はMCMCであったとしても重要だと思う。また,”一意に解を得るのが難しい”という言葉の意味もちょっと不明。MCMCで定常分布からのサンプリングがなされたとしても,サンプリングする個数で推定値が微妙に変わったりするわけでMCMCだからといって一意に推定値が定まるわけではない。

心理学でベイズ統計モデリングが活用できる点として,以下の 3 つの可能性を挙げる。ひとつめは,行動データの生起メカニズムを確率モデルで記述することで,理論の精緻化と行動の予測に役立つ可能性である。もうひとつは,心理的な特徴を測定・推定するときに,より理論的に妥当になる可能性である。

といっているが,これらはすべて”統計モデリングによる利点”だと思う。ベイズ的に,パラメタの事後分布の一部を推定値に用いることの利点ではないと思う。また,測定・推定を並列に書いているが,全く別のものなので並べて書くことには違和感がある。


具体例はすっ飛ばして,

推定の区間推定ができる点も,ベイズ統計モデリングの利点

とあるが,推定の区間推定ってなに?と思ったのと,区間推定ができることは何も事後分布をパラメタの推定に用いたときのみの利点ではない。

モデルを書いたプログラムのコードの共有が簡単であることも意味する

とあるが,その理由がわからない。RなりPythonなりMplusなりSPSSなりSASなりStataなりなんでもいいけど分析に使ったコードをそのまま共有すればいいのでは???あとの段落でRだとモデルと推定が分離できないとか書いてる気がするが,応用系の人が自分で推定プログラムを書くことなんてほぼ無いわけで,例えばRのlm関数に何を書いてあるかとか見れば報告されてるモデルと分析で想定したモデルが正しいかどうかは見てわかると思うのだけれど。あと,検定批判の部分は今回の論文の本筋ではないといいつつ,結局再現性などの問題などは検定批判の文脈で出てきたことな訳で,しっかり議論を入れているので,議論が混乱する。統計モデルを推奨することと,再現性の問題は別問題だと思うのだが私の誤読だろうか。

その他気になったこととして,関数の形がわかるほど精緻に研究されている研究をより精緻にしていく研究も大事だけれど,新しい現象を発掘していくような研究も大事だと思う。なんか出てきたモデリングの例がその関数の形がわかったからどうなんだろう,というような印象がある。ほぼ関数型がわかっていて,それを階層化して当てはまりが良くなったとしても,そりゃ個人差あるから当然だよねくらいの気持ちしか無いし,あまり面白みがよくわからない。検定の文脈で変量効果項を無視するとタイプIエラーが増えるというのは大問題だが。なんか,出された例が魅力的に見えない。


めんどくさくなったのであとは省略するが,気になった点としては以下の点。

  1. ベイズ統計(ベイズ推定)=MCMCという図式。
  2. 用語の使い方があまり正確でない事によるよくわからなさ。
  3. 統計モデリングベイズ推定の関係のなさ。MCMCを使えば不適解を出さないからよい的な感じに見える。結果的にMCMCで適当にパラメタ推定しておけばOKというような感じ?
  4. 再現性の問題とモデリングとの関係のなさ。
  5. 確率的プログラミング言語モデリングの関係のなさ。
  6. 結局,ベイズ的なものを押すために検定批判的な紋切り型になっている点

仕事の仕方メモ

仕事の仕方について自分なりやっていること工夫など自分のためにメモを残しておく。言語化することによって,より精緻化したり,不要な部分を削除したり工夫ができるようになると思われる。

論文を書く手順

  1. イデアを思いつく
  2. 先行研究を調べる(このときに関連論文をまとめておくとそのままイントロに使えるので,その手間を惜しまないようにしておく)
  3. 先行研究を参考にしながら,新しいアイデアを形にまとめる。特にcontributionについて意識をしながら,何が主張できるのかを言語化する。
  4. プログラムを実装する(必要に応じてシミュレーション,実データ解析を行う)。
  5. 結果の図表のイメージを作る。このとき,自分の主張に適切な形になるように留意。ただし,簡単な条件だけでなく,難しい条件なども含めておく。
  6. 論文全体のアウトラインを作る。とくにイントロのストーリーをきちんと考える。あとで目的が大きく変わらないように,目的の部分はぶれないように注意。
  7. イントロを書きながら引用文献を追加してく。
  8. 理論的な部分(数式を書いたりなんだり)は淡々と作業として行う。
  9. シミュレーションなどの条件,結果,考察も順に淡々と書く。
  10. 一旦ドラフトを書き上げることを目指す。
  11. ドラフトが書けたら,全体の推敲や記述を追加するべきところなどを考える。
  12. ほどよいところで,校閲に出す
  13. 校閲の結果を確認して良さそうであれば,投稿する。

研究全体の捉え方・研究の実行方法について

  1. 論文執筆やアイデアの具体化をする際には,何も考えなくても実行できるくらい細かくタスクを分割する(アルゴリズムに落とし込めるように)。
  2. 簡単でつまらない作業を積み重ねるといつの間にか一連の塊になっているようにする。
  3. 大前提として,論文執筆はつまらないしアイデアを考えることも面白い作業であるわけではないことを受け入れる。
  4. つまらない作業を積み重ねることが研究だと割り切る。
  5. 作業それ自体はつまらないかもしれないが,研究をすることで,自分の人間性を磨き,世の中に貢献することができることを理解する。
  6. 研究と執筆は習慣(歯磨きと同じ)。やる気があるときにやるのではなく,やる気とは関係なくやる。とにかく書く。考えることは書くことと同じ。
  7. 面白いからやるのではなく,解かなければならない問題があるからやる。
  8. 最も大事なことの一つが,解くべき問題を同定すること。この問題設定の良し悪しによって,研究の質が決まってしまう。
  9. アドバイスを求める相手を常に持つ。
  10. 現象を捉えるために言語が必要になるので,言葉に敏感になる。
  11. 効率的な論文執筆は難しいが,非効率な論文執筆を避けることはできる。非効率な論文執筆方法は複数の論文を掛け持ちして書くこと。プロジェクトは一つに限定して,論文投稿するまでは基本的に他のプロジェクトを開始しない。
  12. 一つの論文のネタに書けていい時間は3~4ヶ月程度を目安にする。時間がかかりすぎないように手早くまとめる。
  13. タイムスパンとしては,テーマを決めることとアイデアを形にして行けそうかどうかを判断するのは1~2ヶ月以内で,論文執筆は2ヶ月以内程度で行う。
  14. ソフトやPCや機材にはあまりこだわらないが,結局認知資源と手間を削減できるものを選択する。正直Word,Excel,Rをきちんと使えれば,バージョン管理ソフトは実質的に必要にならない気がする。TeXは認知資源が必要すぎる。プレゼン資料もbeamerよりPower Pointで作ったほうが,疲れない。
  15. 執筆は始めるのが一番めんどうなので,何も考えずにソフトを立ち上げ,”ああああ”など打ち込んで,書き始めた気になることが大事。
  16. 気分が優れなかったら,すぐに寝る。昼寝大事。
  17. 寝不足だと何もできないので,寝る。とにかくすぐに寝る。

研究のアイデアの出し方のヒント

  1. イデアを思いつくためには,現状と理想のギャップに注目する
  2. モデルを拡張するためには,パラメタを足す,データの相を増やす,今まで無視されていた項を入れる,モデルを混ぜる,
  3. 推定アルゴリズムを作ることも研究
  4. 既存のモデルを比較することも研究
  5. 実データに最先端のモデルをいくつか適応してみることも研究。このときに,特にどのような困りが生じるのか,狙ったとおりにならない点を注意深く観察する。”もっとこうだったらいいのに”という点を感じる。
  6. 言葉の定義をキチンと自分なりに消化する。研究者が必ずしも,同じ言葉を同じ意味で使っているとは限らない。
  7. 定式化が意味不明なものを理解できる形にするのも大事。拡張の可能性を探る。
  8. 既存のモデルにこだわり過ぎない。そもそもなんのためにそのモデルが存在したのかを調べる。
  9. 自分がしっくりこないものを大事にする。どういうふうに表現するとしっくり来るのか,気持ちのよい表現を考える。
  10. 他分野の理論研究ですでに利用されている知見が利用できないか考える。学問がつながっていることを意識する。
  11. 大きな文脈で考える。すでに流行しているものは決して追わない。誰かがやるから,自分がやれるようになるころには,研究のネタはなくなる。
  12. 地味だけど大事な研究もあることを意識する。
  13. 誤解,誤用されているモデルについて考える。
  14. 疑問や興味をひたすら言語化する。異なったことでも意外とつながっていることがあるので,自分の真の興味に敏感になる。
  15. テクノロジーの発展に注目する。新しく得られるデータや分析方法などがあるかどうか。
  16. 誤った問題設定を行わないように気をつける。意味がないと感じられることは決してやってはいけない。自分にとって意味があっても,他人に意味がないと思われることもやってはいけない。これは,ウケが良い研究を行うということを意味しているのではない。意味がある研究ができるかはわからないが,意味がない研究を避けることは絶対に必要。

Stanでマルチレベル(マニュアルコピペ)

参考:
https://mc-stan.org/docs/2_18/stan-users-guide/multivariate-hierarchical-priors-section.html
1.13 Multivariate Priors for Hierarchical Modelsより
高速なコードにはなっていないが,これで基本的にはOK。
添字が慣れたものではないが,読めばわかるような感じになっている。
cross classifiedモデルへの拡張や,適当なlink functionを仮定して従属変数が連続量以外でもOKな状況に持っていけるはず。

#
# Stanでマルチレベル lkj分布を使ったりなんだり
#
library("rstan")

# データ
library(haven)
mlmdata <- read_dta("https://stats.idre.ucla.edu/stat/examples/imm/imm10.dta")
write.csv(mlmdata,file="imm.csv",row.names = F)
dat <- read.csv(file = "imm.csv")
head(dat)

# データの用意
N <- nrow(dat)
x <- dat[,c("ses","homework","sex","white")]
K <- ncol(x)
J <- length(unique(dat$schid))
u <- dat[,c("meanses","public" )]
u <- u[!duplicated(u[,1]),]
L <- ncol(u)
y <- dat$math

uniq_sch_id <- unique(dat$schid)
jj <- rep(NA, N)
for (j in 1:J){
  jj[dat$schid==uniq_sch_id[j]] <- j
}
dat_stan <- list(N = N,
            K = K,
            J = J,
            L = L,
            jj = jj,
            x = x,
            u = u,
            y = y)
# int<lower=0> N;              // num individuals
# int<lower=1> K;              // num ind predictors
# int<lower=1> J;              // num groups
# int<lower=1> L;              // num group predictors
# int<lower=1,upper=J> jj[N];  // group for individual
# matrix[N, K] x;               // individual predictors
# row_vector[L] u[J];          // group predictors; uはJ×Lの行列でデータを指定する。
# vector[N] y;    
help(package="rstan")

cat("data {
  int<lower=0> N;              // num individuals
  int<lower=1> K;              // num ind predictors
  int<lower=1> J;              // num groups
  int<lower=1> L;              // num group predictors
  int<lower=1,upper=J> jj[N];  // group for individual
  matrix[N, K] x;               // individual predictors
  row_vector[L] u[J];          // group predictors
  vector[N] y;                 // outcomes
}
parameters {
  corr_matrix[K] Omega;        // prior correlation
  vector<lower=0>[K] tau;      // prior scale
  matrix[L, K] gamma;           // group coeffs
  vector[K] beta[J];           // indiv coeffs by group
  real<lower=0> sigma;         // prediction error scale
}
model {
  tau ~ cauchy(0, 2.5);
  Omega ~ lkj_corr(2);
  to_vector(gamma) ~ normal(0, 5);
  {
    row_vector[K] u_gamma[J];
    for (j in 1:J)
      u_gamma[j] = u[j] * gamma;
    beta ~ multi_normal(u_gamma, quad_form_diag(Omega, tau));
  }

//  for (n in 1:N)
//    y[n] ~ normal(x[n] * beta[jj[n]], sigma);
// Speed UP!!!
 {
    vector[N] x_beta_jj;
    for (n in 1:N)
      x_beta_jj[n] = x[n] * beta[jj[n]];
    y ~ normal(x_beta_jj, sigma);
  }
}
",file="multilevel.stan")

fit1 <- stan(
  file = "multilevel.stan",  # Stan program
  data = dat_stan,    # named list of data
  chains = 2,             # number of Markov chains
  warmup = 1000,          # number of warmup iterations per chain
  iter = 2000,            # total number of iterations per chain
  cores = 2,              # number of cores (using 2 just for the vignette)
  refresh = 1000         # show progress every 'refresh' iterations
)

stan_trace(fit1, pars=c("gamma"))
stan_trace(fit1, pars=c("Omega"))
stan_trace(fit1, pars=c("tau"))
stan_ac(fit1, pars=c("gamma"))
stan_ac(fit1, pars=c("Omega"))
stan_ac(fit1, pars=c("tau"))
stan_hist(fit1, pars=c("gamma"))
stan_hist(fit1, pars=c("Omega"))
stan_hist(fit1, pars=c("tau"))

print(fit1, pars=c("gamma", "Omega", "tau"), probs=c(.1,.5,.9))
plot(fit1, pars=c("gamma", "Omega"))

以下のチュートリアルも非常に参考になる
https://web.stanford.edu/class/psych201s/psych201s/papers/SorensenEtAl.pdf