■１１．１　適用と乱数の発生方法

表11-1  モンテカルロ法による円の面積計算
from tkinter import *
import random, math
root=Tk(); root.title("円の面積")
canvas=Canvas(root, width=460,height=500, bg="white");canvas.pack()
for i in range(0,41):#縦軸目盛
    Y0=400-i*8
    if (i%10)==0:
        canvas.create_line(40,Y0, 440,Y0,fill="black",width=1)
        canvas.create_text(25,Y0,text="%4.1f" % (i/10),
                           font=("ＭＳ　ゴシック",8))
    else:canvas.create_line(40,Y0, 440,Y0,fill="gray",width=1)
for i in range(0,11):#横軸目盛
    X0=40+i*40
    canvas.create_line(X0,400, X0,80,fill="black",width=1)
    canvas.create_text(X0-4,410,text="%4.0f" % (i*1000),
                       font=("ＭＳ　ゴシック",8))
def randDT(): return 2*(random.random()-0.5)
PI=400-math.pi*80
canvas.create_line(40,PI, 440,PI,fill="blue",width=2)
canvas.create_text(25,PI,text="π",font=("Tims",8))
random.seed(0.5)
T=0; X2=40;Y2=400
for i in range(1,10000):#10,000解の試行
    X=randDT();Y=randDT()
    if (X*X+Y*Y)<1: T+=1#円の内部のときカウントアップ
    X1=X2; Y1=Y2; X2=40+i/25 ;Y2=400-(float(T)/i)*4*80
    canvas.create_line(X1,Y1, X2,Y2,fill="red",width=3)

表11-2  乗積合同法計算
import datetime
seed=[0,0,0]
def initRandom():#■初期化：最初のシード値として時刻データを用いる
    date=datetime.datetime.now()
    Y=date.year; M=date.month ;D=date.day
    H=date.hour;MI=date.minute;S=date.second
    MS=date.microsecond 
    seed[0]=(MS%30000)+1
    seed[1]=((MS+(H*60+MI)*60+S))%30000+1
    seed[2]=((MS+ Y*365+M*31+D)%30000)+1
def Random():#■擬似乱数発生(16ビット整数でも可能になるように記述)
    #以下２行 : Seed1 =171 *Seed1 mod 30269 と同じ
    seed[0] = (seed[0]%177)*171 - int(seed[0]/177)*2
    if seed[0]<0: seed[0] = seed[0] + 30269
    #以下２行 : Seed2 =172 *Seed1 mod 30307 と同じ
    seed[1] = (seed[1] % 176) * 172 - int(seed[1]/176)*35
    if seed[1]<0: seed[1] = seed[1]+ 30307
    #以下２行 : Seed3 =170 *Seed3 mod 30323 と同じ
    seed[2] = (seed[2] % 178) * 170 - int(seed[2]/178)*63
    if seed[2]<0 :seed[2] = seed[2]+ 30323
    #上記３系統の乱数を0～1の実数に変換
    R =(float(seed[0])/30269.0 + float(seed[1])/30307.0+
        float(seed[2])/30323.0)
    return R-int(R)

表11-3 機械エプシロンを求める
def machineEPS() : #機械エプシロンを求める
    E = 1
    while (E+1)!=1: ED = E; E /= 2
    E = 1 + ED; return E - 1
print(machineEPS())

表11-4  機械エプシロンを最小単位とする方法
RandomIndex=0;
RandomTable=[0 for i in range(55)]
def randomTableSet(): #乱数表設定
    for I in range(24):
        R = RandomTable[I] - RandomTable[I + 31]
        RandomTable[I] = R
        if R < 0 : RandomTable[I] = R + 1
    for I in range(24,55):
        R = RandomTable[I] - RandomTable[I - 24]
        RandomTable[I] = R
        if R < 0 : RandomTable[I] = R + 1
def initRandom2(): #乱数初期化
    global RandomIndex
    date=datetime.datetime.now()
    Seed=date.microsecond*1E-6
    RandomTable[54] = Seed; R = machineEPS()
    for I in range(54):
        J = (21 * I) % 55; RandomTable[J] = R
        R = Seed - R
        if R < 0 : R +=1; Seed = RandomTable[J]
    randomTableSet(); randomTableSet(); randomTableSet()
    RandomIndex = 54
def Random2():
    global RandomIndex
    RandomIndex +=1
    if RandomIndex > 54 :randomTableSet(); RandomIndex = 1
    return RandomTable[RandomIndex]

表11-5  カイ二乗検定
def array(N1, N2=0):
    if N2 == 0 :return [0 for i in range(N1)]
    else: return [[0 for i in range(N2)] for j in range(N1)]
def freqRandom(fn, M,N):
    X=array(10)
    for k in range(N):ID=int(fn()*M);X[ID]+= 1
    return X
initRandom2();KAI=array(100)
M = 10; N = 100000; E=N/M
for J in range(100):
    X=freqRandom(Random2, M,N);Kai=0
    for I in range(M): EE = E-X[I]; Kai+=EE*EE/E
    KAI[J]= KaiT=0
for J in range(100): T+=KAI[J]
print("手法2 χ2乗平均 ",T/100)

以下は，表11-5の以下の3行（赤色，下線部）を変更して実行した結果，すなわち乗積合同法のχ2検定です。

---(上略)---
initRandom();KAI=array(100)
---(中略)---
    X=freqRandom(Random, M,N);Kai=0
print("手法1 χ2乗平均 ",T/100)

以下は，表11-5の以下の4行（赤色，下線部）を変更して実行した結果，すなわちPythonに標準的に備わっている乱数のχ2検定です。

---(上略)---
import random
KAI=array(100)#初期化不要
---(中略)---
    X=freqRandom(random.random, M,N);Kai=0
---(中略)---
print("Python乱数 χ2乗平均 ",T/100)

表11-6  連検定
initRandom()
A = Random(); NS = 10000;X=array(100)
for i in range(NS):
    N = 0
    if A >= 0.5 :
        while A >= 0.5: N +=1; A = Random()
    else:
        while A < 0.5: N += 1; A = Random()
    if N <= 100 : X[N-1] += 1
