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

#統計最尤法で最も気軽に使える情報量規準はAICとBICの2つなのですが、それらのベイズ推定法版のWAICとWBICを実装したつもりの #JuliaLang Jupyter notebook が次のリンク先にあります。

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

タグ： JuliaLang 統計

posted at 00:39:42

#統計すでにこの発言が繋がるスレッドの中で述べたように、確率密度函数は「測度」の「次元」を持つので、確率密度函数 p(x) のサンプル X_1,…,X_n の尤度　p(X_1)…p(X_n) は座標系 x に依存します。しかし、尤度比は座標系によらない量になる。続く

twitter.com/h_okumura/stat...

タグ：統計

posted at 01:25:14

#統計続き。尤度や対数尤度は座標系 x に依存するのですが、尤度比や対数尤度比は座標系 x に依存しない量になります。特に、尤度や対数尤度の大小関係は座標系に依存しない関係になります。ゆえに、最尤法は座標系に依存していないこともそのことからわかる。

タグ：統計

posted at 01:27:59

#統計すでに書いたことですが、座標系 x を座標系 y に x=y(x), y=x(y) によって座標変換すると、確率密度函数 p(x) のサンプルX_1,…,X_nに関する尤度は座標系 y にうつると |x'(y(X_1))|…|x'(y(X_n))| 倍されることになります。尤度比ではこれは分子分母でキャンセルする。

タグ：統計

posted at 01:31:04

#JuliaLang でもScikitLearnで簡単に遊べるんですね。実際に試してみたら本当に簡単でした。一部のコードは使用パッケージの側がまだ0.6に対応してなくて、0.5.2でしか動きませんでしたが。ほとんどの例はそのまま動いた。

github.com/cstjean/Scikit...

タグ： JuliaLang

posted at 01:34:12

Mara Averick @dataandme

ICYMI, 🌈 Incredible animations! "Sorting Algorithms Revisualized" buff.ly/2zWsW2X #visualization pic.twitter.com/UOp9gEB0qH

タグ： visualization

posted at 04:46:05

鰹節猫吉 @sunchanuiguru

分かるところから手をつけて、徐々に全体像が明らかになるという経験を積むことが大事と思われる。 twitter.com/tomatoha831/st...

タグ：

posted at 10:18:37

#統計もう一度まとめておくと、

* サンプルX_1,…,X_nに関する確率密度函数p(x)の尤度p(X_1)…p(X_n)は座標系xの取り方を変えると変わる。

* 2つの確率密度函数の尤度比は座標系xの取り方によらない。

* 特に尤度の大小関係は座標系xの取り方によらない。

タグ：統計

posted at 10:35:10

#統計単位系を変えることは、座標をスケール変換することに対応しているので、単位系の取り方への依存性を見ることは、座標系の取り方への依存性を見ることの特別な場合になっている。尤度の大小関係は座標系の取り方によらないので、特に単位系の取り方によらない。

タグ：統計

posted at 10:37:20

#統計真の分布q(x)に関する予測分布p(x)の汎化損失

G(q||p) = - ∫ q(x) log p(x) dx

は座標系xの取り方を変えると変わる。しかし、Kullback-Leibler情報量(＝相対情報量＝相対エントロピーの-1倍)

D(q||p) = ∫ q(x) log(q(x)/p(x)) dx

は座標系xの取り方を変えても不変である。

タグ：統計

posted at 10:41:04

#統計尤度比やKullback-Leibler情報量で処理できる場合には「次元を持つ量」の対数を考える必要はなくなる。

尤度比は測度論的には Radon-Nikodym derivative なので、尤度比を実用的に使ったことがある人は全員 Radon-Nikodym derivative ユーザーであったと言えると思います。

タグ：統計

posted at 10:44:46

#統計 KL情報量に関するSanovの定理で処理できる場合は論理的に色々クリアになっている点もうれしいです。「等確率の原理」とか言われると「どうしてそれでいいのか？」という難しい問題を考える必要が生じますが、「ベースになる確率分布は何でもよい。相対エントロピーを考えればよい」ならば簡単。

