黒木玄 Gen Kuroki(@genkuroki)/2017年10月11日

#数楽

 #JuliaLang で新たな確率分布を定義し、それをMCMCで利用する方法
以下の後者は最新版。前者にも直に反映されるはず。画像に続く

nbviewer.jupyter.org/gist/genkuroki...
gist.github.com/genkuroki/4653...

タグ： JuliaLang 数楽

posted at 00:38:16

#数楽 #JuliaLang

外れ値のあるサンプルの線形回帰

1つ目の画像
単純な最小二乗法を使うと右上の外れ値に影響されて傾きが大きめになる

2つ目の画像
残差～t分布×定数で回帰。ベイズ推定の予測分布はうまく外れ値を処理している感じ。右上の外れ値はプロットされていない pic.twitter.com/FvulimSSn9

タグ： JuliaLang 数楽

posted at 00:42:18

#数楽 #JuliaLang

本当は最小二乗法の方も予測分布をプロットするべきだったかも。
そうしないと同じものを目で比較できない。

タグ： JuliaLang 数楽

posted at 00:43:20

#数楽 #JuliaLang

二つ前のツイートの添付画像のケースでの事後分布の様子

y = a + bx + 残差

a,bの最小二乗法(LSM)での推定結果(左の方の黄色い四角)とt分布×定数モデルでの推定結果は異なる。 pic.twitter.com/CHvhyIctyZ

タグ： JuliaLang 数楽

posted at 00:47:19

#数楽 #JuliaLang

別のサンプルの例

こちらは、書くべき数値と異なる数値を記載したことが原因で生じた外れ値の処理

残差＝定数×t分布のモデルをベイズ推定で解いた結果(最終結果は、事後分布ではなく、予測分布だと考えた方がよい)はうまく外れ値を処理している。 pic.twitter.com/RWoLEWrD8Z

タグ： JuliaLang 数楽

posted at 00:50:06

@yuzu_96wiz_only 【1本ずつならば計30本、2本ずつならば計60本。誤りではない】
という考え方で　３０×２　とする事はOK、というご意見ですね。
では、
「３×２だと、耳が３本生えたウサギが２羽、ということになるよ」というのは誤りですか？ #掛算
 www.asahi.com/edu/student/te...

タグ：掛算

posted at 00:51:41

#数楽事後分布を眺めるだけで、(事後)予測分布の様子がよくわかるわけではないので、事後分布から予測分布を計算して、何らかの方法でプロットするべきだ。

「パラメーターの推定」という発想は中途半場で誤解の原因になる。「確率分布の推定」というシンプルな考え方の方がよいと思う。

タグ：数楽

posted at 00:53:19

#数楽階層モデルでないケースについては、予測分布の計算は簡単で、計算も比較的軽い。しかし、階層モデルでは数値積分したものをさらに大量に足し上げる計算が必要になる。この問題をうまく解決しないと、階層モデルでのWAICの計算で困ってしまうことになる。

誰か解決して！😆

タグ：数楽

posted at 00:55:40

#数楽 WAICは予測分布が未知の真の分布をどれだけの精度で予測しているかを比較するときの指標の一つです。たぶん、LOOCVよりも優れている。

問題は階層ベイズモデルの場合のWAICの(ケースバイケースではない自動的な)計算。テキトーにやると計算量が爆発する。

タグ：数楽

posted at 00:57:47

#数楽訂正「中途半場」→「中途半端」

ベイズ統計に関する誤解は「パラメーターの推定」という発想から出て来ていることが多いと思う。

twitter.com/genkuroki/stat...

タグ：数楽

posted at 01:04:03

非公開

タグ：

posted at xx:xx:xx

#数楽「パラメーターの推定」という発想から抜け出せないから、「パラメーターも確率分布する」ということについて何か大げさに疑問を持ったり、解釈に苦しんだりしなければいけなくなる。最尤法とベイズ推定法は数学的には明らかに似たような方法であり、根本的に思想が異なるわけではない。

タグ：数楽

posted at 01:06:19

#数楽最尤法でパラメーターを推定するときには、ローカルマキシマムではない真のマキシマムを求まったかどうかを知るために、結局のところ尤度函数のグローバルな様子を見なければいけなくなります。ベイズ推定法で事後分布全体の様子を見なければいけないこととそう変わらない。

タグ：数楽

posted at 01:26:54

#数楽 #R