for i in range(100): X[i]/=NS
V=0.5;D=0.5
for i in range(100):
    print(X[i],",",V,",",X[i]-V); V*=D

■１１．２　色々な分布の乱数

表11-7  分布を表示する処理
from tkinter import *
import random, math
root=Tk()
canvas=Canvas(root, width=460,height=160, bg="white");canvas.pack()
#■目盛描画
def scale(XST,YST,XW=2,YW=20,NX=100,NY=6,stepX=10,stepY=1,
          GXST=-1, GYST=0, XS=50, YS=100,Xfm="%4.1f",Yfm="%4.2f"):
    XED=XST+NX*XW;YED=YST-NY*YW
    for i in range(0,NY+1,stepY):
        Y0=YST-i*YW
        canvas.create_line(XST,Y0, XED,Y0,fill="blue",width=1)
        canvas.create_text(XST-15,Y0,text=Yfm % (GYST+i/YS),
                           font=("ＭＳ　ゴシック",8))
    for i in range(0,NX+1,stepX):
        X0=XST+i*XW
        canvas.create_line(X0,YST, X0,YED,fill="blue",width=1)
        canvas.create_text(X0-4,YST+10,text=Xfm % (GXST+i/XS),
                           font=("ＭＳ　ゴシック",8))
def array(N1, N2=0):#配列
    if N2 == 0 :return [0 for i in range(N1)]
    else: return [[0 for i in range(N2)] for j in range(N1)]
#■グラフ描画
def drawDT(R, XST,YST,GXST=-1,GYST=0,MST=10,MED=100,XS=4,YS=10,
           color="red"):
    X0=XST+(R[0][0]-GXST)*XS;Y0=YST-(R[0][1]-GYST)*YS
    for i in range(MST,MED):
        X1=XST+(R[i][0]-GXST)*XS;Y1=YST-(R[i][1]-GYST)*YS
        canvas.create_line(X0,Y0, X1,Y1,fill="red",width=2)
        X0=X1; Y0=Y1
#■分布を求める
def randHisto(fn, M=100,MX=100,N=1000000,RST=6,MD=10,MM=-1):
    R=array(M*2+1,2)
    for k in range(N):
        RR=fn(); ID=int((RR+RST)*M/MD)
        if ID>=0 and ID<M*2: R[ID][1]+=1
    for i in range(M*2+1): R[i][0]=float(i)/MX-MM; R[i][1]/=N;
    return R
#■データ表示
def printDT(R):
    for i in range(len(R)):print("%4.2f, %15.5f" %(R[i][0],R[i][1]))

① 正規分布
def randNomal():#正規分布乱数発生関数
    T=-6
    for i in range(12):T+=random.random()
    return T
scale(40,140)
R=randHisto(randNomal,MM=1,MX=100,RST=20,MD=20)
drawDT(R,40,140,MST=10,MED=201,XS=100,YS=1000)

（ボックス・マラーの方法：結果は正規分布と同じ）
keepRandom=[True,0]#乱数保存用[偶数番目，前回計算値]
def BoxMuller():#Box-Muller 法
    global keepRandom
    if keepRandom[0]:
        T = math.sqrt(-2*math.log(1- random.random()))
        U = 2 * random.random()*math.pi
        keepRandom=[False,T * math.sin(U)]
        return T*math.cos(U)
    else: keepRandom[0]=True; return 
scale(40,140)
R=randHisto(BoxMuller,MM=1,MX=100,RST=20,MD=20)
drawDT(R,40,140,MST=10,MED=201,XS=100,YS=1000)

② 三角分布
def Triangle():#三角分布
    return random.random()-random.random()
scale(40,140)
R=randHisto(Triangle,MM=1,MX=100,RST=1,MD=1)
drawDT(R,40,140,MST=10,MED=201,XS=100,YS=10000)

③ 累乗分布
def powerDist():#累乗分布
    LR=math.log(random.random()); N=2# 二乗
    return math.exp(LR/(N+1))
scale(40,140,YW=10,GXST=0,XS=100,YS=10,NY=10,
      Xfm="%3.1f", Yfm="%3.1f")#累乗分布用