タグ：統計

posted at 10:47:52

#統計そして、所謂「指数型分布族」は「等確率の原理を仮定しているとは限らない場合のカノニカル分布」に一致しているので、応用対象を物理に限定したくなければ、「ベースになる確率分布q(x)は何でもよい。カノニカル分布p(x)はp(x)∝exp(-βH(x))q(x)の形になる」としておいた方が便利です。

タグ：統計

posted at 10:51:33

#統計例(ガンマ分布)　スケール1のガンマ分布の確率密度函数は

p(x|a) ∝ exp((a-1)log x) q(x), q(x) = exp(-x)

の形をしています。逆温度はβ=-(a-1)で、分配函数はガンマ函数です：

Z = ∫_0^∞ x^{a-1} e^{-x} dx = Γ(a).

ベースになる確率分布は指数分布 q(x) = exp(-x). 続く

タグ：統計

posted at 11:21:25

#統計以下、a>1と仮定します。このとき、逆温度β=-(a-1)は負になります。負の逆温度は普通によく出て来ます。

各々が指数分布 q(x)=exp(-x) に従う独立な確率変数列(X_1,…,X_n)の確率分布を相乗平均に関する不等式

log(X_1…X_n)^{1/n} ≥ c > -γ=-0.5772…

で制限すると～続く

タグ：統計

posted at 11:29:37

アップグレードカラーズ @UP_TKG

みなさんタイラップって買った時に袋の上を開封して、工具箱の中で全部出ちゃってたりいちいち開封口をテープで留めたり面倒な事してません？？
タイラップの袋の便利な開封方法はこうです。これなら勝手に袋から飛び出る事なく、必要なだけ取り出せて、工具箱のなかでも散らかりません！ pic.twitter.com/yWNiyhJD9j

タグ：

posted at 11:32:13

#統計続き～、nが大きなとき、実数直線上でのX_1,…,X_nの分布はガンマ分布でよく近似されるようになります。ただし、ψ(a)=c でパラメーターaを決める。ここで、

ψ(a) = (1/Γ(a))∫ e^{-x} x^{a-1} log x dx = (ガンマ分布での log x の平均)

です。ψ(1)=(指数分布でのlog xの平均)=-γ.

タグ：統計

posted at 11:32:54

#統計ガンマ分布はその特別な場合であるカイ二乗分布として、正規分布と同じくらいよく出て来る。大学新入生は、Gauss積分とガンマ函数をしっかり勉強しておくべきだと思います。しっかり勉強しておかないと、後で、統計学について勉強するときに困る。

タグ：統計

posted at 11:37:06

#統計学部生向けの「数学を避ける方針」で書かれている教科書では、ガウス積分やらガンマ函数やらベータ函数に関する説明も当然避けているので、大学1年生のときに勉強しておかないと、後で勉強する機会が無くなってしまう危険性さえあると思う。

タグ：統計

posted at 11:38:43

#統計実際に統計がらみの数値計算を色々やってみると、「コンピューターで計算すればいいから、ガウス積分やガンマ函数について知らなくても大丈夫だよね」とはならないと思いました。計算量を減らすために大学1年レベルの微積分の計算を使う機会は結構あるように思えた。

タグ：統計

posted at 11:40:59

#JuliaLang whichマクロが便利過ぎ。{at}which hoge(x) で実際に実行される函数 hoge(x) のソースコードのありかを教えてくれる。Julia言語ではhoge(x)の中身はxの型に依存して変わることが多いので、これ結構必須かも。

タグ： JuliaLang

posted at 11:49:22

#JuliaLang #統計

 nbviewer.jupyter.org/gist/genkuroki...
t-distribution linear regression by Mamba

にLOOCVの計算も追加した。WAIC, LOOCV, WBIC の簡単な計算例。計算法の解説付き。確率モデルp(x|w)、入力のサンプルX_k、事後分布のサンプルw_lに関する log p(X_k|w_l) から全部計算できる。

