黒木玄 Gen Kuroki(@genkuroki)/2020年03月08日

#統計「ベイズ統計モデリング」の勉強をしたくてかつ #Julia言語でお気楽に数値計算(two language problem 無しに高速計算可能)を楽しみたい人には、 Turing.jl は結構おすすめかも。 Julia v1.4.0-rc2.0 で Turing.jl を使うには以下のリンク先の情報が必要。

twitter.com/genkuroki/stat...

タグ： Julia言語統計

posted at 23:37:14

#統計 #Julia言語添付画像は既出の

nbviewer.jupyter.org/gist/genkuroki...

より。絶対温度ゼロのベイズ統計が最尤法に一致することのガンマ分布モデルにおける数値的確認。 pic.twitter.com/vxWcxJ6RDS

タグ： Julia言語統計

posted at 23:34:28

#Julia言語 Turing.jl でも任意の対数尤度函数を使用できます。大体において、確率分布の型を自前で定義して、logpdf (など)を定義しておけばよいようです。以下のノートはその応用例。

nbviewer.jupyter.org/gist/genkuroki...
逆温度入りの対数尤度函数のTuring.jlでの使用法

これを使えばWBICも計算可能。

タグ： Julia言語

posted at 23:27:36

#統計 20世紀に「頻度主義vs.ベイズ主義」の論争があったようだが、「虐げられていたベイズ統計が現代において復活した！」の類の解説は誤り。

実際には、ベイズ主義陣営も主観確率云々のような非合理的な考え方をしていて、非難されて当然のことを言っていた。

論争の双方がダメな20世紀の黒歴史。 twitter.com/genkuroki/stat...

タグ：統計

posted at 21:09:36

#Julia言語数値積分の精度を rtol=1e-3 まで落としたら、18秒ちょっとかかっていたのが、15秒ちょっとまで短くなった。積分区間を縮めてもよさそうだが、面倒なのでこの辺でやめにしておく。 (MCMCの結果はどうせ誤差が混じるので精密な数値積分には大して意味がない。)

twitter.com/genkuroki/stat...

タグ： Julia言語

posted at 20:58:50

#Julia言語訂正。数値積分の範囲が狭かったことに気付いて、やり直した。

数値積分無して「公式」を使う方法ならば0.24秒で済むWAICの計算が、18秒以上かかった。7×2万回=14万回数値積分が実行されている。

訂正：13秒→18秒

twitter.com/genkuroki/stat...

タグ： Julia言語

posted at 20:50:14

僕はよく学生さんに「ヒストグラムはモデルであってデータではないのだよ」と口癖のように言っています。

タグ：

posted at 20:37:38

#統計赤池弘次さん本人によれば、AICのAは "an" を意味します。AIC = an information criterion なのだそうです。証拠は

www.jstage.jst.go.jp/article/butsur...
エントロピーとモデルの尤度(<講座>物理学周辺の確率統計)
赤池弘次
1980

のp.612にある(添付画像)

twitter.com/jp011100/statu... pic.twitter.com/s8B6BPE0Gy

タグ：統計

posted at 20:37:00

階級幅を変えてプロットしたり、場所を変えたりと無限の選択肢がある中でデータの見え方の構造が変わるという体験で、統計的推測がそもそも不良設定問題であり、ヒストグラムはデータではなくモデルなのだ、と認識することは大切に思います。

タグ：

posted at 20:36:33

モデル分布や分布そのものについて考えるには、ヒストグラムモデルを最初に学ぶのが良さそうだなと思っています。Sanovの定理を考えるときも。個人的には神経スパイクヒストグラムの幅をどのようにとるか、という問題でこのことを考える機会を得ました。

タグ：

posted at 20:29:09

英語版は、1986です(※注意：手書きのメモで1987と書いてしまっていますが、1986の間違いです)
twitter.com/katzkagaya/sta...

タグ：

posted at 20:25:25

情報量統計学はこちらの本です。twitter.com/katzkagaya/sta...

タグ：

posted at 20:21:52

#Julia言語スレッド違いを回収

