黒木玄 Gen Kuroki
- いいね数 389,756/311,170
- フォロー 995 フォロワー 14,556 ツイート 293,980
- 現在地 (^-^)/
- Web https://genkuroki.github.io/documents/
- 自己紹介 私については https://twilog.org/genkuroki と https://genkuroki.github.io と https://github.com/genkuroki と https://github.com/genkuroki/public を見て下さい。
2018年01月29日(月)
#Julia言語 最大公約数が1になる確率を利用したπの計算を、最大公約数を求めずにやってみました。
（同じ計算を3回やったらなぜか速くなった） pic.twitter.com/7StS4vzMx1
タグ： Julia言語
posted at 00:07:15
function count_relatively_prime_pi(n)
φ = fill(1, (n))
for p in 2:n
φ[p] != 1 && continue
for i in p:p:n
φ[i] *= p - 1
end
pの巾乗 = p^2
while pの巾乗 <= n
for i in pの巾乗:pの巾乗:n
φ[i] *= p
end
pの巾乗 *= p
end
end
√(6n^2 / (2 * foldl(+, φ) - 1))
end
タグ：
posted at 00:07:16
ちなみにJScheme版。
こっちを先にやったら予想外に速かったので Julia でもやってみたのでした。
n=20000 のとき Scheme Droid で1秒以上かかります。
1 と 1 の最大公約数を2回数えてる計算になるので結果が少し不正確(？)です。 pic.twitter.com/kADjzuW1Z6
タグ：
posted at 00:07:20
昨日数学カフェの予習会で話して、その準備で久しぶりにテンソル計算をガリガリやったら、大学出て7年くらい経つけれど意外と計算力落ちてなかった。サーベイの方法とかプレゼン方法とかを会社で学んで、数学やるための総合的な技術はむしろ上がってるのかもしれないとふと思った（だといいな...）
タグ：
posted at 07:28:58
#Julia言語 空の配列 [] を用意してそれにpushしまくる書き方ってもしかして普通？わたしゃ、頭が古いせいか、配列は必要なだけ確保したいと感じるのだが。pushのオーバーヘッドが気になる。
L = 10^6
a = []
for i in 1:L
push!(a, rand())
end
は非常に遅い。
nbviewer.jupyter.org/gist/genkuroki...
タグ： Julia言語
posted at 07:50:55
念のため補足すると，大学でサーベイやプレゼン技術が身につけられないと言いたいのではなく（指導教員の先生にはお世話になりました）今から思うと大学生の自分はこれらの能力に改善の余地があったという意味です．逆に企業で数学の能力が必要ないかと言えば，そんな事はなく大いに必要だと思います
タグ：
posted at 08:00:26
#Julia言語 本気で配列を使う計算を効率化したければ、全てをin-place計算で行う根性が必要。コードが見難くなる。この点はもっとどうにかならないかと思います。しかし、効率化の手段があるのはとてもありがたい。
nbviewer.jupyter.org/gist/genkuroki...
タグ： Julia言語
posted at 08:02:52
Julia、三次元プロットをPyPlot使わずにやるのが難しい。未だに
x y z
と横に並んだデータをそのままプロットする方法がわからない。単純な格子に並んでなくてもgnuplotみたいにメッシュ表示してほしい
タグ：
posted at 08:13:05
#仮面ライダービルド
昨日の量子誤り訂正のツイートを沢山の方に見ていただいたみたいですが、そういえば番組開始前に下のツイートをしてました。
ちなみにこのプロジェクトのHPは↓トーリックコードの実装により量子状態を長時間守り続けることはメインの挑戦的課題の一つ。
www.jst.go.jp/erato/nakamura... twitter.com/makoto0218ne56...
タグ： 仮面ライダービルド
posted at 08:36:09
twitter.com/cometscome_phy...
#Julia言語 PyPlotなら
github.com/gizmaa/Julia_E...
Plotsなら、それのPyPlotをPlotsに、surfの行をsurface(x,y,z)に置き換えれば行けます。
タグ： Julia言語
posted at 09:21:27
思考実験：「1000本に1本の割合で当たりが出る」（仮説A）とされているクジを40回引いたら1本当たった。Aが正しいなら、40回中の当たりの本数Xは二項分布B(40,0.001)に従い、P(X=0)=0.999^40≒0.961なので、確率0.04以下の現象が起こったことになる。以上の考察に問題はあるでしょうか？（続く）
タグ：
posted at 09:29:37
仮に問題ないとして、これを仮説検定に当てはめてみると、帰無仮説p=0.001（仮説A）が有意水準5％で棄却されることになる。よってAは間違っており、当たりくじの比率は0.001より高いと考えられる。この判断は正しい？正しくないとしたら、その根拠は？（続く）
タグ：
posted at 09:30:58
考察1
『上の検定では対立仮説を「p>0.001」としていることになり片側検定。p<0.001の可能性を排除する根拠が示されていない以上、両側検定の考え方に基づいてp値を2倍すべきなので、有意水準5％では棄却されない。よってAを否定するのは間違い』（続く）
タグ：
posted at 09:32:33
考察2
『Xは非負整数なので、両側検定における左側の棄却域は設定不可能。よって（むしろ両側検定は不自然で）片側検定で判断して良い。よってAを否定する判断は正しい』（続く）
タグ：
posted at 09:33:51
さてどっちが正しい？
ちなみに『0.001が95%信頼区間に入るので仮説Aは否定できない』といった見方もあると思いますが、区間推定と両側検定の食い違いは概ね計算方法の違いによるもので、考察1の主張と本質的な違いは無いだろうと思います（この点についてもご意見ありましたらどうぞ）。（続く）
タグ：
posted at 09:35:13
その上で
考察3
『そもそも有意水準や信頼係数で機械的に判断するのが間違い。両方の考え方を踏まえた上で、状況に応じて「確率が低過ぎるかどうか」「高い確率であり得るかどうか」を検討して判断すべき』
タグ：
posted at 09:36:58
Masayuki Ohzeki＠雑談方程 @mohzeki222
手でトーラス符号を解析できる稀有な理論です。小数点何桁までっていうマニアックな戦いをしていたのは良い思い出だ。またここら辺の話が気になってきた。
タグ：
posted at 10:12:30
#Julia言語 gcd計算の効率を測る主旨の意図的に効率の悪いコードだったのですが、こんなに速くなるとは！
Julia言語での函数の1回目の実行時間はコンパイル時間を含みます。2回目以降はコンパイル結果を使います。だから @ time で測った時間が短くなる。
twitter.com/brv00/status/9...
タグ： Julia言語
posted at 10:36:43
中澤 港％人類生態学者＠神戸大学 @MinatoNakazawa
@yoriyuki 「網羅」というと微妙ですが（とくに実務は分野によって違うので），CRANのDocumentationのContributed
cran.r-project.org/other-docs.html
が，かなり充実しています。
cran.r-project.org/doc/contrib/Sh...
cran.r-project.org/doc/contrib/us...
など素晴らしいです。
タグ：
posted at 10:59:01
やり方がわからない人は卒論ファイルを自分のGmailアドレスに送りつけろ．情報系なら全部BitBucketにつっこめ．今すぐ！ twitter.com/hyuki/status/9...
タグ：
posted at 11:02:26
非公開
タグ：
posted at xx:xx:xx
尤度函数の様子を少数の数値で代表させずに見せて欲しいと感じることが実に多い。
統計分析の再検証を可能にするデータかダウンロード可能になっていればさらによい。
あと、予測分布の様子も見せて欲しい。
タグ：
posted at 12:47:56
データを見せて欲しいと感じるのは、尤度函数(ベイズ推定の場合は事後分布)や予測分布の様子を見せてくれずに、少数の数値しか見せてくれないから。
少数の数値で代表する習慣も、科学的な要請ではなく、歴史的社会的経済的な理由で定着しているのだと思う。
タグ：
posted at 12:52:28
ベイズ推定の統計パッケージには事後分布の様子をプロットする機能がついているので、事後分布の様子を見せるのは簡単。最尤法(最小二乗法を含む)でも尤度函数を見せて欲しいと思う。
タグ：
posted at 13:14:09
ベイズ統計の解説を選ぶときには「予測分布」の正しい定義を書いてあるものを選ぶべきだと思う。
そして、モデルや予測分布の精度を比較するための指標に触れていない教科書を読んでも、
単にフィッティングしただけ
になってしまい、科学的な信頼性は失われる可能性が高い。
タグ：
posted at 13:22:52
実質的に事後分布を通して「尤度函数の様子を見ているだけ」になっているベイズ統計の入門書が多いのではないか？
モデル全体=確率モデル+事前分布
の選択法にふれないものだから、事後分布も選択される対象にできなくなっている。
尤度函数を特定の事前分布経由で見ているだけになってしまう。
タグ：
posted at 13:29:03
あれ、Juliaで[1;2]って列ベクトルにならんのか。2×1 Array{Float64,2}:じゃなくて、2-element Array{Int64,1}:になる
タグ：
posted at 13:45:50
ベイズ統計の予測分布の定義は
p^*(x) = ∫ p(x|w)φ(w|D)dw
です。ここで、p(x|w)はパラメーターwを持つ確率モデルで、φ(w|D) はサンプルDに対する事後分布です。MCMCで事後分布のサンプルw_kを作っていた場合には
p^*(x) ≒ (1/L)Σ_{k=1}^L p(x|w_k)
で計算できる。続く
タグ：
posted at 13:48:14
最尤法の予測分布の定義は
p^*(x) = p(x|w^*)
です。ここで w^* は尤度函数を最大にするパラメーター。
以上の式を比較すればベイズ推定と最尤法の予測分布の形はそっくりなことがわかります。最尤法はベイズ推定法に逆温度βというパラメーターを入れたときのβ→∞(絶対零度)の極限になっています。
タグ：
posted at 13:50:12
モデルが簡単な場合(単純な正則モデルの場合)には最尤法とベイズ推定はほぼ同じような結果を与えます。最尤法の方が計算量が小さくて優れている場合が少なくないものと思われます。そのようなケースにまでベイズ推定法を使う必然性はない。必然性がない場合に力を入れたベイズ推定法の解説は苦しい。
タグ：
posted at 13:52:09
最尤法によるパラメーターの推定が不安定になるケースであっても、ベイズ推定法の予測分布が安定してより信頼できる結果を与えることは数学的に証明されています(渡辺澄夫さんによる)。ベイズ推定の正しい予測分布の定義を採用しないと、そういうベイズ推定法のメリットを享受できません。
タグ：
posted at 13:55:54
あと、予測分布と実際に得られた分布の比較をきちんとすると、事後分布を見ていただけではわからないようなこと(例えば全然うまくフィットできていないとか)がわかります。予測分布はベイズ推定法に限らず、最尤法でも重要だと思います。
タグ：
posted at 13:57:26
「パラメーターの推定」という発想は人為的で、歴史的社会的経済的に権威付けられた発想だと思う。それに対して「確率分布の推定」は普遍的でクリアでわかりやすい。
タグ：
posted at 13:59:50
パラメーターの推定の結果は最尤法では尤度を最大にするパラメーターであり、ベイズ推定では事後分布になります。
それに対して、確率分布の推定結果は、最尤法であってもベイズ推定法であっても、どちらも予測分布になります。予測分布どうしは比較可能です。
タグ：
posted at 14:00:49
サンプル(データ)があるときに、それを再現すると期待する確率分布が2つあるとき、どちらの確率分布がよりよくデータを再現するかの指標として、尤度が使えます。実際の計算では対数尤度を使います。確率分布p(x)のサンプルX_1,…,X_nに関する対数尤度の定義は
Σ_{i=1}^n log p(X_i)
です。続く
タグ：
posted at 14:03:07
これの値が大きい方の確率分布の方が相対的に信頼できる。その理由は大数の法則とSanovの定理を使った簡単な議論でわかります。その辺のアイデアの開拓者である赤池弘次さん自身の解説をインターネット上で読めます。
twitter.com/genkuroki/stat...
タグ：
posted at 14:06:46
リンク切れしていたので再リンク
次の２つは必読。
www.jstage.jst.go.jp/article/butsur...
エントロピーとモデルの尤度
赤池弘次
ismrepo.ism.ac.jp/index.php?acti...
統計的推論のパラダイムの変遷について
赤池弘次
前者には
この統計量は AIC (an information criterion) と呼ばれ
とある😅。ダウロードしておくべき。
タグ：
posted at 14:13:16
どちらも1980年出版なのですが、21世紀の現代になっても、それらの赤池さんによる解説を超える解説は出ていないのではないか？
大学の授業で学生が「頻度主義 vs. ベイズ主義」のようなおかしなことを言い出す先生に出会ったら、赤池さんの解説を読んで猛爆撃をしかけるべき。
タグ：
posted at 14:16:19
赤池さんの論説であれば、権威の側面に負けることはありえない。
権威があってかつ正しいことを言っている人が我々の大先輩として存在しているのだから、どんどん利用するべきだと思う。
タグ：
posted at 14:22:58
で、予測分布の選択を対数尤度を見るだけで行うことが不適切な場合にどうするかの一つの答えが赤池情報量規準AIC (赤池さん自身が an information criterion の略だとしている文献はすでに示した(笑)) なわけです。そしてベイズ推定への適切な拡張はWAICです。
タグ：
posted at 14:25:52
予測分布の精度の指標はAICやWAIC。
それとは違って、もとのモデル全体＝確率モデル＋事前分布がサンプルを生成した確率分布にどれだけ近いかを測る指標がBICやベイズ自由エネルギーとその近似値達。
タグ：
posted at 14:27:50
正規分布＋共役事前分布の場合にはすべてをexactに手計算できます。正確な公式を使って数値計算できるようにした #Julia言語 カーネルの Jupyter notebook が
nbviewer.jupyter.org/gist/genkuroki...
正規分布の共役事前分布
にあります。これ、初学者にとって極めて有用な基本的な話が書いてあるノートのはず。
タグ： Julia言語
posted at 14:30:27
正規分布＋共役事前分布の場合に任意の逆温度βについて、分配函数、ベイズ自由エネルギー、事後分布、予測分布、WAIC、LOOCV、WBICなどの閉じた公式が書いてあり、しかも #Julia言語 の函数として定義されています。
nbviewer.jupyter.org/gist/genkuroki...
正規分布の共役事前分布
タグ： Julia言語
posted at 14:32:53
ベイズ推定について、正規分布の簡単な場合にさえ理解していないのに、もっと複雑な場合についてMCMCを回してみたりしていたのですが、なまいきでした。最初に正規分布＋共役事前分布の場合を試してみるべきだった。全部、手計算で閉じた公式を出せる。大学1年レベルの地道な計算をするだけ。
タグ：
posted at 14:34:32
「ベイズ推定の統計パッケージを使えば事後分布の様子を見せるのは簡単」という話に戻る。それでは、予測分布を見せるのは簡単か？ これは結構面倒です。特に「階層ベイズモデル」の場合、すなわち、確率モデルが内部パラメーターzによる積分
p(x|w)=∫p_1(x|z)p_2(z|w)dz
となっているときが面倒。
タグ：
posted at 14:43:02
p(x|w)=∫p_1(x|z)p_2(z|w)dz に関する明示的な公式がないので、p_1(x|z) と p_2(z|w) を使った階層モデルにして、モデルを数値的に解いているで、p(x|w) を扱うためには z に関する数値積分が別に必要。
階層ベイズモデルを理解している人は必然的に数値積分にも興味を持つことになります。
タグ：
posted at 14:45:14
確率モデル p(x|w) を p(x|w)=∫p_1(x|z)p_2(z|w)dz で階層化して数値的に解いている場合には、予測分布の計算にその積分の数値計算が必要になり、芋づる式にWAICやLOOCVやWBICなどの計算にも数値積分が必要になってしまう。
この問題は個人的に重要だと思う。
タグ：
posted at 14:47:06
Juliaだとf2(x) = x^2がラムダ式（と言うかまだ知らないけど）なのね。Pythonだとf2 = lambda x: x**2と同じと考えていいのかな～
タグ：
posted at 14:48:35
非公開
タグ：
posted at xx:xx:xx
これ本当に何処かが音頭をとって何とかして欲しさはある。言語の生まれ育ったその環境，歴史，思想すべてぶち込んで表すことができればいい(突然の向井秀徳)んだけど絶対まとまらない。
タグ：
posted at 16:31:08
非公開
タグ：
posted at xx:xx:xx
@antimon2 @genkuroki
関数呼び出しもコンパイルされるんですね。（そりゃそうか）
メモリが関係ないことは、1回目に渡す値だけを2に変えてみても結果が変わらなかったので、気づいていたのですが…。
twitter.com/antimon2/statu...
twitter.com/genkuroki/stat...
タグ：
posted at 19:02:28
Anaconda Navigatorがすこぶる使いやすくなってる…。
いつから仮想環境の管理ができるようになったのか。
チャンネル追加してパッケージ管理までGUIでできる。
Julia専用Python環境を作りたいが可能だろうか。
タグ：
posted at 19:07:57
@brv00 @genkuroki #Julia言語 のJITコンパイルは関数単位です。定義されてから最初に呼び出される時にコンパイルが走ります。2回目以降は〜は @genkuroki さんの解説の通り
タグ： Julia言語
posted at 19:20:31
@antimon2 @genkuroki あ、定義時にはコンパイルされないで、最初の呼び出し時にコンパイルされるんですね。全然わかってませんでした。ありがとうございます。
タグ：
posted at 19:26:37
ごまふあざらし(GomahuAzaras @MathSorcerer
#julia言語 またJuliaにすぴーどめまけた (´・ω・｀) vs Python with Numba pic.twitter.com/BtmIPhVWzu
タグ： julia言語
posted at 19:56:33
ごまふあざらし(GomahuAzaras @MathSorcerer
Python のほうは while pow_p <= n: じゃないといけない気がしてきた．
タグ：
posted at 20:11:11
ごまふあざらし(GomahuAzaras @MathSorcerer
なんでこのプログラムでええんだろうかって考えた結果Euler関数（ja.wikipedia.org/wiki/%E3%82%AA...）の計算をしているということだった・・・．
タグ：
posted at 20:32:46
#Julia言語 自分で使う範囲内だけでもPyPlotからPlotsへの翻訳をできるようになりたいというモチベーションで次のノートブックを作成しました。
nbviewer.jupyter.org/gist/genkuroki...
PyPlot と Plots の比較
1つ目はPyPlot、2つ目はPlots pyplot backend pic.twitter.com/g4jjpKT2gt
タグ： Julia言語
posted at 20:32:51
#Julia言語
添付画像:
1つ目はPyPlot
2つ目はPlots pyplot() backend
3つ目はPlots gr() backend pic.twitter.com/EAbo5HoODg
タグ： Julia言語
posted at 20:34:43
#Julia言語 小ネタ：Julia言語では
A = [
11 12 13
21 22 23
]
のとき、A[2,3] の値は 23 になるのだが、
A[1] == 11
A[2] == 21
A[3] == 12
A[4] == 22
A[5] == 13
A[6] == 23
となる。 pic.twitter.com/3x56rdrjVO
タグ： Julia言語
posted at 21:07:38
私は授業で「列ベクトル」「行ベクトル」と言うときには必ず口頭で「縦ベクトル」「横ベクトル」と言い直してかつ、板書上の「列」の上に「たて」と、「行」の上に「よこ」とルビをふっている。
タグ：
posted at 21:09:35
【第 6 回統計物理学懇談会のお知らせ】
プログラムが確定しています。また、いくつかの講演のタイトルも掲載しました。
それにしても初日最後の渡辺澄夫さんの（←同期の親友なんで普段は渡辺って呼んでるけど）タイトルにはもうめちゃめちゃ期待が高まることであるよ！
www.gakushuin.ac.jp/~881791/spm/20...
タグ：
posted at 22:07:13
julia proですが、アンインストールして（自動ではちゃんとできなかった）再インストールすれば簡単にインストールできました。
セキュリティソフトが変なことをするかもなので、セキュリティソフトは止めておいたほうがいいと思います。
タグ：
posted at 23:42:36
最近Juliaで計算機科学するのが流行ってるけど、zeusやplutoをcからJuliaに移植するとかしたら需要あるかなぁ…
条件設定はスクリプト化されてるけど、やっぱコードいじる部分があったし、今風のコーディングできたら研究も楽だったなと。
タグ：
posted at 23:42:55