github.com/kosugitti/BSJ_...
行動計量学会第19回春の合宿セミナー　Aコース
今日から始めるベイジアンモデリング

Rの勉強をしたい人には役に立つかも。

タグ： R 数楽

posted at 01:33:39

PsycheRadio @marxindo

佐藤実｜『プリンセス・フィリシア　物理の @sato_minoru

「研究代表者は長年にわたり科研費を申請していなかったが、このたび本務校が科研費を申請しない教員にペナルティを課すことになり、科研費を申請せざるを得なくなった。それが本研究の着想につながった」。

タグ：

posted at 05:52:59

胸骨圧迫を絶え間なく続けることが最優先．気道確保や人工呼吸は可能ならやるが，そのために胸骨圧迫が続けられないならやらなくてもいい．とにかく，脳への血流を止めないことが重要，とのこと twitter.com/konamih/status...

タグ：

posted at 06:34:50

グレッグ @glegory

Dr. Chris Rackauckas @ChrisRackauckas

安倍政権前は「ゼロ金利だ。金融政策にもうできることはない。円高を甘受して、前提にした企業戦略を考えなければならない」なんて平気で言われていたんですよ。どうですか？この５年間、ほぼ金融政策しかしてません。この程度の回復はするんですよ。いかにマスコミや評論家がいい加減かわかりますねｗ

タグ：

posted at 11:51:13

非公開

タグ：

posted at xx:xx:xx

@ipnosimmia @randyzwitch @rmflight @dataandme And I'm programming it through VNC on my iPhone.

タグ：

posted at 14:56:11

しかしCVの理論解析はAICとかBICに比べて全然ないのだ...

タグ：

posted at 17:06:37

積分定数 @sekibunnteisuu

@yuzu_96wiz_only 「掛け算の単元だから出てきた数字をそのまま掛けた」となるのを避けるには、足し算や引き算や、あるいは「割り算」の問題も混ぜて出せばいいだけです。6個の蜜柑を3人で分けたら1人何個か？これは、割り算どころか足し算も学んでいなくても、数の概念を習得していれば答えられます。

タグ：

posted at 17:13:50

積分定数 @sekibunnteisuu

@yuzu_96wiz_only 掛け算の単元で、文章問題に2つの数値しか出てこなくて、掛ければ答えが求まるならような問題しか出さないなら、「掛け算の単元だから掛け算にする」というのは当然のことです。そこで、本来全く同じ意味の３×４と４×３を違うものと教えるなどと言うのは馬鹿げています。

タグ：

posted at 17:15:17

Takashi Miyamoto @tmiya_

@dekosuke 渡辺澄夫先生の教科書辺りの議論じゃダメなんでしょうか？

タグ：

posted at 17:16:53

@tmiya_ W理論そこらへん言及しているんですね。理解するにはどの本を読めばよいでしょう

タグ：

posted at 17:19:12

トデス子'\ @todesking

@dekosuke 「ベイズ統計の理論と方法」でCVとWAICの関連について解析されてました

タグ：

posted at 17:47:14

gezigezi_本体 @_gezigezi_

夜泣き原因を探るため、真冬に子供の布団に熱電対を入れて取ったデータなんだけど、最初に寝たあと沢山の波が現れるのは子供が暑くて寝返りを打つから。そして起きてしまい寝かしつけに１時間半かかってる。２時半以降気温が下がったので5時近くまで寝たがまた起きてしまった…という図 #夜泣き対策 pic.twitter.com/TW7BhniJLv

タグ：夜泣き対策

posted at 17:55:55

Takashi Miyamoto @tmiya_

@dekosuke コロナ社の「ベイズ統計の理論と方法」か、BDA3のWAICについて言及されてる辺りとかでしょうか。

タグ：

posted at 18:22:49

#数楽添付画像は

StanとRでベイズ統計モデリング読書会　Chapter 7（7.6-7.9）　回帰分析の悩みどころ　～統計の力で歌うまになりたい～
www.slideshare.net/nocchi_airport...

のp.45より。そのグラフから目で数値を読み取って私もやってみた。 pic.twitter.com/Ym1EMt8ofa

タグ：数楽

posted at 20:30:13

#数楽 #JuliaLang 続き。私もやってみた結果は以下のリンク先にある。

gist.github.com/genkuroki/0125...
カラオケにおける抑揚と点数の関係

slideshareの方のもとの分析でもt分布を使ったモデルをベイズ推定で解いています。続く