タグ： JuliaLang 統計

posted at 11:59:24

#JuliaLang #統計 MCMCをやってくれるパッケージには、確率モデルp(x|w)、入力のサンプルX_k、事後分布のサンプルw_lに関する

log p(X_k|w_l)

を自動的に計算してくれる機能が欲しいよな。これはあまりにも基本的な量なので、ユーザーに自前で計算させるようなものではないと思う。

タグ： JuliaLang 統計

posted at 12:04:58

Hiroyasu Kamo @kamo_hiroyasu

「邪道とは？」との疑問がありましたので、具体例をあげます。「はじき」とか「モルグリコ」とかです。

タグ：

posted at 14:15:49

Hiroyasu Kamo @kamo_hiroyasu

このあと、「はじき」や「モルグリコ」の擁護が出てきて、邪道を邪道と理解できない人が可視化されるでしょう。

タグ：

posted at 14:19:20

#JuliaLang #統計最尤法の計算のコードが一ヶ所バグっていたので直した。ついでにGadflyじゃなくて、自前でPyPlotを使ってプロットするように書き換えた。

nbviewer.jupyter.org/gist/genkuroki...
t-distribution linear regression by Mamba

Julia言語でNUTSでMCMC。WAIC、LOOCV、WBICの計算。

タグ： JuliaLang 統計

posted at 15:48:31

#JuliaLang まだ

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

の方には反映されていないが、

gist.github.com/genkuroki/1c9f...

の方はもう置き換わっている。

Julia言語でNUTSでMCMC。WAIC、LOOCV、WBICの計算。MCMCの結果を初期条件として使って最尤法。AICとBICの計算。

タグ： JuliaLang

posted at 15:51:58

ぷらぎあ @plastic_gear

中華3Dプリンタ、思ってたよりはしっかりしてる
必要工具やら日本用電源全部ついてくるし動画付きの説明書だし、この梱包の綺麗さ
こんなのが二万円代とか勝ち目ないな…と思って組み立て中悲しくなった pic.twitter.com/Od2NrSZB8D

タグ：

posted at 17:59:30

とにぃAntonio Camillo @Antonio_Camillo

#JuliaLang 自由エネルギーを数値積分で直接計算する函数を追加した。WBICは逆温度1/log(n)での被∫函数の値で積分を近似したものである。これで正しいかについてまだ確信はない。

gist.github.com/genkuroki/1c9f...

nbviewer.jupyter.org/gist/genkuroki... pic.twitter.com/LLtBtXe7wF

タグ： JuliaLang

posted at 20:17:52

Permalink - 2017年11月12日

2017年11月13日(月)

すごいもんですねぇ
匠の技ですな @Limportant_fr

pic.twitter.com/KtCdkGjeoX

タグ：

posted at 01:50:46

Matthew Wingate @WingateMatthewB

Just read this at link below. Funny cuz it's true [tho I'm still learning]. #JuliaLang ucidatascienceinitiative.github.io/IntroToJulia/H... pic.twitter.com/dZOMVYSY9g

タグ： JuliaLang

posted at 03:44:14

ネタ再掲 #統計
 www.jstage.jst.go.jp/article/butsur...
エントロピーとモデルの尤度
赤池弘次

で赤池さんは

AIC (an information criterion の略記)

と書いていて、それに編集委員会が脚注を付けている(笑)。 pic.twitter.com/e6v49w9jlC

ismrepo.ism.ac.jp/?action=pages_...
もおすすめ。

タグ：統計

posted at 10:46:59

Hiroki MORI (森裕紀) モ @HirokiMori

#統計赤池弘次さんの論説

ismrepo.ism.ac.jp/?action=pages_...
統計的推論のパラダイムの変遷について

www.jstage.jst.go.jp/article/butsur...
エントロピーとモデルの尤度

「尤度ってどうしてもっともらしさなの？」「客観確率 vs. 主観確率という対立図式はおかしいよね？」のような疑問を持つ人達には特におすすめ。

タグ：統計

posted at 10:48:27