R=randHisto(powerDist,M=101,MD=1,MM=0,
RST=0,MX=10) 
drawDT(R,40,140,GXST=0,MST=0,MED=51,XS=40,YS=13000
（もうひとつの方法）
def powerDist2():#累乗分布
    P=random.random();N=2#二乗
for i in range(N):
R=random.random()
 if R>P: P=R
 return P

④ 指数分布
def expDist():#指数分布
    V=random.random()
    while V==0: V=random.random()
    return -math.log(V)
scale(40,140,YW=10,NY=10, YS=10,GXST=0,
XS=10,Xfm="%3.0f")
R=randHisto(expDist,M=500,MX=500,MD=100,MM=0,RST=0)

⑤ コーシー分布
def CauchyDist():#コーシー分布
    X=1;Y=1
    while X*X+Y*Y>1:
        X=random.random();Y=2*random.random()-1
    return Y/X
scale(40,140,GXST=-6,XW=2,XS=8,NX=100,Xfm="%3.0f")
R=randHisto(CauchyDist,M=600,MX=600,MD=100,MM=0,RST=6)

⑥ ロジスティック分布
def LogisticDist():#ロジスティック分布
    R = random.random(); T=0; P3=math.sqrt(3)/math.pi
    if R > 0 : T = P3* math.log((1-R)/R)
    return T
scale(40,140,GXST=-6,XW=2,XS=8,NX=100,Xfm="%3.0f")
R=randHisto(LogisticDist,M=600,MX=600,
MD=100,MM=0,RST=6)
drawDT(R,40,140,GXST=0,MST=0,MED=73,
XS=5000/3,YS=1000)

⑦ ニ項分布
def binomialDist():#二項分布
    N=100; P=0.5; R = 0
    for i in range(N):
        if random.random()<=P: R+=1
    return R
scale(40,140,GXST=0,XW=2,XS=1,YW=10,
NY=10,NX=100,Xfm="%3.0f")
R=randHisto(binomialDist,M=50,MX=1,
MD=50,MM=0,RST=0)
drawDT(R,40,140,GXST=0,MST=0,MED=101,XS=2,YS=1000)

⑧ ガンマ分布
def GammaDistFun(A):
    A=1.5; E = 2.7182818284590452354
    if A>1 :
        T = math.sqrt(2*A - 1)
        eps=-1
        while eps<=random.random():
            U=-100
            while U<=-50:
                X=0
                while X<=0:
                    XY=2
                    while XY>=1 or X==0:
                        X=random.random()
                        Y=2*random.random()-1
                        XY=X*X+Y*Y
                    Y /= X; X = T*Y + A - 1
                U = (A-1)*math.log(X/(A-1))-T*Y   
            eps=(1 + Y*Y) * math.exp(U)
    else:
        T= E/(A+E);Y=-1; Rnd=random.random()
        while Rnd>=Y:
            if Rnd < T :
                X=0;Y=1; R=random.random()
                if R>0 :
                    X = math.exp(math.log(R)/A)
                    Y = math.exp(-X)
            else:
                R = random.random(); X = 1; Y = 0
                if R > 0 :
                    X = 1-math.log(R)
                    Y = math.exp((A-1)*math.log(X))
            Rnd=random.random()
    return X
def GammaDist():return GammaDistFun(1.5)
scale(40,140,YW=20,NY=5, YS=1000,GXST=0,XS=10,
      Xfm="%3.0f",Yfm="%6.3f")
R=randHisto(GammaDist,M=100,MX=20,MD=5,MM=0,RST=0)
drawDT(R,40,140,GXST=0,MST=0,MED=201,XS=20,YS=4000)

⑨ ベータ分布
#  ガンマ分布乱数を利用する方法
def BetaDistFun(A, B):
    T = GammaDistFun(A)
    return T / (T + GammaDistFun(B))
def BetaDist():return BetaDistFun(1.5,2.0)
scale(40,140,YW=20,NY=5, YS=200,GXST=0,XS=100,
      Xfm="%3.1f",Yfm="%6.3f")
R=randHisto(BetaDist,M=100,MX=200,MD=1,MM=0,RST=0)
drawDT(R,40,140,GXST=0,MST=0,MED=101,XS=400,YS=4000)

#  ガンマ分布乱数を利用しない方法
def BetaDistFun2(A, B):
    XX=0
    while XX>=1 or XX==0:
        X = 0; R = random.random()
        if R > 0 : X=math.exp(math.log(R)/A)
        Y = 0; R = random.random()
        if R > 0 : Y=math.exp(math.log(R)/B)
        XX = X + Y
    return X / (X + Y)
def BetaDist2():return BetaDistFun2(1.5,2.0)

⑩ ポアソン分布
def poissonDistFun(M):#ポアソン分布関数処理
    MM=math.exp(M)*random.random(); R=0
    while MM>1: MM*=random.random();R+=1
    return R
def poissonDist():#ポアソン分布関数
    return poissonDistFun(50)
scale(40,140,GXST=0,XS=1,Xfm="%3.0f")
R=randHisto(poissonDist,MD=100,MM=0,RST=0,MX=1)
drawDT(R,40,140,GXST=0,MST=0,MED=101,XS=2,YS=2000)

⑪ カイ二乗分布
#正規分布を利用する方法
def Kai2DistFun(N):
    T = 0
    for i in range(N):S = randNomal();T +=S*S
    return T
def Kai2Dist():return Kai2DistFun(3)
scale(40,140,YW=20,NY=5, YS=200,GXST=0,XW=1.5, 
NX=120,XS=10, Xfm="%3.0f",Yfm="%6.3f")
R=randHisto(Kai2Dist,M=120,MX=20,MD=12,MM=0,RST=0)
drawDT(R,40,140,GXST=0,MST=0,MED=121,XS=30,YS=4000)

#ガンマ分布の乱数を利用する方法
def Kai2DistFun2(N):
    return 2*GammaDistFun(N/2)
def Kai2Dist2():return Kai2DistFun2(3)

⑫ F分布
def fDistFun(N1,N2):# F分布
    return Kai2DistFun(N1)*N2/(Kai2DistFun(N2)*N1)
def fDist():return fDistFun(3,4)
scale(40,140,YW=15,NY=7, YS=1000,GXST=0,XW=1.5,
 NX=120,      XS=10, Xfm="%3.0f",Yfm="%6.3f")
R=randHisto(fDist,M=120,MX=20,MD=12,MM=0,RST=0)
drawDT(R,40,140,GXST=0,MST=0,MED=121,
XS=30,YS=1500)

⑬ ｔ分布
def tDistFun(N):# t分布
    return randNomal()/math.sqrt(Kai2DistFun(N)/N)
def tDist():return tDistFun(3)
scale(40,140,YW=15,NY=7, YS=1000,GXST=-6,XW=1.5,
NX=120,
      XS=10, Xfm="%3.0f",Yfm="%6.3f")
R=randHisto(tDist,M=120,MX=20,MD=12,MM=0,RST=6)
drawDT(R,40,140,GXST=0,MST=0,MED=121,XS=30,YS=1500)

⑬ ワイブル分布
def weibullDistFun(M):#ワイブル分布
    R = random.random(); RR=0
    if R !=0:
        R = -math.log(R)
        if R !=0: RR= math.exp(math.log(R)/M)
    return RR
def weibullDist2():return weibullDistFun(2)
def weibullDist1():return weibullDistFun(1)
def weibullDist05():return weibullDistFun(0.5)
scale(40,140,YW=10,NY=10, YS=500,GXST=0,XW=2,
      XS=50,Xfm="%3.1f",Yfm="%6.3f")
R=randHisto(weibullDist2,M=100,MX=20,MD=10,MM=0,RST=0)
drawDT(R,40,140,GXST=0,MST=0,MED=21,XS=200,YS=500)
R=randHisto(weibullDist05,M=100,MX=20,MD=10,MM=0,RST=0)
drawDT(R,40,140,GXST=0,MST=0,MED=21,XS=200,YS=500,color="green")
R=randHisto(weibullDist1,M=100,MX=20,MD=10,MM=0,RST=0)
drawDT(R,40,140,GXST=0,MST=0,MED=21,XS=200,YS=500,color="indigo")

⑬ 二次元正規分布
import random, math
def array(N1, N2=0):#配列
    if N2 == 0 :return [0 for i in range(N1)]
    else: return [[0 for i in range(N2)] for j in range(N1)]
def normal2D(M):#二次元正規分布
    S=-1
    while S<=0:
        S=2
        while S>1:
            R1 = 2*random.random() - 1
            R2 = 2*random.random() - 1
            S = R1*R1 + R2*R2
    S = -math.log(S)/S
    R1 = math.sqrt((1+M)*S)* R1
    R2 = math.sqrt((1-M)*S)* R2
    return [R1+R2, R1-R2]
def distribution(M):
    X=array(51,51);N=1000000;NM=50
    for i in range(N):
        P = normal2D(M)
        ID1 = int((P[0] + 5) * NM/10)
        ID2 = int((P[1] + 5) * NM/10)
        if (ID1 >= 0 and ID1 <= NM and
            ID2 >= 0 and ID2 <= NM): X[ID1][ID2]+=1
    for i in range(51):
        for j in range(51):
            X[i][j]/=N
    return X


表11-8  色マップおよび等高線図
from tkinter import *
import random, math
root=Tk()
canvas=Canvas(root, width=540,height=540, bg="white");canvas.pack()
def maxValue(X):
    IX=0;IY=0; MDT=X[IX][IY]
    for i in range(len(X)):
        for j in range(len(X[0])):
            if X[i][j]>MDT:MDT=X[i][j]
    return MDT
def codeCH(CD):
    if CD<0: return 0
    if CD>255: return 255
    return int(CD)
def colorCode(H,MaxV):#値による色設定
    V=H/MaxV
    if   V<0.25: R=0;G= int(127*V*4); B=int(127*(0.25-V)*4)
    elif V<0.5 : R=0; G=(127*(V-0.25)*4+127); B=0
    elif V<0.75: R=int((V-0.5)*4*255); G=255; B=0
    else:        R=255; G=int((1-V)*4*255);B=0
    return "#%02X%02X%02X" %(codeCH(R),codeCH(G),codeCH(B))
def colorMap(X):# 色地図
    maxV=maxValue(X);NX=len(X);NY=len(X[0])
    for i in range(NX):
        ii=i*10+20
        for j in range(NY):
            jj=j*10+20; CD=colorCode(X[i][j],maxV)
            canvas.create_rectangle(ii,jj,ii+10,jj+10, 
fill=CD,outline="")
    canvas.create_text(250,10,text=("最大確率 %9.5f" %maxV))     
def interp(H,Z1,Z2): return ((H-Z1)/(Z2-Z1))
def DrawLine(X1,Y1,X2,Y2, color):
    canvas.create_line(X1+25,Y1+25,X2+25,Y2+25,fill=color )
def contour(Z,j, k, H,color="black"):# １区画の等高線
    Z0 = Z[j][k]; Z1 = Z[j+1][k];
    Z2 = Z[j+1][k+1]; Z3 = Z[j][k+1]
    P01 = (Z0> H) ^ (Z1 > H); P12 = (Z1> H) ^ (Z2 > H)
    P23 = (Z2> H) ^ (Z3 > H); P30 = (Z3> H) ^ (Z0 > H)
    dx=10; dy=10         
    if P01 or P12 or P23 or P30:
        x0=0; y1=0; x2=0; y3=0
        if P01:x0 =dx*(j+interp(H, Z0,Z1))
        if P12:y1 =dy*(k+interp(H, Z1,Z2))
        if P23:x2 =dx*(j+interp(H, Z3,Z2))
        if P30:y3 =dy*(k+interp(H, Z0,Z3))
        yk =dy*k; yk1=dy*(k+1);xj =dx*j;xj1=dx*(j+1)
        if P01 and P12: DrawLine(x0 , yk , xj1, y1 ,color)
        if P12 and P23: DrawLine(xj1, y1 , x2 , yk1,color)
        if P23 and P30: DrawLine(x2 , yk1, xj , y3 ,color)
        if P30 and P01: DrawLine(xj , y3 , x0 , yk ,color)
        if P01 and P23: DrawLine(x0 , yk , x2 , yk1,color)
        if P12 and P30: DrawLine(xj1, y1 , xj , y3 ,color)
def contours(Z, v1, v2,dv):#等高線
    for j in range(50):
        for k in range(50):
            v=v1; ii=0
            while v<(v2+dv):
                color="black"
                if ii %10 ==0 : color="white"
                contour(Z,j, k, v,color=color); v+=dv;ii+=1
def allContours(Z):
    maxV=maxValue(X);dv=maxV/50
    contours(Z, dv, maxV,dv)
X=distribution(0.8)#引数を0.0～1.0の間で指定する
colorMap(X)
allContours(X)

■１１．３　モンテカルロ法の応用

表11-9  重積分へのモンテカルロ法の適用
from tkinter import *
import random, math
root=Tk()
canvas=Canvas(root, width=540,height=160, bg="white");canvas.pack()
---(中略：array, scale,drawDTについては表11-7参照)---
def F(P, Q, R):return 6*(2*P+1)*(3*Q-1)*R
def Rnd():return random.random()
def multiInte():
    V = 0; N = 100000;K=0; X=array(int(N/100),2)
    for i in range(1,N+1):
        V += F(Rnd(), Rnd(), Rnd())
        if (i % 100) == 0: X[K][0]=i; X[K][1]= V/i;K+=1
    return X
random.seed(1.0)
X=multiInte()
scale(40,140,YW=30,NY=4, YS=2,GXST=0,GYST=2,XW=5,NX=50,
      XS=0.0005,Xfm="%3.0f",Yfm="%6.1f")
drawDT(X,40,140,GXST=0,GYST=2,MST=0,MED=1000,
XS=0.0025,YS=60,color="red")

表11-10  混合物質の電気伝導性
from tkinter import *
from tkinter import messagebox
import random, math, time
def array(N1, N2=0, N3=0):
    if N2 == 0 and N3==0: return [0 for i in range(N1)]
    elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)]
    else: return [[[0 for i in range(N3)] for j in range(N2)]
                  for k in range(N1)]