Julia v1.4.0-rc2.0 で Turing@0.9.0 を動かしたい人は Libtask.jl のインストールで工夫する必要がある。(私の環境特有の問題の可能性もありますが。誰かがそうであることを確認してくれたら、このスレッドでコメントします。)

twitter.com/genkuroki/stat...

タグ： Julia言語

posted at 20:21:08

#統計 n = 7 でMCMCのchainの長さが 20000 なので、WAICを計算するために数値積分を7×20000回行っている。計算に13秒ちょっとかかった。

タグ：統計

posted at 20:17:03

#統計 #Julia言語添付画像を見ればわかるように、階層ベイズモデルも Turing.jl でシンプルに記述できます。

全部 Julia なので、WAICを計算するために数値積分が必要になる場合には Julia の数値積分パッケージを普通に使用することになります(配列 Y にサンプルを入れて渡す)。 pic.twitter.com/gRbfx9ePZg

タグ： Julia言語統計

posted at 20:15:03

#統計 #Julia言語

簡単な例題である The Seven Scientists の例題を Turing.jl で解いてみました。4種のモデルで推定して、WAICを比較しています。

添付画像は階層ベイズのWAICを数値積分を使って計算しているところ。

nbviewer.jupyter.org/gist/genkuroki...
測定スキルの異なる7人の科学者たち (Turing.jl版) pic.twitter.com/rklj8X6Jaa

タグ： Julia言語統計

posted at 20:10:56

積分定数 @sekibunnteisuu

@genkuroki @sekibunnteisuu @metameta007 @sunchanuiguru @kuri_kurita あなたのツイートをトゥギャりました。読んでいただけると幸いです。 togetter.com/li/1478635

タグ：

posted at 19:13:11

Masa Yamamoto予測誤差が大き @mshero_y

#Julia言語 v1.4.0-rc2.0 で Turing.jl を使う方法。

]add Libtask#358fa5637233aab396644b0945252e65630b9a07
build Libtask
add Turing@0.9.0

#49 のコミットまで戻したら、build Libtask に成功した。

github.com/TuringLang/Lib...

タグ： Julia言語

posted at 14:54:46

これ本当に試すべき。いろんなこと確認できる。 twitter.com/genkuroki/stat...

タグ：

posted at 14:22:57

#統計本を読むだけはAICやBICについて理解できない。それらを実際に使うのはn→∞ではなく、有限のnでなので、コンピュータなどで実際に計算してみないと様子がわからない。

大数の法則や中心極限定理でも同じ。

n→∞での漸近論は理想化されたフィクション。

有限のnでの様子の予備知識が必須。

タグ：統計

posted at 14:17:53

#統計 AICやBICなどを理解するためには、コンピュータで

(1)真の分布q(x)に従う乱数列(=サンプル)を生成

(2)そのサンプルを用いて、モデル0とモデル1のAICやBICなどを求める。

(3)以上の(1)と(2)を100~1万回繰り返して、AICやBICによるモデル選択が確率的にどのように揺らぐかを確認。

タグ：統計

posted at 14:12:53

#統計ベルヌーイ分布モデルでのBIC、自由エネルギー、WBIC

やはりグラフ中で0未満になると真のモデル選択に失敗する。

しかし、n→大で左右のモデル選択失敗領域がどんどん狭くなって行く(確率が0に収束する)。

twitter.com/genkuroki/stat...

タグ：統計

posted at 14:08:23

#統計ベルヌーイ分布モデルでのAIC, WAIC, LOOCV

グラフ中でそれらが0未満になると真のモデル選択に失敗する。

x軸はサンプル。サンプルが生じる確率でx軸を非線形にスケールしてある(プロットでがんばった！)。

サンプルが左右に大きく偏るとモデル選択に失敗する。

twitter.com/genkuroki/stat...

タグ：統計

posted at 14:04:39

#統計「AICとBICのどちらが優れているか」という問いは不毛。目的が違う道具でかつ性質も違うだけ。

20世紀の統計学史はそういう不毛な問いに関するくだらない論争の歴史(黒歴史)なので要注意。20世紀の統計学における論争は本当にひどいものが多い。あれを肯定的に後世に伝えてはいけないと思う。