@HirokiSayama 佐山先生、この前言っていた話です。 twitter.com/genkuroki/stat...

タグ：

posted at 10:50:45

#統計以上は
twitter.com/genkuroki/stat...
の一部分の再掲。

「尤度をもっともらしさだ」と言いながらその理由を説明しないとか、ベイズ統計について「客観確率 vs. 主観確率」という対立図式に基いて説明したりする人達は「反赤池の立場」に立っていることになる。それって相当に怖いことだと思う。

タグ：統計

posted at 10:52:07

Jun Окаdа @ojunn

「AIC (an information criterion の略記)」笑 twitter.com/genkuroki/stat...

タグ：

posted at 10:56:43

#統計 WBICの原論文の定理3の証明中にある定積分で自由エネルギーを表示する公式で自由エネルギーを計算しても、1変数函数の数値積分なので待ち時間はほとんどありません。添付画像は #JuliaLang での実装例。

nbviewer.jupyter.org/gist/genkuroki... pic.twitter.com/oHhoDkTMni

タグ： JuliaLang 統計

posted at 12:12:21

#統計 #JuliaLang 被積分函数は E2nLn(β)=2E^β_w[nL_n] (2倍してあるのはBICの伝統的なスケールに合わせるため)で、quadgkが数値積分の函数。被積分函数の計算には逆温度βでの事後分布は必要なくて、逆温度1の通常の事後分布だけで十分です。

nbviewer.jupyter.org/gist/genkuroki...
にも解説があります。

タグ： JuliaLang 統計

posted at 12:15:44

#統計理論的にも言えていることなのですが、モデルが正則ではなく、特異な(場合に近い)ときには、サンプル数が小さくても収束し**易く**なります。特異モデルでは収束が**速く**なる。事後分布が広がったままで、事後予測分布が収束してしまう。

タグ：統計

posted at 12:21:16

やねうら王 @yaneuraou

弱かろうと将棋ソフトをフルスクラッチで書いている学生さんとかが本当のガチ勢であって、ライブラリ使ってお手軽に勝つことを目指しているのはエンジョイ勢だと私は思う。前者の人達が大会で「二週間で出来ますた！」みたいなライブラリ勢のソフトにボコボコにされて「自分エンジョイ勢なんで…」と…

タグ：

posted at 13:17:24

やねうら王 @yaneuraou

エクスキューズしてる姿に私は心を痛めている。「君たちこそが、本当のガチ勢なんだよ。もっと胸を張っていいんだよ。」と言ってやりたいが、やねうら王ライブラリ勢の親玉だと思われている私が言っても「何言ってんだ、こいつ」みたくなるので言わないし、言えない。

タグ：

posted at 13:23:45

@kuri_kurita

斉藤ひでみ・現職教師（西村祐二） @kimamanigo0815

「復興増税」を言い出したのは13日だったそうだし、絶対にゴジラが暴れてる最中に増税を言い出したはず。

タグ：

posted at 14:04:19

非公開

タグ：

posted at xx:xx:xx

【教員記者会見の映像公開】

11/6に行ないました、現職審議会による記者会見を、ほぼノーカットで公開しました。
現場の教員による訴えを、どうか聴いて下さい。
僕たちは引き続き、中教審・文科省等に働きかけを続けます。

www.kyodo-bukatsu.net/gs/kaiken/
※前半もですが、後半の質疑応答も見所です pic.twitter.com/JQqWRKLRw9

タグ：

posted at 17:29:23

三村智保 @igomimu

内田良：新刊『だれが校則を決めるのか』　 @RyoUchida_RIRIS

超期待！ twitter.com/asahi_igo/stat...

タグ：

posted at 17:37:07

どゎゎゎ・・・
伝説の記者会見，動画公開に踏み切った・・・

国の審議会が始まる直前に，隣の会議室で，現職の声をアピール。
しかも，迅速に動画を編集，審議会の議事録をアップロードと，すごすぎる。 twitter.com/kimamanigo0815...

タグ：

posted at 18:08:26

