模造３（状態の変化）
このページはマス旅の一部です。
今回は「状態の遷移」のシミュレーションです。
課題１：シェルピンスキーのギャスケットをセルオートマトンで実現してください。
シェルピンスキーのギャスケットはフラクタルの1種と言えます。漸化式のロジックで書けます。
パスカルの三角形をかいて、奇数を１、偶数を０と数字で塗り分けることでもかけますね。
1次元セルオートマトンは、画面にセルが並んでいるとしたときに自分と左右の隣の３セルの状態が
次の時刻の自分の状態を決めるというルールをセルに連続適用するものです。
1次元セル列を更新する代わりに、
セル列を時刻が１つ進むたびに次の行に追加すると2次元図がかけますね。
2進法で左、自分、右の状態を０か１で表示すると、
2進数の000から111までの8通りの3セル状態があります。
これを10進数でかくと、0から7のルールインデックスにできますね。
3状態(000,001,010,011,100,101,110,111)は10進数(0,1,2,3,4,5,6,7)となり、
ルール(0, 1, 0, 1, 1, 0, 1, 0)を対応させると
シェルピンスキーのギャスケットの次の行の状態を作れますね。
ちゃんと、
左右同じもので挟まれると０状態、
異なるものに挟まれると１状態という仕組み
になってますね。
＜コード化＞
#初期セルｓ１をはさんで左右に5個の0を並べる
cnt=0
C1={0,0,0,0,0,1,0,0,0,0,0}
C2={0,0,0,0,0,0,0,0,0,0,0}
C3={0,0,0,0,0,0,0,0,0,0,0}
C4={0,0,0,0,0,0,0,0,0,0,0}
C5={0,0,0,0,0,0,0,0,0,0,0}
C6={0,0,0,0,0,0,0,0,0,0,0}
C7={0,0,0,0,0,0,0,0,0,0,0}
C8={0,0,0,0,0,0,0,0,0,0,0}
C9={0,0,0,0,0,0,0,0,0,0,0}
C10={0,0,0,0,0,0,0,0,0,0,0}
C11={0,0,0,0,0,0,0,0,0,0,0}
Cs={C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11}
FormulaText(Cs)
Rule={0,1,0,1,1,0,1,0}
#更新ボタンのイメージはiを2から10まで動かす。
#Cs(i)=Rule(4 Cs(i-1)+2 Cs(i)+Cs(i+1)+1)これをgeogebraの1スタートインデックスで更新する。
#論理的にはもっと簡単にかけるが、
オブジェクトのメモリ配置などの暗黙ルールを推定して順次的コードにする
cnt=cnt+1
c2=Rule(4 Element(Element(Cs,cnt),1)+2 Element(Element(Cs,cnt),2)+Element(Element(Cs,cnt),3)+1)
c3=Rule(4 Element(Element(Cs,cnt),2)+2 Element(Element(Cs,cnt),3)+Element(Element(Cs,cnt),4)+1)
c4=Rule(4 Element(Element(Cs,cnt),3)+2 Element(Element(Cs,cnt),4)+Element(Element(Cs,cnt),5)+1)
c5=Rule(4 Element(Element(Cs,cnt),4)+2 Element(Element(Cs,cnt),5)+Element(Element(Cs,cnt),6)+1)
c6=Rule(4 Element(Element(Cs,cnt),5)+2 Element(Element(Cs,cnt),6)+Element(Element(Cs,cnt),7)+1)
c7=Rule(4 Element(Element(Cs,cnt),6)+2 Element(Element(Cs,cnt),7)+Element(Element(Cs,cnt),8)+1)
c8=Rule(4 Element(Element(Cs,cnt),7)+2 Element(Element(Cs,cnt),8)+Element(Element(Cs,cnt),9)+1)
c9=Rule(4 Element(Element(Cs,cnt),8)+2 Element(Element(Cs,cnt),9)+Element(Element(Cs,cnt),10)+1)
c10=Rule(4 Element(Element(Cs,cnt),9)+2 Element(Element(Cs,cnt),10)+Element(Element(Cs,cnt),11)+1)
Next={Element(Element(Cs,cnt),1),c2,c3,c4,c5,c6,c7,c8,c9,c10,Element(Element(Cs,cnt),11)}
If(cnt==1,SetValue(C2,Next))
If(cnt==2,SetValue(C3,Next))
If(cnt==3,SetValue(C4,Next))
If(cnt==4,SetValue(C5,Next))
If(cnt==5,SetValue(C6,Next))
If(cnt==6,SetValue(C7,Next))
If(cnt==7,SetValue(C8,Next))
If(cnt==8,SetValue(C9,Next))
If(cnt==9,SetValue(C10,Next))
If(cnt==10,SetValue(C11,Next))
#リセットボタン
cnt=0
C1={0,0,0,0,0,1,0,0,0,0,0}
C2={0,0,0,0,0,0,0,0,0,0,0}
C3={0,0,0,0,0,0,0,0,0,0,0}
C4={0,0,0,0,0,0,0,0,0,0,0}
C5={0,0,0,0,0,0,0,0,0,0,0}
C6={0,0,0,0,0,0,0,0,0,0,0}
C7={0,0,0,0,0,0,0,0,0,0,0}
C8={0,0,0,0,0,0,0,0,0,0,0}
C9={0,0,0,0,0,0,0,0,0,0,0}
C10={0,0,0,0,0,0,0,0,0,0,0}
C11={0,0,0,0,0,0,0,0,0,0,0}
Cs={C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11}
＜振り返り＞
次に何になるかはＬとＲが同じが別かで決まるので排他的論理和(ＸＯＲ，⊕）を使ってＬ⊕Ｒで計算することもできそうです。a⊕b はmod(a+b,2)でも代用できます。
(a,b)=(0,0)なら0だし、(a,b)=(1,1)のとき0にするからです。
また、ＬとＲにあたるセルを毎回読み取るのではなく、もとのＣｓの１つ手前を読むリストをＬにし、
１つ先を読むリストをＲにして、先に作っておくと、Ｌ+Rの成分xをMod(x,2)にZipなどで変換すれば、
Nextを数行で作れるかもしれないね。
ただ、スマートでスリムにすればするほど、コードの可読性は落ちます。
だから、わかったあとの進化版として作るにはよいでしょう。
シェルピンスキーのギャスケット（セルオートマトン版）
5×5のライフゲームは以前geogebraで作りましたね。（くわしくはこちら）
これを改良しましょう。
課題２：表示サイズ8×8以上のライフゲームを作ろう。
ライフゲームは2次元セルオートマトンです。
自己の周りの8隣接和Sum8でvが生(1)死(0)が変化します。
基本ルールの例
・v=0状態セルをかこむ8隣接和Sum8が3なら、1状態（再生、誕生）に変更。next_v=1
・v=1状態セルをかこむ8隣接和Sum8が2か3なら。1状態を持続する。next_v=1
・それ以外なら、0状態。next_v=0
5×5=25のサイズでは隣接行列Aを手動で入力してA.v=countを25要素のベクトルとして求めることで、
各セルのSum8とし、それからnext_vを更新していました。
しかし、8×8＝64のサイズで隣接行列Aを手動入力するのも大変です。
ならばと、コードでAを作ってみましょう。
A = Sequence(Sequence(If(i != j && Min(abs(mod(i-1,8) - mod(j-1,8)), 8 - abs(mod(i-1,8) - mod(j-1,8))) <= 1 && Min(abs(floor((i-1)/8) - floor((j-1)/8)), 8 - abs(floor((i-1)/8) - floor((j-1)/8))) <= 1, 1, 0), j, 1, 64), i, 1, 64)
こうすれば、なんとかできますが重いです。
そこで、発想を変えます。
vの1番目の左はｖのラストです。
ｖに対して左隣りはvL=Join({Last(v),Sequence(v(k),k,1,Le-1)})ということです。
ｖを右にスライドしたリストです。テレビゲーム式にはみ出したら反対側から取るとも言えます。
vの右は1番目の次です。
vに対して右隣はvR=Join({Sequence(v(k),k,2,Le),First(v)})ということですね。
vを左にスライドしたリストです。
同じようにして、vの上、下、右上、右下、左上、左下は
行列の辺の長さを考えれば、同じように事前に用意できるはずだね。
すると、Sum8＝左＋右＋上＋下＋右上＋右下＋左上＋左下
というvと同サイズのベクトルが簡単にできます。
また、メモリーをさらに節約するにはSequenceコマンドで切り取るのではなく、
Take(n,m)でビューのように切り取ればいいね。
＜コード化＞
# vは状態ベクトル（セル全部の状態）
# vの表示disp
# 不透明度の設定をいじるのではなく、「v の中身が 1 の場所だけに自動で四角形を描画する（0 の場所は描画しない）」にします。
# vはｒ行ｃ列の行列を1次元のベクトルにしているので、読み取った結果をｖの(r - 1) * 5 + c番目とします。そこが１のときだけ描画します。
# gは世代カウンター
# A.dot(v)をすると、隣接8成分とｖの生死データの内積から、すべてのセルの8隣接和がベクトルcountに反映されます。
N= Slider(5,30,1)
Le = N*N
v = Sequence(RandomBetween(0,1),k,1,Le)
disp=Sequence(Sequence(If(Element(v, (r - 1) * N + c) == 1, Polygon((c, -r), (c + 1, -r), (c + 1, -r - 1), (c, -r - 1))), c, 1, N), r, 1, N)
g = 0
# 「更新」ボタン
SetValue(g, g + 1)
vL = Join({Last(v), Take(v, 1, Le - 1)})
vR = Join({Take(v, 2, Le), {First(v)}})
vD = Join({Take(v, Le - N + 1, Le), Take(v, 1, Le - N)})
vU = Join({Take(v, N + 1, Le), Take(v, 1, N)})
vRU = Join({Take(vR, N + 1, Le), Take(vR, 1, N)})
vRD = Join({Take(vR, Le - N + 1, Le), Take(vR, 1, Le - N)})
vLU = Join({Take(vL, N + 1, Le), Take(vL, 1, N)})
vLD = Join({Take(vL, Le - N + 1, Le), Take(vL, 1, Le - N)})
Sum8 = vL + vR + vD + vU + vRU + vRD + vLU + vLD
SetValue(v, Sequence(If((Element(v, i) == 1 && (Element(Sum8, i) == 2 || Element(Sum8, i) == 3)) || (Element(v, i) == 0 && Element(Sum8, i) == 3), 1, 0), i, 1, Le))
# 「ランダムな初期値」ボタン
SetValue(g, 0)
SetValue(v, Sequence(RandomBetween(0,1),k,1,Le))
#数式に次の行をはりつける。できたテキストラベルを適当に移動してください。
""+g +"世代"
サイズの大きいライフゲーム
過去の有限個の記号発生が次の記号発生に影響する情報源を
マルコフ（Markov)情報源といいい、その状態遷移がマルコフ連鎖だ。
定常確率が可能な情報源をエルゴード性があるという。
くわしくはこちら（くわしくはこちら）
これをさらに探求しよう。
課題３：マルコフ情報源と定常確率（お天気シミュレーション）「天気（晴れ(0)、曇(1)、雨(2)）の確率w0,w1,w2の和が１。
iの翌日がjになる条件付き確率をP(j|i)をpijと書き、
その行列P={pij}が{{12/18,5/18,1/18},{4/7,2/7,1/7},{1/4,1/4,2/4}}のとき、
定常確率w0,w1,w2を求めよう。」
＜解析解＞
定常確率は未知数として、連立方程式
p0=p0p00+p1p10+p2p20=w0*12/18+w1*4/7+w2*1/4=w0
p1=p0p01+p1p11+p2p21=w0* 5/18+w1*2/7+w2*1/4=w1
p2=p0p02+p1p12+p2p22=w0* 1/18+w1*1/7+w2*2/4=w2
w0+w1+w2=1
これを手動で解いてみる。
つまり、
-6/18w0+4/7w1+1/4w2=0...a式
5/18w -5/7w1+1/4w2=0...b式
1/18w0+1/7w1-2/4w2=0...c式
1/4 w0+1/4w1+1/4w2=1/4...d式
w2を消す。
a-b: -11/18w0+9/7w1=0...e式
b-d:(5/18-1/4)w0-(5/7+1/4)w1=-1/4 ...f式
e*7*18:-77w0+9*18w1=0
w1=77/(9*18) w0
これをfに代入。
(5*4-18*1)/18*4 w0-(5*4+1*7)/7*4 (77/(9*18) w0=-1/4
(1/9-27*11/(9*18)w0=-1
(-4+66)w0=36
w0=36/62=0.5806451612
w1=36/62*77/(9*18)=0.27598566
w2=1-(w0+w1)=0.14336918
＜反復法＞
確率w=(w0,w1,w2)を行列Pをかけると、翌日の確率が計算できる。これを反復していけば、確率wは徐々に定常的になるだろう。
geogebraでは、行ベクトルwと行列Pをかけるとき
P・t(w)は
P Transpose(w)で、正しくできます。
しかし、
w・ Pは
w Pでは、ｗのサイズではなくPのサイズの行列が表示されて、
Dot(w,P)では未定義になります。
Sequence(Dot(w, {P(1,k),P(2,k),P(3,k)}),k,1,Length(w))でやっとw・Pを計算して
行ベクトルもどきの
リストを返してくれます。
行列の扱い方があいまいです。リストの内積でリストを返す構造が安全です。
行列と認識されると表示が丸カッコに変わりますが、
中カッコ{}のままはリストとして認識されているということです。
＜コード化＞
メインのメニューの設定で小数点以下10桁にします。
g=0
P={{12/18,5/18,1/18},{4/7,2/7,1/7},{1/4,1/4,2/4}}
w={0.2,0.1,0.7} #適当な初期確率（リスト）
next_w={} #変化する確率(リスト）
"確率推移：晴れw0="+w(1)+"くもりw1="+w(2)+"雨w2="+w(3)+""
#「更新」ボタンに貼り付ける。転置に注意
SetValue(g, g + 1)
next = Sequence(Dot(w, {P(1,k),P(2,k),P(3,k)}),k,1,Length(w))
SetValue(w,next)
#「リセット」ボタンにはりつける。
SetValue(g,0)
SetValue(w, {0.2,0.1,0.7})
27回ほどで解析解とほぼ同じになります。