タグ： JuliaLang 数楽

posted at 20:32:36

#数楽 #JuliaLang GitHubでの閲覧が重いと感じたら、

nbviewer.jupyter.org/gist/genkuroki...
カラオケにおける抑揚と点数の関係

を利用して下さい。

そこで扱われているモデルの解(予測分布)は事前分布によって結果が大きく変わります。

タグ： JuliaLang 数楽

posted at 20:34:25

#数楽 #JuliaLang 残差が「t分布×定数」に従っているとするモデルの尤度函数を眺めたことがあれば、そのモデルの尤度函数はローカルマキシマムを持つ場合があったり、t分布の自由度が大きな方に分布が細く長く伸びている場合があることを知っているはずです。

タグ： JuliaLang 数楽

posted at 20:36:23

#数楽 #JuliaLang 「t分布×定数」モデルの尤度函数の形については

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

を見て下さい。添付画像のように尤度函数は左右反対のL型に伸びた分布になります。大きなνはt分布がほぼ正規分布になることを意味しています。 pic.twitter.com/MMCHrbBiOX

タグ： JuliaLang 数楽

posted at 20:40:43

#数楽「事後分布＝尤度函数×事前分布」なので、事前分布の設定によってモデルで使用しているt分布が正規分布に近付くことを許したり、抑制したりすることができます。「外れ値」の問題を扱っている場合にはこれが重要な問題になるはずなことが、尤度函数の形を知っていればわかる。

タグ：数楽

posted at 20:42:46

#数楽「ρ×(自由度νのt分布)」はν=∞では平均0分散ρの正規分布に一致します。ρは正規分布の分散に対応するパラメーターであり、νは正規分布からどれだけ離れるか(ν>0が小さいほど離れたことになり、たくさん離れるほどすそが重い分布になる)を表わすパラメーターです。

タグ：数楽

posted at 20:45:24

#数楽「外れ値」の処理は誤差項が正規分布に従っていると仮定するとフィッティングがうまく行かない場合に必要になります。正規分布から離れるためのパラメーターνを用意して、外れ値が出易い分布までモデルで試用する確率分布の範囲を広げれば、外れ値を自動的にうまく処理できる可能性がある。

タグ：数楽

posted at 20:47:19

#数楽しかし、現実にはそう簡単な話ではないです。変形パラメーターνをt分布の自由度だとみなせば、t分布を用いたモデルができあがるのですが、その尤度函数はνが大きな側に細く長く伸びている場合が多い。ベイズ推定でもνが大きくなることを抑制しないと、νが巨大になることがありえる。

タグ：数楽

posted at 20:49:43

#数楽 t分布の自由度νの値が巨大であると推定された場合には、モデル＋事前分布の組み合わせが「このサンプルは正規分布に従っている疑いが強い」と判断していることを意味しています。

もしもモデル＋事前分布の組み合わせが適切なら、これは別に悪いことではありません。

タグ：数楽

posted at 20:51:39

#数楽しかし、事前分布の取り方によってt分布の自由度が巨大になることを容易に抑制できるという事実は、「外れ値だとみなしたい数値を無視して直線を引くこと」を事前分布をうまく選ぶことによってやりやすくなることも意味しています。事前分布を自分好みの結論を出すために使える！😱

タグ：数楽

posted at 20:54:01

#数楽ネタもとの www.slideshare.net/nocchi_airport... では添付画像のようにt分布の自由度νの事前分布を設定しています。

ν = 1 + 5 T(4)、 T(4)～(自由度4のt分布)

これの確率密度函数が2つ目の添付画像の青の実戦です。橙の破線は正規分布。 pic.twitter.com/TWqzb2iUb6

タグ：数楽

posted at 21:09:28

#数楽続き。さらにν>1と仮定しているようなので、この設定はνは0から30程度までのあいだに分布することにしてしまっているのとあまり違わないことになります。このことに気付いたので、νが幾らでも大きくなれる設定だとどうなるかを自分で試してみました。

タグ：数楽

posted at 21:11:20

#数楽添付画像はνが0.5～10^10の範囲を動きやすいような事前分布のもとでのMCMCで近似的に求めた事後分布の様子の一つです。

ρとνのスケールが対数になっていることに注意。νの値は見事に巨大になっています。事前分布を変えると劇的に推定結果が変わる。 pic.twitter.com/8K3JgaHt93