#統計 loglik[l,i]=log p(X_i|w_l)です(X_iは真の分布が生成したサンプル、w_lは逆温度1の普通の事後分布のサンプル)。WBICや0から1までの定積分による自由エネルギーの計算には、逆温度βの事後分布のサンプルをMCMCで作る必要はないです。

nbviewer.jupyter.org/gist/genkuroki...
pic.twitter.com/oHhoDkTMni

タグ：統計

posted at 19:38:06

#統計 #JuliaLang pic.twitter.com/oHhoDkTMni
loglik[l,i]=log p(X_i|w_l)でsize(loglik)[1]=size(loglik,1)はlの個数(右辺のように書くべきだったか)。

sum(loglik,2)＝Σ_i loglik[:,i]

E2nLn(β)=E^β_w[nL_n]は逆温度1の事後分布のサンプルを使って逆温度βでの平均を計算する函数。

タグ： JuliaLang 統計

posted at 19:42:18

#統計 #JuliaLang {at}sum マクロは for ループで和を取るためのマクロです。Julia言語では for ループで「ダサく」和を取った方がメモリが節約されてかつ計算も速いことが結構あります。E2nLn(β)はβで数値積分するので繰り返し何度も呼ばれるのでforループで和を取って計算することにしました。

タグ： JuliaLang 統計

posted at 19:45:22

#統計 #JuliaLang 以上のような方法である場合に計算したFreeEnergyはWBICより少し小さな値になってかつ、MCMCのサンプルの取り方に関する分散も小さいです。私はWBICを使うよりももっと素朴に数値積分した方がメリットがあるのではないかと思っているのですが、まだ確証はありません。

タグ： JuliaLang 統計

posted at 19:47:26

歩行者は右側通行 @JikanBae

自己相関関数やらケプストラムやらの図を作るのにJulia（とPyPlot）を使っていたら、それらの関数がどんどん蓄積されてゆく。車輪の再発明だしオレオレ分析関数ばかりになるけど、使いやすくなっていくなぁ。 #JuliaLang

タグ： JuliaLang

posted at 20:55:01

#JuliaLang

using KernelDensity
function makefunc_pdfkde(X)
local ik = InterpKDE(kde(X))
local pdfkde(x) = pdf(ik, x)
return pdfkde
end

1次元のサンプルXからKDEで推定した確率密度函数を返す函数

使い方：
kde_sample = makefunc_pdfkde(sample)
kde_sample(1.5)

タグ： JuliaLang

posted at 21:29:22

#JuliaLang これを以下のように「短く」しては**いけない！**

using KernelDensity
function makefunc_pdfkde(X)
local pdfkde(x) = pdf(InterpKDE(kde(X)), x)
return pdfkde
end

こうやって作った確率密度函数は激遅になる。ik = InterpKDE(kde(X))を1回計算して使い回さないとダメ。

タグ： JuliaLang

posted at 21:31:29

#JuliaLang KDEは離散的なサンプルから連続的なプロットを行うときによく使われる。サンプル(配列)が二次元(2個)の場合には

function makefunc_pdfkde(X,Y)
local ik = InterpKDE(kde((X,Y)))
local pdfkde(x, y) = pdf(ik, x, y)
return pdfkde
end

密度による色付けによく使っている。

タグ： JuliaLang

posted at 21:34:08

おお！英数字をたくさん書けるようになったことは、ミニコードを貼り付けるときに便利！

タグ：

posted at 21:39:47

#JuliaLang 以上の準備のもとで

import PyPlot; plt=PyPlot
n=10^3
X,Y=randn(n),2randn(n)
X,Y=X+Y,X-Y
f=makefunc_pdfkde(X,Y)
plt.scatter(X,Y,c=f.(X,Y),s=5,cmap="CMRmap")
plt.axes()[:set_facecolor]("k")
plt.axes()[:set_aspect]("equal")

using PyPlotすればplt.は不要になる。 pic.twitter.com/vks0QcUfZ7

タグ： JuliaLang

posted at 21:55:08

