黒木玄 Gen Kuroki Profile picture
Dec 6, 2021 16 tweets 8 min read Read on X
以前書いた #Julia言語 版HMC(Hamiltonian Monte Carlo)のサンプルコード

ポテンシャル函数φ(x)から、確率分布p(x)=exp(-φ(x))/Zのi.i.d.サンプルを生成する方法の1つ。

60行程度しかない。 Image
#Julia言語 leapfrog法でHamiltonの正準方程式を解くので、leapfrog法のために必要な情報をLFProblem型の変数に格納し、それをHMC函数に渡すと分布p(x)=exp(-φ(x))/Zのサンプルを返してくれる。

そういうシンプルなコードになっています。

nbviewer.org/github/genkuro… Image
#Julia言語 この手の問題では、ポテンシャル函数φ(x)がパラメータに依存している場合が多いので、ポテンシャル函数はφ(x, param)の形式の函数で与える仕様になっています。

だから、HMC函数および関連の函数にはポテンシャル函数を決めるためのparamを渡す必要があります。

nbviewer.org/github/genkuro… Image
#Julia言語 Juliaでは(constも含めて)グローバル変数を引数を経由せずに函数の中で使うことは損になります。計算速度が落ちたり、柔軟な試行錯誤の可能性が潰されてしまう。

それを防ぐ最も単純な方法は「グローバル変数を使うときには毎回函数に引数として渡す」です。トートロジカルな方法が良い!
#Julia言語 しかし、複数の函数に毎回引数として渡される変数達が多いと面倒になります。

だから、複数の函数に毎回引数として渡される情報達を1つ(もしくは高々2〜3個)の変数にまとめるとよいです。

そのために自前でstructを定義して利用することが多いです。
#Julia言語 コードは

* 問題を記述するための情報を格納するとオブジェクトのコンストラクタ

* 問題を記述する情報を渡すと問題を解いてくれる函数

の形式で整理すると分かりやすくなることが多いです。
#Julia言語 日本に限らず、

* 問題を記述する情報を沢山のグローバル変数にべた書き

* 長大なmain函数で問題を解く

のようなスタイルのコードを書くスタイルの悪しき教育を受けた人達は多く、上で解説したスタイルなどで整理されたコードの書き方が普及するとよいと思います。
#Julia言語 HMCの短いサンプルコード

nbviewer.org/github/genkuro…

で参考になると思われる点は他にもあります。

まず、JuliaにおけるStaticArrays.jlを用いた計算の効率化。Juliaには多彩な配列の型があり、使いこなすと楽に見易く高速なコードを書けます。続く Image
#Julia言語 次に、自動微分!

Hamiltonの正準方程式ではポテンシャル函数φ(x, param)のxに関する導函数(gradient)が必要になります。

手計算でgradientを計算して、それを函数として実装して、問題を記述する情報が格納された変数に格納するのは面倒です。続く#

nbviewer.org/github/genkuro… Image
#Julia言語 φ(x, param)を与えるだけで、そのgradientを丸め誤差を除いて正確にかつ高速に自動的に計算してくれるならば、人間側はgradientを手計算する作業から解放されます。

nbviewer.org/github/genkuro… のサンプルコードでは実際にそれをForwardDiff.jlによる自動微分で実現しています! Image
#Julia言語 シンプルなサンプルコードだけでどれだけのことができるのか?

例1. φ(x) = (x₁² + x₁x₂ + x₂²)/2の場合の分布p(x)=exp(-φ(x))/Z (これは2次元の正規分布になる)のサンプルの生成。

問題を記述するデータを格納した変数lfを

lf = My.LFProblem(2, φ)

で作っています。 Image
#Julia言語 例2. φ(x) = a(x₁² - 1)² の場合の分布p(x)=exp(-φ(x))/Zのサンプルの生成。a=3,4,5,6,7,8の場合を計算しています。

aが大きいほど、2つの山が強く分離されるようになり、偏りのないサンプル生成が難しくなります。 Image
#Julia言語 例3. 正規分布モデルでのベイズ統計