N=100; M=100
A=array(N+2,M+2);D=array(N+2,M+2)
def setCell():
    for i in range(2*N*M):
        X=int(N*random.random()+1);Y=int(M*random.random()+1)
        if A[X][Y]==0 : break
    A[X][Y]=1
def direction(i, j):
    if   j > M                  : return True
    elif j < 1 or i < 1 or i > M: return False
    elif D[i][j]!=0             : return False
    else:
        D[i][j] = 1  #以前来たことがあることを示すフラグ
        if A[i][j]==0 : return False
        R = direction(i, j + 1)
        if R: D[i][j+1]=2; return R
        R = direction(i, j - 1)
        if R: D[i][j-1]=2;return R
        R = direction(i + 1, j)
        if R: D[i+1][j]=2;return R
        R=direction(i - 1, j)
        if R: D[i-1][j]=2;return R
        return R
def judge():
    for i in range(1,N+1):
        for j in range(M+1): D[i][j]=0
    for i in range(1, N+1):
        if A[i][1]!=0:
            if direction(i,1): D[i][1]=2;return True
    return False
def percolation(root):
    for i in range(N*M+1):
        if i % 100 ==0: root.title("粒子数=%d" % i);root.update()
        setCell()
        if judge():
            root.title("粒子数=%d" % i);root.update();return i
    return 0