#JuliaLang Plotsパッケージじゃなくて、PyPlotパッケージを使うようになってしまいました。Python情報が大量に出回っているのでそれを「真似」すればいいので楽。ただし、axes()[:set_aspect]("equal")のような使い方ができないとPython→Juliaの翻訳は困難。

タグ： JuliaLang

posted at 21:59:47

#JuliaLang Julia経由でPythonのmatplotlib.pyplotを使うときに、日本語フォントを使う方法のまとめが次のリンク先にあります。

gist.github.com/genkuroki/e748...
日本語フォントをPyPlotで使用する方法

この手の「数学」と無関係なところで時間がかかってしまう。

タグ： JuliaLang

posted at 22:04:17

#JuliaLang の確率分布パッケージDistributions.jlの考え方は色々素直でいいです。例えば、ガンマ分布に従う乱数の100個生成は

d=Gamma(25,0.04)
X=rand(d,100)

だし、確率密度函数や累積確率函数は

pdf(d,x)
cdf(d,x)

です。確率分布dを各種確率がらみの函数に代入すると欲しい結果が得られる。

タグ： JuliaLang

posted at 22:09:43

#JuliaLang 混合正規分布のような混合分布もDistributionsパッケージで定義できます。後、結構最近のことなのですが、1次元の確率分布を平行移動とスケール変換してできる確率分布を定義する仕組みも整備されました。

タグ： JuliaLang

posted at 22:11:18

#JuliaLang 自分で確率分布を定義することも簡単です。自分で定義した確率分布dの乱数の発生はquantile()を実装しておけばquantile(d,rand())で勝手にやってくれるようです。もちろんより高速なrandを直接実装してもよい。ソースコードを見れば大抵のことがわかります。

タグ： JuliaLang

posted at 22:13:51

#JuliaLang を使うときに、最初に覚えた方がよいことは「重い計算はすべて函数の中に入れること」です。函数の外にナマで書くと非常に遅くなります。

あと、zero(x)とかone(x)のような書き方。xと同じ型の0や1です。xと同じ型の最大値もしくは無限大はtypemax(x)です。続く

タグ： JuliaLang

posted at 22:20:12

#JuliaLang では函数に何かを代入して計算させようとすると、代入したものの型に応じて、可能な限り適切にすぐさまコンパイルして実行してくれます。様々な型のものを代入しても高速に計算してくれる可能性を高くするためにはone(x)とかzero(x)のような書き方をしておくといいかも。

タグ： JuliaLang

posted at 22:21:49

#JuliaLang 私が非常に感心したのは、そのようなスタイルで書かれた連分数展開を利用させた数値計算用の函数に「文字変数」を代入してLaTeX的に綺麗に成形された連分数の式を出力させることができることを知ったときです。数値計算のコードはまず「紙の上で」最適化するものだと思っていました。続く

タグ： JuliaLang

posted at 22:23:37

#JuliaLang 続き。考え方が古過ぎだった。数値計算用の最適化されたコードの生成は高級言語の各種マクロを使って自動化した方がよい。以前にも紹介しましたが、そういう講義がMITではあるようです。以下のリンク先にJupyter notebookがあります。

nbviewer.jupyter.org/github/steveng...

タグ： JuliaLang

posted at 22:26:33

#JuliaLang の魅力は、言語そのものの魅力だけではなく、その周辺で様々な解説をしてくれる人達の解説が非常に面白いところにもあると思います。Julia言語使いの人の解説は面白いことが多いです。

例えば、次のリンク先には単純な多項式計算の高速化を扱っています。
nemocas.org/benchmarks.html

タグ： JuliaLang

posted at 22:30:14

#JuliaLang nemocas.org/benchmarks.html を見ると、Julia言語で書かれたNemo.jlの多項式計算が全般的に速そうなことがわかります。巨大多項式の高速計算がもしも必要になったら、試してみる価値があると思います。高速多項式計算で困難に出会ったら、解決してくれそうな雰囲気もあるように思える。

タグ： JuliaLang

posted at 22:32:24