問題文
座標平面上で $x$ 座標と $y$ 座標がいずれも整数である点を格子点という。格子点上を次の規則に従って動く点 $\mathrm{P}$ を考える。
(a) 最初に、点 $\mathrm{P}$ は原点 $\mathrm{O}$ にある。
(b) ある時刻で点 $\mathrm{P}$ が格子点 $(m,n)$ にあるとき、その $1$ 秒後の点 $\mathrm{P}$ の位置は、隣接する格子点 $(m+1,n),$ $(m,n+1),$ $(m-1,n),$ $(m,n-1)$ のいずれかであり、また、これらの点に移動する確率は、それぞれ $\frac{1}{4}$ である。
- 点 $\mathrm{P}$ が、最初から $6$ 秒後に直線 $y=x$ 上にある確率を求めよ。
- 点 $\mathrm{P}$ が、最初から $6$ 秒後に原点 $\mathrm{O}$ にある確率を求めよ。
（2017 東京大学 理系第２問）
- $\frac{5}{16}$
- $\frac{25}{256}$
※ この答えは学校側が公表したものではありません。個人が作成した非公式の答えです。
コード
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as ani
from matplotlib.colors import Normalize
#定数の準備
left = -7
right = 7
bottom = -7
top = 7
delta = 0.1
cm_name = "Wistia"
fdic = dict(color="darkgreen",ha="center",va="center",animated=True)
#紙の準備
fig = plt.figure(facecolor=plt.get_cmap(cm_name)(0))
fig.subplots_adjust(bottom=0.0,top=1.0)
fig.suptitle("2017東大数学 理系第２問",
color="0.5",ha="right",x=0.96,y=0.96)
ax = fig.add_subplot(111)
ax.set_xlim(left,right)
ax.set_ylim(bottom,top)
ax.set_aspect("equal")
ax.axis("off")
#xy座標の準備
x = np.arange(left,right+delta,delta)
y = np.arange(bottom,top+delta,delta)
X,Y = np.meshgrid(x,y)
#時間ごとの確率分布を作成
Pn = []
for n in range(7):
P = np.zeros((13,13))
for i in range(n+1):
for j in range(n+1):
p = math.comb(n,i)*math.comb(n,j)/4**n
a = i-j
b = n-i-j
P[a+6,b+6] = p
Pn.append(P)
#確率マップを描くための関数
def get_pmap(F):
pmap = ax.contourf(X,Y,F,cmap=cm_name,
norm=Normalize(vmin=1/3.0,vmax=1/2.5))
return pmap
#0.0~5.9秒後
imgs = []
for t in range(60):
tq = t//10
tr = (t%10)/10
Ft = np.zeros_like(X)
for i in range(13):
for j in range(13):
if Pn[tq][i,j] != 0:
a = i-6
b = j-6
sigma_l = (1-tr)*Pn[tq][i,j]+tr*Pn[tq+1][i-1,j]
sigma_r = (1-tr)*Pn[tq][i,j]+tr*Pn[tq+1][i+1,j]
sigma_b = (1-tr)*Pn[tq][i,j]+tr*Pn[tq+1][i,j-1]
sigma_t = (1-tr)*Pn[tq][i,j]+tr*Pn[tq+1][i,j+1]
Ft += np.exp((-(X-a+tr)**2-(Y-b)**2)/(2*sigma_l**2))
Ft += np.exp((-(X-a-tr)**2-(Y-b)**2)/(2*sigma_r**2))
Ft += np.exp((-(X-a)**2-(Y-b+tr)**2)/(2*sigma_b**2))
Ft += np.exp((-(X-a)**2-(Y-b-tr)**2)/(2*sigma_t**2))
P = get_pmap(Ft)
#整数秒のときには場合の数も表示します
if tr==0:
for k in range(10):
if k == 5:
for l in range(30):
if tq == 0:
txt = []
else:
txt = [ax.text(-1.5-tq/2,1.5+tq/2,chr(0xFF10+tq)+"秒後",fontsize=16,fontdict=fdic)]
txt+= [ax.text(1.5+tq/2,-1.5-tq/2,"$\\times\\,\\frac{1}{"+str(4**tq)+"}$",fontsize=16,fontdict=fdic)]
for i in range(tq+1):
for j in range(tq+1):
a = i-j
b = tq-i-j
latex = '$'+str(math.comb(tq,i)*math.comb(tq,j))+'$'
txt += [ax.text(a,b,latex,fontsize=15-tq,fontdict=fdic)]
imgs.append(P.collections+txt)
else:
imgs.append(P.collections)
else:
imgs.append(P.collections)
#6.0秒後
for k in range(70):
F = np.zeros_like(X)
for i in range(13):
for j in range(13):
if Pn[6][i,j] != 0:
a = i-6
b = j-6
sigma = Pn[6][i,j]
F += 4*np.exp((-(X-a)**2-(Y-b)**2)/(2*sigma**2))
P = get_pmap(F)
if k >= 5:
txt = [ax.text(-4.5,4.5,chr(0xFF10+6)+"秒後",fontsize=16,fontdict=fdic)]
txt+= [ax.text(4.5,-4.5,"$\\times\\,\\frac{1}{"+str(4**6)+"}$",fontsize=16,fontdict=fdic)]
for i in range(7):
for j in range(7):
a = i-j
b = 6-i-j
latex = '$'+str(math.comb(6,i)*math.comb(6,j))+'$'
txt += [ax.text(a,b,latex,fontsize=9,fontdict=fdic)]
imgs.append(P.collections+txt)
else:
imgs.append(P.collections)
#上演
mov=ani.ArtistAnimation(fig, imgs, 100)
plt.show()