Mathtodon
#gist

たまに誰かに自分が書いたものを送り付けると有用なバックアップになったりします。

GitHub Gistにシークレットで(URLを知らないとアクセスできない形式で)書いたものを投稿しまくるという手もあると思う。

最近のエディターだとワンタッチでGistに編集中のファイルをアップロードしてくれたりします。

gist.github.com/genkuroki/509a
カイ二乗分布のGIFアニメーション

これを見れば

・Julia言語で確率分布とそれに従う乱数や確率密度函数を扱う方法

・GIFアニメーションの作り方

の両方がわかります。

余計な装飾抜きにシンプルに書いたつもりです。(その代わり数学的に非自明な部分はありません。)

mathtod.online/media/wTSGu0qGe

Julia言語の確率分布のパッケージ Distributions.jl と mpmath の一般超幾何函数を使った素朴な方法(最適化していない)で計算した超幾何分布のCDFでは、前者が後者よりも約1000倍速かったです。

Distributions.jl のソースを見ると速い理由がわかります。(もちろん超幾何函数は使っていない。)

github.com/JuliaStats/Distribu

gist.github.com/genkuroki/e58a
超幾何分布の分布函数を一般超幾何函数 ${}_3F_2$ で表わす公式の数値的確認

en.wikipedia.org/wiki/Hypergeo
のCDFの計算をJulia言語の Distributions.jl パッケージと一般超幾何函数の両方の方法で行って一致することを確認しました。

こういう作業をするときに問題になるのは、一般超幾何函数の計算です。

Pythonのmpmathライブラリに一般超幾何函数 ${}_pF_q$ が見付かったのでそれを使いました。

メジャーな他言語のどれかに必要な函数があればそれをそのまま利用できる感じになっていることはJulia言語の非常に便利な点だと思います。

nbviewer.jupyter.org/gist/genk

上のリンク先に、Julia言語でKdV方程式の数値計算するためのJupyter notebookが置いてあります。

下の方の kdv(t,x,u0) という函数だけを見てくれれば十分です。それ以外の部分は試行錯誤の跡です。

N = 2^8
x = collect(linspace(-10,10,N))
t = collect(linspace(0,9,180))
u0 = -5*sin.(pi*x/10)
u = kdv(x,t,u0);

のようにすれば初期条件 u0 の解 u が得られます。

nbviewer.jupyter.org/gist/genk
複数の確率分布でカイ二乗検定とG検定とFisherの正確検定を比較

私の計算によれば(バグっていないことを期待しているのですが、自信がないので、誰か再検証して!)、

2×2の分割表での独立性の検定において

・(補正無しの古典的な)カイ二乗検定は意外なほど頑健である。

・(補正無しの)G検定は結構不安定で、p値がx以下になる確率がxより大きくなる場合が多く、不当に有意差を出し易くなっているように見える。

・Fisherの正確確率検定は、サンプルを生成している確率分布が超幾何分布の場合以外では、誤差が大きい。p値がx以下になる確率がxよりかなり小さくなる。

添付画像はG検定を使うと有意差を不当に出し易くなっている場合。

mathtod.online/media/v15HHF8hk

2×2の分割表の各行各列の独立性に関するカイ二乗検定、G検定、Fisherの正確検定の比較をシステマティックにやってみました。 Jupyter notebookがリンク先にあります。

nbviewer.jupyter.org/gist/genk
2×2の分割表のカイ二乗検定とG検定とFisherの正確検定の比較

ツイッターの方でも宣伝しておきました。
twitter.com/genkuroki/status/9

私の結論:2×2の分割表における各行各列の独立性の判定には、特別な場合を除いて、G検定やFisherの正確確率検定よりも(補正無しの古典的な)カイ二乗検定を使うべきである。

この結論にはまだ自信がないので、正しい考え方を解説してくれる人がいるとうれしいです。

あと私の計算結果の再検証をしてくれる人がいるとうれしい。私の計算にはバグが残っているかもしれない。

KdV方程式の数値解法で遊んでなかったことを思い出したので、さくっとやっつけました。

twitter.com/genkuroki/status/9


nbviewer.jupyter.org/gist/genk
にソースコードがあります。

全然整理されていない。今回のソースコードはかなりひどい。

添付画像は2ソリトン解。

真の2ソリトン解では1ソリトンの非線形な重ね合わせを初期条件にしないといけないのですが、面倒なので単なる足し算にしてあります。 mathtod.online/media/-iOXyDvHO

同期させる力の大きさを意味するパラメーター $K$ の大きさがある一定の値を超えた瞬間に同期が始まるんですね。 (これが蔵本予想。)

既出ですが、 にJulia言語で気楽に書いたJupyter notebook が置いてあります。

nbviewer.jupyter.org/gist/genk

蔵本予想の数値計算をやってみました。
(「蔵本」の漢字を間違っていた。)

Julia言語のJupyter notebook を に投稿しておきました。

nbviewer.jupyter.org/gist/genk

参考にした文献は

www2.math.kyushu-u.ac.jp/~chib

の最初の方だけ。

添付動画は $N=256$ で $K=0.7 K_c$ の場合です。臨界点より $K$ の値が小さいので同期は起こっていません。 mathtod.online/media/5ebY8if9v

問題は次の行のリンク先に
mathtod.online/@genkuroki/5630

● 100 prisoners problem の解答

に投稿しておきました。