def dsp(canvas):
    for i in range(1,N+1):
        ii=i*2+10
        for j in range(1, M+1):
            jj=j*2+10
            if A[i][j]!=0:
                CL="#7777FF"
                if D[i][j]==2: CL="#FF0000"
                canvas.create_rectangle(jj,ii,jj+2,ii+2,fill=CL,outline="")    
#■メイン処理
root=Tk(); root.title("パーコレーション(浸透)問題")
canvas=Canvas(root, width=220,height=230, bg="white");canvas.pack()
if percolation(root)!=0: dsp(canvas); root.update()
 
表11-11  酔っ払った人が橋を渡れる確率
from tkinter import *
from tkinter import messagebox
import random, math, time
def array(N1, N2=0, N3=0):
    if N2 == 0 and N3==0: return [0 for i in range(N1)]
    elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)]
    else: return [[[0 for i in range(N3)] for j in range(N2)]
                  for k in range(N1)]
XB=0; YB=0
def dsp(canvas, X,Y):
    global XB,YB
    XX=X*50+10; YY=Y*50+10
    canvas.create_line(XB, YB, XX,YY,fill="red",width=2)
    XB=XX; YB=YY
def dspInit(canvas,X,Y): global XB,YB; XB=X*50+10; YB=Y*50+10
def UnitRandom(root, canvas):#　1回試行
    W = 4; L = 8; X = 0; Y = W / 2; K = 2
    canvas.delete(ALL); PI2=2*math.pi
    canvas.create_rectangle(10, 10, L*50+10,W*50+10,fill="#FFFF77",
                            outline="#7777FF",width=2)
    dspInit(canvas, X,Y)
    while X < L:
        TH = PI2 * random.random()
        X += 0.5 * math.cos(TH); Y += 0.5 * math.sin(TH)
        K+=1; dsp(canvas, X,Y)
        if X < 0 or Y < 0 or Y > W : return False
        root.update()
    return True
def RandomWalk(root, canvas):
    N=10000;M=0; DM=100 # DM:成功したときの終了回数(画面キャプチャ用)
    for i in range(1,N+1):
        if UnitRandom(root, canvas):
            M +=1
            root.title("試行回数 %d, 成功回数 %d(%8.4f)"%(i, M, (M/i)))
            if M==DM: return M/i
        else: root.title("試行回数 %d, 成功回数 %d(%8.4f)"%(i, M, (M/i)))
    return M / N
root=Tk()
canvas=Canvas(root, width=420,height=230, bg="white");canvas.pack()
RandomWalk(root,canvas)

■	条件の変更
X += 0.5 * math.cos(TH); Y += 0.5 * math.sin(TH)
　　　　　　　　　　　↓
X += 0.5 * math.cos(TH); Y += 0.2 * math.sin(TH)

表11-12  落石／土砂流出影響範囲シミュレーション
from tkinter import *
from tkinter import messagebox
import random, math, time
def array(N1, N2=0, N3=0):
    if N2 == 0 and N3==0: return [0 for i in range(N1)]
    elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)]
    else: return [[[0 for i in range(N3)] for j in range(N2)]
                  for k in range(N1)]
KX=80; KY=80; MX=KX+1; MY=KY+1
Myu    =0.12# 粘性抵抗(様々な障害物の影響を代表させるので大きくする)
Grad   =1   # 重力加速度
MyuRand=0.0 # 粘性抵抗の揺らぎ(0～0.25程度:粘性抵抗を負にしない)
DT=1        # 時間刻み
Z=array(MX,MY)
MaxB=40; nBall = -1; Ball=array(MaxB, 5)
for j in range(KY):#■初期条件の設定
    for i in range(KX):
        if i>25: Z[i][j]=0
        else: Z[i][j]=10-i*20/50
def genBall(canvas,KX,KY, CS):#岩石／土砂粒子の生成
    global nBall,MaxB, Ball
    if nBall<(MaxB-1):
        nBall+=1; ii=CS+10; jj=int(KY*CS/2+20)
        Ball[nBall][0]=canvas.create_oval(ii,jj,ii+CS,jj+CS,fill="#FFFFFF")
        Ball[nBall][1]=1           #位置
        Ball[nBall][2]=int(KY/2)
        Ball[nBall][3]=0           #速度
        Ball[nBall][4]=0
def chkC(R):                       #■R/G/B値のチェック
    if R>255: return 255
    if R<0: return 0
    return int(R)       
