Juliaでブラウン運動① - Bioinformathematics

Pythonより速く、文法的にも数学に近く書きやすい言語として人気急上昇中のJuliaの練習がてらにブラウン運動を発生させてみた。julia 1.0.1を使用している。

正規分布を元にBrown運動を発生させようと思ったのでこのパッケージを使用。

まずはパッケージのinstallから。

julia> Pkg.add("Distributions")
 Resolving package versions...
 Installed QuadGK ─────────── v2.0.2
 Installed PDMats ─────────── v0.9.5
 Installed Rmath ──────────── v0.5.0
 Installed StatsFuns ──────── v0.7.0
 Installed Compat ─────────── v1.4.0
 Installed BinaryProvider ─── v0.5.3
 Installed JSON ───────────── v0.20.0
 Installed Distributions ──── v0.16.4
 Installed Arpack ─────────── v0.3.0
 Installed SpecialFunctions ─ v0.7.2
…

上記のように表示される。また、Compatも同様に

julia> Pkg.add("Compat")
 Resolving package versions...
  Updating `C:\Users\***\.julia\environments\v1.0\Project.toml`
  [34da2185] + Compat v1.4.0
  Updating `C:\Users\***\.julia\environments\v1.0\Manifest.toml`
 [no changes]

とインストールする。さらにusingを用いて必要なパッケージを宣言しよう。

julia> using Compat, Random, Distributions
[ Info: Recompiling stale cache file C:\Users\***\.julia\compiled\v1.0\Compat\GSFWK.ji for Compat [34da2185-b29b-5c13-b0c7-acf172513d20]
[ Info: Precompiling Distributions [31c24e10-a181-5473-b8eb-7969acd0382f]

ではまず簡単な一次元ブラウン運動から考えてみよう。t=0にx=0にいる粒子がT秒間一次元ブラウン運動を続けた際のt=Tにおける位置はN(0,T)に従うとされる。そのような乱数を1000個生成してみよう。

Random.seed!(123) # 乱数のシードを与える
T=100 #時間
d = Normal(0,sqrt(T)) #正規分布 N(0,T)、引数として与えるのはσなのでsqrtによって平方根をとる
x = rand(d,1000) #分布dに従う乱数を1000個生成し配列xとする

さてここで生成されたデータをGadFlyを用いてヒストグラムにしてみてみよう。

using DataFrames,Gadfly
df = DataFrame( 
    x = x
) #先ほどのxをデータフレームに入れる
plot(df, x="x", Geom.histogram(bincount=50)) #プロット

省略するがDataFramesやGadflyもinstallが必要ならPkg.addを使ってinstallしよう。さてこれで以下のような図が表示されれば成功である。

f:id:rikein12:20181128202108p:plain — N(0,T)

もちろん乱数を使用しているので厳密で同じである必要はないがおおまかに正規分布になる（当たり前）。

もうちょっと高度な内容は後日書くことにする