タグ：数楽

posted at 21:15:45

#数楽 t分布の自由度νの動く範囲を0.5から20程度までに制限して、正規分布に近付くことを許さない事前分布だと添付画像のようになります。回帰直線の傾きbの推定値が一つ前の添付画像では0.1を切りそうですが、こちらでは0.1を超えている感じ。 pic.twitter.com/yVaq1xm1ew

タグ：数楽

posted at 21:19:21

#数楽 1つ目はt分布が正規分布化することを許した場合の予測分布で、2つ目は正規分布化することを禁止した場合の予測分布です。点線は最小二乗法で引いた直線です。正規分布に近付くことを許すと最小二乗法の結果の方に近くなります。近付くことを許さないとうまく外れ値を処理するようになる。 pic.twitter.com/OhVvA4FrA3

タグ：数楽

posted at 21:23:32

#数楽さて、あなたは以上の2つの推定結果(および最小二乗法の結果)のどれを最も信用するでしょうか？以上のような「推定ごっこ」をすると常にこのような問題に悩まされます。

そこで、WAICを計算してみました。続く

watanabe-www.math.dis.titech.ac.jp/users/swatanab...

タグ：数楽

posted at 21:25:22

#数楽各々3000×8のchainsが得られたので、それらを全部合わせて、事後分布とみなして、WAIC(と最小二乗法のケースのAIC)を計算してみました(添付画像)。AICやWAICの値が小さいほど事後予測分布の予測精度が高い可能性が高い。続く pic.twitter.com/Oyzt9eDhXb

タグ：数楽

posted at 21:29:19

#数楽続き。この計算結果を見て、このケースでは、「事前分布をνが大きくなれないように設定すること」は十分正当化され得ると思いました。

AICやWAICの導出とその使い方は渡辺澄夫著『ベイズ統計の理論と方法』に詳しく書いてあります。

タグ：数楽

posted at 21:32:13

#数楽階層モデルでないモデルのケースでベイズ推定をMCMCで行ったのであればWAICは気軽にかつ簡単に計算できます。しかし、階層モデルの場合には重い数値積分が必要になるので、気軽に計算できる感じでは全然無くなるのだ。これは実用上相当に重要な問題だと思います。

タグ：数楽

posted at 21:34:11

#数楽実は、こういう場所で話題にするために画像を作るときには「見た目のきれいさ」をかなり優先しているところがあります。欲しい情報が正確にわかるようにしたい人はソースコードを全公開しているので、自分でやってみて下さい。

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

タグ：数楽

posted at 21:36:57

#数楽最尤法でt分布を使ったモデルを解く場合に、最適化函数に与える初期値を変えると真のマキシマムではないローカルマキシマムの方に収束してしまうことがよくあります。

大雑把には、最適化函数に与える初期値はベイズ推定で用いる事前分布に類似していると思いました。

タグ：数楽

posted at 21:42:51

#数楽続き。最適化函数に与える初期値を変えることによって収束先を変えることができるのと同じように、事前分布の取り方を変えることによって予測分布の収束先を大きく変えることができる。

タグ：数楽

posted at 21:44:21

#数楽続き。予測分布の予測精度の相対比較のための指標であるWAICは計算できて欲しい。「複数の予測分布のどれの予測を信頼した方が勝ち目が高いか」を知りたくなることは、現実世界への統計学の適用においては多いでしょう。

しかし、階層ベイズモデルだと全然簡単じゃない。

タグ：数楽

posted at 21:47:04

@tmiya_ ありがとうございます。見てみます。

タグ：

posted at 21:56:20

@todesking ありがとうございます。明日見てみます。

タグ：

posted at 21:56:46

#数楽訂正：「0から30くらい」→「1から30くらい」

twitter.com/genkuroki/stat...

タグ：数楽

posted at 22:17:10

トデス子'\ @todesking

@dekosuke 渡辺先生のクロスバリデーションの話ここで読めました watanabe-www.math.dis.titech.ac.jp/users/swatanab...

タグ：

posted at 22:18:58

#数楽個人的な意見では、外れ値の処理のためにt分布を使ったモデルを考えるときに、AICやWAICなどの情報量基準を計算せずに、「自由度νは大きくない」という仮定をモデル(事前分布の設定を含む)に組み込んでしまうのは*反則*だと思う。予測分布の精度を客観的な指標で比較して欲しい。