def dspCell(canvas, Z,KX,KY,CS):   #■地形図をカラーマップで表示
    ZMax = Z[1][1]; ZMin = Z[1][1] #最大・最小値を求める
    for i in range(KX):
        ZMax=max(ZMax, max(Z[i]))
        ZMin=min(ZMin, min(Z[i]))
    ZDX= ZMax - ZMin; B=0
    for j in range(1,KY):
        jj=j*CS+20
        for i in range(1,KX):
            ii=i*CS+10
            V=int(((Z[i][j]-ZMin)/ZDX)*765)#表示色の決定
            if   V<256: G=int(55+200*V/255); R=0
            elif V<512: V-=256; R=int(V)   ; G=255       
            else      : V-=512; R=255      ; G=255-int(V)
            C="#%02X%02X%02X" % (chkC(R),chkC(G),B)
            canvas.create_rectangle(ii,jj,ii+CS,jj+CS,
                                    fill=C, outline="")
def randNomal():
    T=-6
    for i in range(12):T+=random.random()
    return T
def simulate(canvas,CS):
    global nBall, Ball,Z,Myu,minMyu,maxMyu
    X=Ball[nBall][1]; Y=Ball[nBall][2]
    i=int(X); j=int(Y)  #運動方程式計算
    DX=Z[i+1][j]-Z[i-1][j]; DY=Z[i][j+1]-Z[i][j-1]
    GRD2=DX*DX+DY*DY; F=math.sqrt(GRD2)/math.sqrt(1+GRD2)
    TH=math.atan2(DY,DX)+math.pi*randNomal()/6;
    Fx=-F*math.cos(TH);Fy=F*math.sin(TH)
    rMyu=Myu*(1+randNomal()*MyuRand)
    #if rMyu<minMyu: minMyu=rMyu #以下3行チェック用
    #if rMyu>maxMyu: maxMyu=rMyu
    #print("%16.10f, %16.10f, %16.10f " % (rMyu, minMyu, maxMyu))
    Ball[nBall][3]=(Ball[nBall][3]+Fx*Grad*DT)*(1-rMyu)#速度計算
    Ball[nBall][4]=(Ball[nBall][4]+Fy*Grad*DT)*(1-rMyu)
    X=X+Ball[nBall][3]*DT; Ball[nBall][1]=X            #位置計算
    Y=Y+Ball[nBall][4]*DT; Ball[nBall][2]=Y
    ii=X*CS+10;jj=Y*CS+20
    #canvas.coords(Ball[nBall][0], ii,jj,ii+CS, jj+CS)
    Ball[nBall][0]=canvas.create_oval(ii,jj,ii+CS,jj+CS,fill="#FFFFFF")
    return Ball[nBall][3]*Ball[nBall][3]+Ball[nBall][4]*Ball[nBall][4]
#■メイン処理
root=Tk(); root.title("斜面上の落石シミュレーション")
canvas=Canvas(root, width=340,height=340, bg="white");canvas.pack()
N=0; minMyu=99999;maxMyu=-999999
S="粘性抵抗基準値 % 8.4f  揺らぎ %8.4f" % (Myu, MyuRand)
while True:
    dspCell(canvas, Z,KX,KY, 4)
    canvas.create_text(160,10, text=S,font=("ＭＳ　ゴシック",10), fill="blue")
    for k in range(MaxB):
        genBall(canvas,KX,KY,4); root.update()
        while True:
            epsi=simulate(canvas,4); root.update()
            if epsi<0.00001:break
            time.sleep(0.01)
    if not messagebox.askyesno('落石','繰り返しますか'): break

表11-13  分子の初期位置設定例と共通処理
from tkinter import *
from tkinter import messagebox
import random, math, time
def array(N1, N2=0, N3=0):
    if N2 == 0 and N3==0: return [0 for i in range(N1)]
    elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)]
    else: return [[[0 for i in range(N3)] for j in range(N2)]
                  for k in range(N1)]
def genMole():
    global N, R
    for i in range(N):
        R[0][i][0] = random.random() * (Xmax - 2.0) + 1.0
        R[0][i][1] = random.random() * (Ymax - 2.0) + 1.0

表11-14  分子間力のシミュレーション
・
・---(前略) 表11-13参照---
・
Xmax = 10; Ymax = 10; DispR = 10
EPSI = 1.0; SIGMA = 1;MG = 1
N=40                               # 分子の数
DT=0.05                            # 時間刻み
timerInt=10                        # 表示時間間隔
R=array(2,N,2); F=array(N,2)       # R:位置，F:力
V=array(2,N,2); Bname=array(N)     # V:速度, Bname:表示円形
E1 = -48 * EPSI * pow(SIGMA,13)    # 計算用係数
E2 = 24 * EPSI * pow(SIGMA, 7)
alfa= 0.1                          # alfa = M/(3kT) とみなす(M: 分子質量合計)
ID1 = 0; ID2 = 1                   # 実行制御用
Nup=0;Vup=0;Nlf=0; Vlf=0           #衝突回数と衝突速度の絶対値合計
Nbt=0;Vbt=0;Nrt=0; Vrt=0
sCtr=0                             #シミュレーション回数
def genMono(canvas, X0, Y0,R):#■画面上の分子の生成
    return canvas.create_oval(X0-R,Y0-R, X0+R, Y0+R,fill="red")
def init():             #■初期設定
    global N, F, V,Bname,DispR
    for i in range(N):
        Bname[i] = genMono(canvas, 0, 0, DispR)
def dsp(ID1):           #■表示
    global N, canvas, R, Bname, DispR
    RD=DispR
    for i in range(N):
        X = R[ID1][i][0]*20
        Y = R[ID1][i][1]*20
        canvas.coords(Bname[i], X-RD, Y-RD, X+RD, Y+RD)
def force(ID1):         #■力の計算
    global N, F, R, alfa
    for i in range(N):
        F[i][0] = 0; F[i][1] = 0
        for j in range(N):
            if i != j :
                dX = R[ID1][j][0] - R[ID1][i][0]
                dY = R[ID1][j][1] - R[ID1][i][1]
                DR = math.sqrt(dX * dX + dY * dY)
                if DR < 1 : DR = 1
                FF = E1 / math.pow(DR,13) + E2 / math.pow(DR,7)
                F[i][0] = F[i][0] + FF * dX / DR
                F[i][1] = F[i][1] + FF * dY / DR