HMC法はベイズ統計での事後分布のサンプルの構成でも役立つ。

面倒なので事前分布はフラット事前分布にしてしまっています(手抜き)。その場合のポテンシャル函数φ(x)はモデルのパラメータxの対数尤度函数の-1倍になります。

nbviewer.org/github/genkuro… Image
#Julia言語 サンプルコード

nbviewer.org/github/genkuro…

はJuliaで効率的に計算するための基礎になる「型安定性」と「アロケーションの削減」もきっちり実現しています。

こういう点も参考になると思います。 Image
#Julia言語 ループの内側でrandn(rng, n)で長さnの配列の乱数を生成すると、毎回その分だけにメモリ割当が発生し、計算効率が悪化してしまいます。

サンプルコードでは randn!(rng, vtmp) の形式で事前割当された配列vtmpに乱数を書き込むことによってそれを防いでいます。

nbviewer.org/github/genkuro… Image
#Julia言語 注意・警告: StaticArrays.jlの使用は長さ30程度の配列を扱う場合程度までは有効ですが、サイズがそれを超えると大幅に遅くなる場合があるので要注意です。

StaticArrays.jlが使えないほど大きな問題では、in-place計算でメモリ割当を削減する必要があります。

• • •

Missing some Tweet in this thread? You can try to force a refresh
 

Keep Current with 黒木玄 Gen Kuroki

黒木玄 Gen Kuroki Profile picture

Stay in touch and get notified when new unrolls are available from this author!

Read all threads

This Thread may be Removed Anytime!

PDF

Twitter may remove this content at anytime! Save it as PDF for later use!

Try unrolling a thread yourself!

how to unroll video
  1. Follow @ThreadReaderApp to mention us!

  2. From a Twitter thread mention us with a keyword "unroll"
@threadreaderapp unroll

Practice here first or read more on our help page!

More from @genkuroki

Jun 13
#統計 いつも言っていることをそのまま書きます。長めのスレッドになります。

以下スクショによるスライドの引用は より。赤字と青字は私による書き込みコメント。

まず、p.12について。詳しい解説に続く。 speakerdeck.com/shuntaros/jia-…

Image
#統計 「違いがない」の型の帰無仮説のP値をnull P値と呼びます。

null P値は「違いは○○である」の型の仮説に関する無数のP値の特別な場合で、null P値へのこだわりは悪しきnullismである云々とGreenlandさんは言っています。

biostat.ucdavis.edu/sites/g/files/…
Image
#統計 平たく言えば、「違いがない」の型の帰無仮説を「null P値<α」という条件によって棄却して「違いはある」という結論を出すためにP値を単純に使うことはP値の誤用の典型例であり、科学のプロセスを害しています。

biostat.ucdavis.edu/sites/g/files/…
Image
Read 36 tweets
Jun 18, 2023
#統計 念の為のコメント

1️⃣「t検定の使用が適切なためには、母集団が正規分布に従っていることが必要である」という考え方は誤り。

2️⃣「Wilcoxonの順位和検定=Mann-WhitneyのU検定であれば、無条件使用は適切である」という考え方も誤り。

以上の誤りを信じている人達をよく見る。続く
#統計

1️⃣「t検定の使用が適切なためには、母集団が正規分布に従っていることが必要である」という考え方は誤り。

これについてはツイッター上で繰り返し非常に詳しく解説して来ました。

ツイログ検索

twilog.togetter.com/genkuroki/sear…
#統計

2️⃣「Wilcoxonの順位和検定=Mann-WhitneyのU検定であれば、無条件使用は適切である」という考え方も誤り。

これについてもツイッター上で繰り返し非常に詳しく解説して来ました。

ツイログ検索