タグ：数楽

posted at 22:21:36

#JuliaLang Jupyterのおかげで「計算結果の保存」についてあんまり考えなくてもよいのですが、保存したい場合には

github.com/simonster/JLD2...

が便利。{at}save "all.jld2" で一挙に保存、all=load("all.jld2")

タグ： JuliaLang

posted at 22:25:54

き［ﾄﾁｭｳﾉｴｷﾆﾊ］よ［ﾄﾏﾘﾏｾﾝ @mexjp

#数楽「外れ値」のある場合にCauchy分布を使うという話があるようですが、私の実験

nbviewer.jupyter.org/gist/genkuroki...
(「比較」をページ内検索、AICを比較している)

ではCauchy分布＝自由度1のt分布を自由度νのt分布に拡張した方がうまく行きました。

タグ：数楽

posted at 22:53:24

作家？教授？矢吹樹氏の放送大学批判 - Togetterまとめ togetter.com/li/260513 @togetter_jpさんから　群大クビ教授がやらかしたのはこれらしいんだが、まさか5年越しでオチが付くとはな‥

タグ：

posted at 22:54:34

#数楽 ρ×(Cauchy分布) というモデルでρの値だけを推定するようにしてしまうと、最初から「裾野が激重で外れ値有り」で決め打ちしまくっていることになります。ρ×(自由度νのt分布)なら正規分布っぽい場合には推定で勝手にνが大きくなってくれます。続く

タグ：数楽

posted at 22:56:05

#数楽ただし、最尤法で ρ×(自由度νのt分布) モデルを解こうとすると、微妙な最適化過程での初期値依存性に苦しむことになります。その問題を緩和するためには「νは大きくならない」と決め打ちして最適化函数を使うといいです。その場合には正規分布のモデルとAICを比較した方がよい。

タグ：数楽

posted at 22:58:31

#数楽同様のことはベイズ推定する場合も同じで、「νは大きくならない」と決め打ちするための事前分布を採用しないと外れ値の処理がうまく行かないかもしれない。この場合にも(最低限の手続きとして)WAICを比較した方がよいです。

タグ：数楽

posted at 23:00:38

@yuzu_96wiz_only 【先生のそれまでの指導や児童の実態が不明なので、判断しかねます】
それでは、先生が「それは白だ」と言えば白になり、先生が「それは黒だ」と言えば黒になる、という事ですね。#掛算

タグ：掛算

posted at 23:02:15

#数楽ベイズ推定も最尤法もどちらも同じ尤度函数を扱っているので、どちらを扱っても結果的に尤度函数に関するなにがしかの知見が得られます。片方での経験はもう一方でも役に立ちます。「頻度主義 vs. ベイズ主義」とか言い出す人を次の世代の若い人達は完全無視するべきだと思う。黒歴史。

タグ：数楽

posted at 23:02:51

#数楽統計学の黒歴史に関わり合わずに数学を地道に勉強した方が健全。

・ベイズの定理はベイズ統計の基礎である←これは誤り

事実はこちら↓

・ベイズの定理は条件付き確率の定義から自明に出て来る結果に過ぎず、ベイズの定理を使ってもベイズ統計がどうしてうまく行くかを理解できない。

タグ：数楽

posted at 23:05:37

#数楽ここまで解説をつけてソースコードまで公開していれば、きっと誰かが、株価や為替や投資信託などの値動きの「外れ値」をCauchy分布よりも一般的なt分布を使って分析する人が出て来て結果を私に教えてくれるに違いない。😆

タグ：数楽

posted at 23:16:52

#数楽「その辺に落ちている数値列」はどれも結構面白い。

タグ：数楽

posted at 23:19:20

#数楽
 nbviewer.jupyter.org/gist/genkuroki... の最後の段落【Cauchy分布のモデルは極限として正規分布モデルを含んでいない～このような理由から、外れ値があっても、実際には正規分布である疑いがある場合にはCauchy分布モデル～を使用するべきではないと考えられる。】

タグ：数楽

posted at 23:21:55

@yuzu_96wiz_only 【先生が解釈の仕方の指導を記事から読み取れず、判断材料が不十分だからです】
結局、
「３匹の兎の耳の総数を求めよ」というとき
「２×３だけが正しい」となったり、
「２×３　と　３×２　のどちらでも良い」となったりするという話ですよね。　#掛算

タグ：掛算

posted at 23:38:18