def velocity(ID1,ID2):  #■速度の計算
    global N, F, V,MG
    for i in range(N):
        V[ID2][i][0] = V[ID1][i][0] + DT * F[i][0]/ MG
        V[ID2][i][1] = V[ID1][i][1] + DT * F[i][1]/ MG
    T = 0 # 温度境界条件
    for i in range(N):
        VX = V[ID2][i][0]; VY = V[ID2][i][1]
        T += VX * VX + VY * VY
    T = math.sqrt(T)* alfa
    if T > 0.00001:
        for i in range(N):
            V[ID2][i][0] /= T; V[ID2][i][1] /= T
def place(ID1,ID2):     #■位置の計算
    global N, V, X
    for i in range(N):
        R[ID2][i][0] = R[ID1][i][0] + DT * V[ID2][i][0]
        R[ID2][i][1] = R[ID1][i][1] + DT * V[ID2][i][1]
def wallHit(ID2):#■壁への衝突
    global Nup,Vup,Nlf, Vlf, Nbt, Vbt,Nrt, Vrt
    for i in range(N):
        if R[ID2][i][0] < (0.5 * SIGMA):
            R[ID2][i][0] = 0.5 * SIGMA
            AV=abs(V[ID2][i][0])
            V[ID2][i][0] = AV; Nlf+=1; Vlf+=AV
        if R[ID2][i][0] > (Xmax - 0.5 * SIGMA):
            R[ID2][i][0] = Xmax - 0.5 * SIGMA
            AV=abs(V[ID2][i][0])
            V[ID2][i][0] = -AV; Nrt+=1; Vrt+=AV
        if R[ID2][i][1] < (0.5 * SIGMA):
            R[ID2][i][1] = 0.5 * SIGMA
            AV=abs(V[ID2][i][1])
            V[ID2][i][1] = AV; Nup+=1; Vup+=AV
        if R[ID2][i][1] > (Ymax - 0.5 * SIGMA):
            R[ID2][i][1] = Ymax - 0.5 * SIGMA
            AV=abs(V[ID2][i][1])
            V[ID2][i][1] = -AV; Nbt+=1; Vbt+=AV
def simulate(ID1, ID2):#■シミュレーション
    global sCtr,Vup,Vlf,Vbt, Vrt, root
    dsp(ID1); force(ID1); velocity(ID1,ID2)
    place(ID1,ID2); wallHit(ID2)
    sCtr+=1;
    print("%d, %7.4f, %7.4f, %7.4f, %7.4f" %
          (sCtr, Vup/sCtr,Vbt/sCtr, Vlf/sCtr,Vrt/sCtr)) 
    root.update()
#■画面コピー用マウスイベント  
def leftMouseDown(event):
    global ID1, ID2
    ID2=1-ID1; simulate(ID1, ID2); ID1=ID2
#■メイン処理
root=Tk(); root.title("分子運動のシミュレーション")
canvas=Canvas(root, width=200,height=200); canvas.pack()
init(); genMole();ID1=0
#左ボタン
#canvas.bind("<Button-1>",leftMouseDown)  #画面コピー用マウス
#root.mainloop()
#画面コピーするときは以下をコメントにする
ID1=0
print("計算回数,上壁, 下壁,左壁, 右壁")
print("0,  0.0000,  0.0000,  0.0000,  0.0000")
while True:
    ID2=1-ID1; simulate(ID1,ID2)
    time.sleep(0.001)
    ID1=ID2

表11-15  ほこりの運動のシミュレーション
・
・---(前略) 表11-13参照---
・
Xmax = 10; Ymax = 10; DispR = 1 #ほこりっぽく見せるために小さくする
EPSI = 1.0; SIGMA = 1;MG = 1
N=100                              # 分子の数
　・
　・---(中略) 表11-13と同じ---
　・
def force(ID1):         #■力の計算
    global N, F, R
    for i in range(N): # 力をランダムな力にする
        F[i][0] = random.random()-0.5
        F[i][1] = random.random()-0.5
　・
　・---(中略) 表11-13と同じ---
　・
def wallHit(ID2):#■壁への衝突
    global Nup,Vup,Nlf, Vlf, Nbt, Vbt,Nrt, Vrt
    for i in range(N):
        if R[ID2][i][0] < (0.5 * SIGMA):
            R[ID2][i][0] = 0.5 * SIGMA
            AV=abs(V[ID2][i][0])
            V[ID2][i][0] = AV * 0.0 # 0.0 以上にすると端に集まりにくい
        if R[ID2][i][0] > (Xmax - 0.5 * SIGMA):
            R[ID2][i][0] = Xmax - 0.5 * SIGMA
            AV=abs(V[ID2][i][0])
            V[ID2][i][0] = -AV * 0.0 # 0.0 以上にすると端に集まりにくい
        if R[ID2][i][1] < (0.5 * SIGMA):
            R[ID2][i][1] = 0.5 * SIGMA
            AV=abs(V[ID2][i][1])
            V[ID2][i][1] = AV * 0.0 # 0.0 以上にすると端に集まりにくい
        if R[ID2][i][1] > (Ymax - 0.5 * SIGMA):
            R[ID2][i][1] = Ymax - 0.5 * SIGMA
            AV=abs(V[ID2][i][1])
            V[ID2][i][1] = -AV * 0.0 # 0.0 以上にすると端に集まりにくい