タグ：統計

posted at 13:59:18

#統計 BICは採用した確率モデル内(ベイズのモデルを考える)での標本分布がどれだけ真の標本分布を近似しているかの相対比較の推測のための道具です。予測分布より、モデル内標本分布の方が扱いが易しい。BICなら対数尤度比検定の状況でもサンプルサイズ→∞でモデル選択に失敗する確率は0に収束する。

タグ：統計

posted at 13:55:45

#統計 AICは予測分布の予測精度の大小の確率的な推測を行うための道具であり、予測分布の予測精度は数学的に難しい対象なので、サンプルサイズ→∞としても予測精度が低い側を選んでしまう確率は0にならない。

タグ：統計

posted at 13:51:55

#統計対数尤度検定の状況での対数尤度比に関するχ²検定は、AICの差がある値を超えたときに帰無仮説側のモデルを棄却することと同じ。対数尤度比検定は、帰無仮説側のモデルが選択され易いように下駄を履かせた場合のAICによるモデル選択と同じ。

対数尤度比検定が使えない状況でもAICは使える。

タグ：統計

posted at 13:49:02

#統計罰則項を、定数から O(log n) に厳しくしたことによって、BIC では対数尤度比検定の帰無仮説が成立する状況であってもモデル選択に失敗する確率は n→∞ で0に収束してくれます。

タグ：統計

posted at 13:40:48

#統計 BICではパラメータ数を増やすことの罰則項は O(log n) のオーダー:

BIC = -2×最大対数尤度 + パラメータ数×log n.

一瞬、ぎょっとするかもしれませんが、最大対数尤度の項は O(n) のオーダーなので、O(log n)の罰則項の影響はサンプルサイズ n を大きくすると相対的に小さくなります。 twitter.com/jp011100/statu...

タグ：統計

posted at 13:37:31

#統計ある意味でAICはモデル選択が機能するぎりぎり小さな罰則項を採用している。

だから、対数尤度比検定の帰無仮説が成立する状況でAICでモデル選択をすると、n→∞としてもモデル選択に失敗する確率が0に収束せずに、正の値に収束してしまいます。

これはAICユーザーにとって必須の予備知識。

タグ：統計

posted at 13:33:21

#統計

AIC = -2×最大対数尤度 + 2×パラメータ数

の最大対数尤度の項はO(n)のオーダーで、パラメータ数は固定されているので、サンプルサイズn→∞では最大対数尤度の項が支配的になります。

AICでは、nを大きくすると、パラメータ数を増やすことへの罰則項の影響が急激に小さくなります。

タグ：統計

posted at 13:29:31

#統計サンプル (X_k, Y_k), k=1,…,n の確率モデル p(y|x,w) の対数尤度は

log p(Y_1|X_1, w) + … + log p(Y_n|X_n, w)

は同程度の大きさの量のn個に和になっているので、O(n) のオーダーの量になります。最大対数尤度では w が動くのですが、同様に O(n) のオーダーの量になります。

タグ：統計

posted at 13:25:48

#統計

AIC = -2×最大対数尤度 + 2×パラメータ数

と2倍する御利益は、統計学を学んだ人にとって馴染みのあるχ²分布のスケールになることです。

対数尤度比比のχ²検定ではχ²分布のスケールに合わせるために、対数尤度比を2倍します。

この習慣には便利な所があるので私もよく使っています。

タグ：統計

posted at 13:14:48

#統計これは鋭い。

2である必要はないです。

AIC = -1×最大対数尤度 + 1×パラメータ数

としても、AICは単に大小比較で使われるだけなので、モデル選択には影響しません。渡辺澄夫『ベイズ統計の理論と方法』では

AIC = -(1/n)×最大対数尤度 + (1/n)×パラメータ数

となっている。続く twitter.com/jp011100/statu...

タグ：統計

posted at 13:11:16

巨大な集団中の小さな割合の人数

が

相対的に小さな集団そのものの人数