gist.github.com/genkuroki/fa34

せっかくなのでモンテカルロシミュレーションを並列処理にして高速化しておきました。

パッケージの読み込みにも at-everywhere を付けておかないといけないんですね。

高速計算を超気楽に行える点は使う人のストレスを大幅に減らします。

● wikipedia でも詳しく解説されていました。(問題がどう呼ばれているかをそれによって初めて知った。)

en.wikipedia.org/wiki/100_pris

● まだ考えている人は答えを見ない方がいいかも。

置換について理解が深まれば必ず自力で答えを出せるはずの問題。

mathtod.online/@tanutarou/5687

私もn次元球体の体積を使って円周率のモンテカルロ計算をやってみました。

以下は、次元 $n$ と次元 $n$ のモンテカルロ計算で球体の中にランダムに与えた点が入る確率の組です。

(13, 0.000111161)
(14, 3.65762e-5)
(15, 1.16407e-5)
(16, 3.59086e-6)
(17, 1.0756e-6)
(18, 3.13362e-7)
(19, 8.89236e-8)
(20, 2.46114e-8)

繰り返しの回数が $10^5$ の場合には $n >13$ でまったくうまく行かなくなる理由がよくわかります。

次元が大きくなると、表面にほとんどすべての体積が集中するようになります。

私は $10^7$ 回の繰り返しで試してみた。

gist.github.com/genkuroki/df3b

大域的な情報も再現されるperplexity=50の動画も作ってみました。

「みんな、あつまってぇー!」
「あつまったらせいれつしてぇー!」
「あれ、そこちょっとちがうよぉー!」
「ううう、とまった」

最後にドラマがおきます。

これは結構かわいいです!❤

nbviewer.jupyter.org/gist/genk

訂正のための再投稿 mathtod.online/media/haBwWsUqE

t-SNEについて

今度はアニメーションを作ってみました。

nbviewer.jupyter.org/gist/genk

添付動画は3つ作ったうちの一つです。

残りの2つは次の場所で見れます。
twitter.com/genkuroki/status/9

オリジナルの点の配位を一度ランダムに配置して、そこから近くの点は近くに、Kullback-Leibler情報量で測ったオリジナルとの違いが小さくなるように、点の位置を動かして行っているようです。

最初にオリジナルで近くの点どうしがくっつき、その後のt-SNEで認識可能な精度でトポロジーが復元される感じです。

mathtod.online/media/DSHF_fxGH

でt-SNEを試してみた。

t-SNEは高次元空間中の点の分布を近い点が近い点にうつるように低次元空間に非線形射影する技術の一つです。以前から遊んでみたかった。

今回試してみたのはそういう通常の使い方とは違って、2次元を2次元に射影しています。

添付画像は射影前と射影後です。

いつものように で使用したJulia言語のコードを公開しておきました。

gist.github.com/genkuroki/5550

ツイッターにも書きました。
twitter.com/genkuroki/status/9

mathtod.online/media/KAIo_QIxY mathtod.online/media/SB5bEYhrx

引き続き、

gist.github.com/genkuroki/0059
f(2x)f(x+1/2)=f(x)を満たす謎の函数

の話。プラストさんの話はとても面白かったので、Julia kernelのJupyter notebookに簡単にまとめてみました。

最後にコンピューターが示してくれている等式

simplify(F(7,2t)*F(6,t+1//2)/F(6,t)) ==1

は、函数等式 $f(2t)f(t+1/2)=f(t)$ の有限積レベルでの恒等式になっています。

コンピューター上に様々な式を実装しておくと、色々楽をできて、気付かなかったことに気付くということが出て来ます。

無料で使えるソフトなら情報が共有し易いのでとてもありがたいです。

数学は情報を共有してこそ現実世界に影響を与えることができる分野だと思います。

mathtod.online/@1_324718/500943
【こちらで定義した関数 $f(m,n)$ の $m$ を実数に拡張する操作はあったのですが、$f(m,n+1)=f(2m-1,n)/f(2m,n)$ で押し進めていたのは見通しが悪かった様です。】

とプラストさんは言っていたのですが、$f(1,n)$ をプラストさんの方法で計算する方が、$1$ から $2^n$ までの数の $\pm1$ の積を順番に計算するよりも丸め誤差が少ないという利点があるかもしれません。

以下のリンク先の最後の F(20,1.0) (プラストさんの方法で計算)と f(2^18,1.0) のどちらが $1/\sqrt 2$ に近いかを見て下さい。

有理函数として F(20,t) と f(2^18,t)$は等しいのですが、プラストさんの方法で計算した方が丸め誤差が小さいように見えます。

gist.github.com/genkuroki/0059

Juliaとgccの比較について、円周率モンテカルロのような単純計算ではなく、実際に役に立つ数値計算で比較したらどうなるかに興味がある人は多いと思います。

インターネット上で見付かる指数積分函数 Ei(x) のCによる数値計算のコードをほとんどそのままJuliaに忠実に翻訳したものと、もとのCのコードをmingwのgccでコンパイルしたものとの速度比較を実行してみました。

結論:計算の速さはほぼ同じ。

gist.github.com/genkuroki/dabc

JuliaをCの代わりに用いても数値計算の速さ的には何も問題がないというのが結論。

github.com/stevengj/18S096-iap
のような最適化は行っていません。

mathtod.online/media/JfVZMKAwD