ランダムウォークの数理

二つ(+1,-1)の状態をとるランダムウォークは二項分布 $_{10} C _r (\frac{1}{2})^r (\frac{1}{2})^{10-r} = _{10} C _r (\frac{1}{2})^{10}$ に従う

参考文献

[1] 「なるほど確率論   著者:村上雅人」
[2] 「ランダム・ウォーク 乱れに潜む不思議な現象   著者:津野義道」

In [ ]:
times = 10
state = 0
ran_num = 0.0
for i = 1:times
    ran_num = rand()
    if ran_num >= 0.0 && ran_num <= 0.5
        state += 1
    elseif ran_num  > 0.5 && ran_num<= 1.0
        state -= 1
    end
    print("$state ,")
end
In [27]:
state
Out[27]:
2

約20%の確率で6の状態にいる.

In [3]:
# 2項分布の計算
((factorial(10)/(factorial(10-6)))/factorial(6)) * (1/2)^(6) * (1/2)^4
Out[3]:
0.205078125
In [4]:
using PyPlot
using Plots
In [5]:
function binomial_dist(n)
    state = 0
    global X = []
    for i = 0:n
        ran_num = rand()
        if ran_num >= 0.0 && ran_num <= 0.5
            state += 1
        elseif ran_num  > 0.5 && ran_num<= 1.0
            state -= 1
        end
        #print("$state ,")
        # 組み合わせの計算をしているところ
        var =  (1 * i)+(-1 * (n - i))
        #print("var: $var , \n")
        # mpmathとかのインポートエラーになるのでダメ
        # p = ((factorial(big(n)))/(factorial(big(n-abs(var))))/factorial(big(abs(var)))) * (1/2)^(abs(var)) * (1/2)^(n-abs(var))
        p = binomial(n,abs(var)) * (1/2)^(abs(var)) * (1/2)^(n-abs(var))
        #p = ((factorial(n)/(factorial(n-abs(var))))/factorial(abs(var))) * (1/2)^(abs(var)) * (1/2)^(n-abs(var))
        #print("p: $p \n")
        push!(X,p)
    end
 end
Out[5]:
binomial_dist (generic function with 1 method)
In [6]:
binomial_dist(15)

きちんと二項分布になっている.  このシミュレーションでは二項分布 $_{15} C _r (\frac{1}{2})^r (\frac{1}{2})^{15-r} = _{15} C _r (\frac{1}{2})^{15}$ に従う

In [7]:
histogram(X)
Out[7]:
0.00 0.05 0.10 0.15 0.20 0 2 4 6 8 10 12 y1

例えば-10の状態と+10の状態では15C10と15C10で同じ確率になるので二つの凸が発生する.

In [8]:
PyPlot.plot(X)
Out[8]:
1-element Array{PyCall.PyObject,1}:
 PyObject <matplotlib.lines.Line2D object at 0x12df06f90>
In [9]:
@show ims = []
@time for i = 1:50
    binomial_dist(i)
    push!(ims,X)
end
ims = [] = Any[]
  0.000512 seconds (1.57 k allocations: 59.375 KiB)
In [10]:
@show ims = Any[]
@time for i = 1:50
    binomial_dist(i)
    push!(ims,X)
end
ims = Any[] = Any[]
  0.000624 seconds (1.57 k allocations: 59.375 KiB)
In [11]:
function binomial_repitition()
    @show ims = Any[]
    @time for i = 1:50
        binomial_dist(i)
        push!(ims,X)
    end
end
Out[11]:
binomial_repitition (generic function with 1 method)
In [12]:
# functionで囲ったほうが早くなる.
@time binomial_repitition
  0.000009 seconds (16 allocations: 1.031 KiB)
Out[12]:
binomial_repitition (generic function with 1 method)

1から50までの二項分布をプロットする.

In [13]:
using PyPlot
using PyCall
@pyimport matplotlib.animation as anim

########## display file
function displayfile(mimetype, filename)
    open(filename) do f
        base64text = base64encode(f)
        display("text/html", """<img src="data:$mimetype;base64,$base64text">""")
    end
end

########## initialize the frame
#xs = [ims[i][1] for j in 1:2, i in 1:length(ims)]
#ys = [ims[i][2] for j in 1:2, i in 1:length(ims)]
#xmin, xmax = minimum(xs)-0.5, maximum(xs)+0.5
#ymin, ymax = minimum(ys)-0.5, maximum(ys)+0.5
function initframe()
    xlim(xmin, xmax)
    ylim(ymin, ymax)
    grid(ls=":")
end

########## update the frame
function updateframe(t)
    # clf() #←全部消して書き直す場合にはコメントを外す
    #---------- フレームのアップデート開始
    if t  1
        X = ims[t]
        PyPlot.plot(X, lw=1.0)
    end
    title("t = $t")
    #---------- フレームのアップデート終了
    # plot() #←念のためのおまじない. []と書いてもよい. いらない場合も多い.
end

########## filename
filename = "randomwalk.gif"

######### Construct Figure and Plot Data
fig = figure(figsize=(6.4, 4.8))

# Create the animation object by anim.FuncAnimaton
frames = collect(0:length(ims))
#@time myanim = anim.FuncAnimation(fig, updateframe, init_func=initframe, frames=frames, interval=100)
@time myanim = anim.FuncAnimation(fig, updateframe, frames=frames, interval=500)
# Convert it to an animated GIF file by Imagemagick.
@time myanim[:save](filename, writer="imagemagick")

sleep(0.1)

# Display the movie.
displayfile("image/gif", filename)

# Don't display the axes of figure().
clf()
  0.744225 seconds (136.02 k allocations: 7.230 MiB)
 11.131688 seconds (91.69 k allocations: 4.609 MiB)