twilog.togetter.com/genkuroki/sear…
Read 40 tweets
Jun 17, 2023
#数楽 ℤ[√2]やℤ[√3]はEuclid整域なのでPIDでUFDになるので、ℤ[√2]やℤ[√3]係数の多項式の √2や√3が出て来る因数分解の問題も既約元の積に分解する問題として意味を持ちます。続く
#数楽 ただし、整数dに関する√dが出て来る場合には、既約元の積への分解は因子の可逆元倍と順序の違いを無視しても一意的でなくなる場合が出て来ます。

実はそういうところに面白い数学が隠れている!
#数楽 整数の平方根が出て来る因数分解もちょっと話題になっていますが、その話はとてつもなく面白い数学の話に繋がっています!

中学生であっても思いつきそうな話の中にも素晴らしい数学が隠れています!
Read 20 tweets
Jun 16, 2023
東工大出身者のような理系の人達が、上野千鶴子が自閉症の母親原因説を唱えるくらい科学的に無能でかつ優しさに欠けた人物であることぐらいは知っておいた方が、我々の社会はよくなる可能性が高まると思います。

有名かつ有力になってしまった人物はたとえク○であっても無視できなくなる。
上野千鶴子は、自閉症の原因について母子密着説を唱えていたのですが、それが誤りであることが定説になっていることを指摘された後には、定説と上野千鶴子的なトンデモ説を平等に扱うという態度を取りました。

上野千鶴子の自分が苦しめた人達への態度は真にあきれるものでした。
上野千鶴子的な活動家は科学的無知と優しさに欠けた態度の両方の力を行使していました。

そういうことを許す伝統が現代においても人々の苦しみの源泉の1つになっているわけです。
Read 6 tweets
Jun 15, 2023
私は、環論を学ぶまで、重根もしくは重解の概念を十分に理解できた感じがしてなかったです。(代数)方程式の概念も同様。

実数体上の方程式x²=0は環

A = ℝ[x]/(x²)

で表現されます。これと方程式x=0に対応する環

ℝ[x]/(x)

は異なる。環論を使えば方程式x²=0とx=0を明瞭に区別できます。
環k上の環Aで表現された方程式のk上の環Bでの解集合はk上の環準同型全体の集合

Hom_{k-ring}(A, B)

で表現されます。例えば、集合として、

Hom_{ℝ-ring}(ℝ[x,y]/(x²+y²-1), ℝ) ≅ {(x,y)∈ℝ²|x²+y²=1}.
そして、以上のような代数方程式の表現になっている環の話について前もって知っておいた方が、環論の勉強はしやすいように思えます。
Read 6 tweets
Jun 15, 2023
以下のリンク先スレッド中にも書きましたが、

* 最初に共通の定数因子を括り出すと、その後の計算が楽になる場合がある。

と教えるようにして、

* 共通の定数因子を括り出していなくても、目くじらをたてない。

という教え方にすればよいと思いました。
教科書に従って「a(3x-6y)は誤りで、3a(x-2y)が正解だ」と安易に教えてしまった中学校の数学の先生は

 数学の先生なのに
 教科書通りにおかしなことを教えて
 ごめんなさい

と言って欲しいです。数学では教科書の内容を正しいと信じてはいけない。数学はそういうものだと大学で習っているはず。
数学を教えていれば、細かい条件を言い忘れるというような失敗は日常茶飯事のはずです。

人間だから仕方がないです。

大したことではないので、よりクリアになるように訂正すればよいと思います。
Read 18 tweets

Did Thread Reader help you today?

Support us! We are indie developers!


This site is made by just two indie developers on a laptop doing marketing, support and development! Read more about the story.

Become a Premium Member ($3/month or $30/year) and get exclusive features!

Become Premium

Don't want to be a Premium member but still want to support us?

Make a small donation by buying us coffee ($5) or help with server cost ($10)

Donate via Paypal

Or Donate anonymously using crypto!

Ethereum

0xfe58350B80634f60Fa6Dc149a72b4DFbc17D341E copy

Bitcoin

3ATGMxNzCUFzxpMCHL5sWSt4DVtS8UqXpi copy

Thank you for your support!

Follow Us!

:(