よりも大きくなることがあるという易しい話なのに、なぜか難しく聞こえるように説明することが行なわれている。

タグ：

posted at 12:56:16

ベイズの定理を使って理解するべきだと思っている人達は、高校数学レベルで「ダメな人達」です。

今回の件で大事なことだと思って「ベイズの定理」について勉強して応用しようとした人達は、伝統的に出回っている「標準的な解説」の被害者。

問題は被害者が被害者を増やす流れになっていること。 twitter.com/genkuroki/stat...

タグ：

posted at 12:45:54

関連

twitter.com/genkuroki/stat...

タグ：

posted at 12:41:41

仲間内ではそういう特殊な専門用語を使ってもよいと思うが、そうでない場合には、「その手の用語を使って難しく聞こえるように説明する必然性は皆無であること」を理解しておかないと、陽性的中率の類の話をベイズの定理を使って説明するレベルが低い迷惑な人達を増やすことになってしまう。

タグ：

posted at 12:38:02

#統計最近知ってあきれたのが、高校数学レベルの易しい話(実際に大学入試問題でも出題されている話)なのに、「感度」「特異度」「事前オッズ」「事後オッズ」「陽性尤度比」のような新たな用語を導入して、さらに「ベイズ」とか言い出す人達が沢山いること。医療統計という分野の話。

タグ：統計

posted at 12:33:56

#統計統計学用語は全般的に権威的に響くものが多くて好きになれない。普及してしまったので仕方無しに使っている。

統計学の応用=ギャンブルを勧める人が学問的に権威を持っているように見えるのはよくない。「こいつは確率的に失敗する方法を勧めて来る危ないやつだ」という扱いが好ましいと思う。

タグ：統計

posted at 12:22:40

☔️ @nardtree

iPadで使えるterminalソフトのblink.sh ってのが凄くて、mosh使えるし、なんかやたら応答いいし、日本語もバグらず入力できるし、256色対応なので、MacBookPro持ち歩かなくなった😊
肩こり減った😊 pic.twitter.com/DStCh7TPb4

タグ：

posted at 12:09:38

#Julia言語ブロードキャストを通したくない配列は Ref() で囲む。これを知らないと、ちょっと面倒な書き方を強いられる。

mean(a .≤ x) と書くと a .≤ x の時点で配列が作られてしまう点が気になるなら

mean(s ≤ x for s in a)

と書く。(s ≤ x for s in a) は配列ではなく、generator になる。

タグ： Julia言語

posted at 11:52:37

#Julia言語解説

a .≤ x で成分ごとに x 以下なら true にそうでないなら false に変換。mean(a .≤ x) で変換後の平均を取ると、true, false を 1, 0 に変換してから平均を取るので、配列 a の成分中で x 以下のものの割合が求まる。

Ref(sample) でブロードキャストから保護される。 pic.twitter.com/YGfoTHMfdn

タグ： Julia言語

posted at 11:52:36

#Julia言語私もやってみた。

Base.Fix1 とか Base.Fix2 というのもある。

@. は多用している。 twitter.com/negi__/status/... pic.twitter.com/IvVaJ6pQna

タグ： Julia言語

posted at 11:41:29

しぶてぃ @takuizum

#Julia言語工夫をすれば，CSVでなくてもテキストデータ(.jlとか.stanでも)をstringとして読み込める。

タグ： Julia言語

posted at 11:15:24

しぶてぃ @takuizum

(｢・ω・)｢ｶﾞｵｰ @bicycle1885

#Julia言語でWeb上のCSVデータをデータフレームとして取り込みたいとき（GitHubに上がっているサンプルデータとか）は，
using HTTP, CSV
res = HTTP.get("https://hoge.csv")
data = CSV.read(IOBuffer(res.body))
で行けると思う。もっと効率良くできるかもしれない

タグ： Julia言語

posted at 11:15:24

非公開

タグ：

posted at xx:xx:xx

Julia全く分からんので実家の笛職人継ぐしかないな。

タグ：

posted at 04:06:33

非公開

タグ：

posted at xx:xx:xx

Josef Heinen @josef_heinen