表11-16  食物連鎖のシミュレーション
# -*- coding: UTF-8 -*-
#----------------------------------------
#■ N巴の食物連鎖(food chain of N cycles)
#----------------------------------------
from tkinter import *
from tkinter import messagebox
import random, math, time
#■配列として用いるリスト宣言
def array(N1, N2=0, N3=0):
    if N2 == 0 and N3==0:return [0 for i in range(N1)]
    elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)]
    else:return [[[0 for i in range(N3)] for j in range(N2)] 
for k in range(N1)]
#■共通データ
MAX = 101   # セル数+1
MD  = 8   # サイクル数
EP  = 0.5   #  敵に食べられる確率。EP=1にすると幾何学模様に近くなる。
dt=array(2,MAX+1,MAX+1)     # 計算用セルデータ
cTab=("#FF0000","#0000FF","#FFFF00","#FF00FF","#FFFFFF","#007755","#FF7700","#770077",
      "#770000","#000077","#777700","#007700","#777777","#00FF77","#FF77FF","#FF0077")
hist=array(201)             # 履歴格納用配列
avg4=array(201)             # 20回平均格納用
dspArea=array(MAX+1,MAX+1)  # 表示セル用配列
MAXMAX=(MAX-1)*(MAX-1)      # 全体個体数
avrN=MAXMAX/MD              # 平均個体数
avLine=MAXMAX/40+30         # グラフ表示平均位置         
cur=0                       # 現在のセルを示す
Vline=199                   # 200回毎のグラフ位置を示す縦線の位置
for i in range(MAX+1):      # 計算用セルデータの初期設定
    for j in range(MAX+1):
        dt[cur][i][j]= int(random.random() * MD + 1)
for i in range(201):        # 履歴、20回平均初期設定
    hist[i]=avrN;avg4[i]=avrN
#■個体数をグラフ表示する生物1の個体数のカウント
def count1():
    C=0
    for i in range(MAX):
        for j in range(MAX):
            if dt[cur][i][j]==1: C+=1
    return C
#■履歴データ設定
def setHist(V):             
    for i in range(1,200):
              hist[i-1]=hist[i]; avg4[i-1]=avg4[i]
    hist[199]=V 	    # 値の登録
    av=0; j=199#過去20回平均
    for i in range(20): av+=hist[j]; j-=1
    avg4[199]= int(float(av)/20.0+0.5)
#■グラフ表示用XY座標設定    
def setHistXY(hist, i): 
    global avrN, avLine
    return [i+20, avLine-(hist[i]-avrN)/50]
#■最初の表示
def initDraw():
    global canvas, MAX,HSgraph, AVgraph, Vline, VlineX,txtDsp
    canvas.delete("DT")
    for i in range(1,MAX):
        for j in range(1,MAX):
            X=i*2+20; Y=j*2+10
            dspArea[i][j]=canvas.create_rectangle(X, Y, X+2, Y+2, fill=cTab[dt[cur][i][j]-1],outline="", tag="DT")
    # 生物1の個体数変化の表示
    V=count1(); setHist(V)
    canvas.create_rectangle(20,220,220,340,
                            fill="#FFFFFF", outline='black',
width=2, tag="DT")
    XY=[]# 過去4回平均
    for i in range(200): XY.append(setHistXY(avg4, i))
    AVgraph=canvas.create_line(XY, fill="#007700", tag="DT") 
    XY=[]# 履歴表示
    for i in range(200): XY.append(setHistXY(hist, i))
    HSgraph=canvas.create_line(XY, fill=cTab[0], tag="DT")
    # 3,000, 5,000, 7,000の線を表示
    canvas.create_line(20,320,220,320, fill="#000077", tag="DT")
    canvas.create_line(20,280,220,280, fill="#000077", tag="DT")
    canvas.create_line(20,240,220,240, fill="#000077", tag="DT")
    VlineX = canvas.create_line(Vline+20,240,Vline+20, 340, 
fill="#000077", tag="DT")
    Vline-=1
    if Vline<0 : Vline=199
    txtDsp=canvas.create_text(120, 230, 
text=("%d 個体数変化(赤) 20回平均(緑) %d" % (V, avg4[199])),
                              font=('sans-serif', '8'), tag="DT")
#■2回目以降のセル表示(セルの色変更)
def draw():
    global canvas, MAX,HSgraph, AVgraph, Vline, VlineX,txtDsp
    for i in range(1,MAX):
        for j in range(1,MAX):
            X=i+20; Y=j+10
            CC=cTab[dt[cur][i][j]-1] 
            if CC != canvas.itemcget(dspArea[i][j], option='fill'):
                canvas.itemconfig(dspArea[i][j], fill=CC)
            j+=1
        i+=1
#■2回目以降のグラフ表示(表示が変わる部分のみ）
def drawGraph():
    global canvas, MAX,HSgraph, AVgraph, Vline, VlineX,txtDsp
    V=count1(); setHist(V)
    XY=[]# 過去4回平均
    for i in range(200): XY.append(setHistXY(avg4, i))
    canvas.delete(AVgraph)
    AVgraph=canvas.create_line(XY, fill="#007700", tag="DT")
    XY=[]# 履歴表示
    for i in range(200): XY.append(setHistXY(hist, i))
    canvas.delete(HSgraph)
    HSgraph=canvas.create_line(XY, fill=cTab[0], tag="DT")
      canvas.coords(VlineX, Vline+20,240,Vline+20, 340)
    Vline-=1
    if Vline<0 : Vline=199
    canvas.itemconfig(txtDsp, 
text=("%d 個体数変化(赤) 20回平均(緑) %d" % 
(V, avg4[199])))
#■食物連鎖のシミュレーション    
def simChain():
    global cur; nxt=1-cur
    for i in range(1,MAX):
        for j in range(1,MAX):
            dt[nxt][i][j]=dt[cur][i][j]
            if random.random()< EP : 
                for k in range(1,MD+1):
                    KF = k + 1
                    if k== MD: KF = 1
                    if dt[cur][i][j]==k :
                        if(dt[cur][i  ][j-1]==KF or
dt[cur][i  ][j+1]==KF or
                           dt[cur][i-1][j ]==KF or
 dt[cur][i+1][j  ]==KF) : dt[nxt][i][j]=KF
cur=nxt
#■メイン処理        
root=Tk();root.title("N巴の食物連鎖(N= %d)" % MD )
canvas=Canvas(root, width=260,height=400)
canvas.pack();simChain();initDraw()
root.update()
while True:
    simChain(); draw(); drawGraph()
    time.sleep(0.01); root.update()