如何計算空方格數的期望？

在一個的網格上有900隻跳蚤，每隻跳蚤一個方格。

從某一時刻開始，每個單位時間，所有的跳蚤隨機地跳到相鄰的方格（一般情況有四種可能，除了邊上的跳蚤）。

求50個單位時間過後，空方格數的期望？

此為Project Euler的213題：https://projecteuler.net/problem=213

---

330.721154

格子i在50步以後是空的。

是每個跳蚤50步後不在i的概率之積。

代碼很醜僅供跑結果， 實際上由於對稱性很多是不用算的，懶得改。

#include &
#include &
#include &
using namespace std;

double board[2][30][30];

void run(int c){
memset(board[!c], 0, sizeof(double)*30*30);
for(int i=0;i&<30;++i){ for(int j=0;j&<30;++j){ int neighbor = 4; double v = board[c][i][j]; if(i==0 || i==29){ --neighbor; } if(j==0 || j==29){ --neighbor; } if(i&>0){
board[!c][i-1][j] += v/neighbor;
}
if(j&>0){
board[!c][i][j-1] += v/neighbor;
}
if(i&<29){ board[!c][i+1][j] += v/neighbor; } if(j&<29){ board[!c][i][j+1] += v/neighbor; } } } } double pnotstay[30][30]; int main(){ for(int i=0;i&<30;++i){ for(int j=0;j&<30;++j){ pnotstay[i][j] = 1; } } for(int i=0;i&<30;++i){ for(int j=0;j&<30;++j){ memset(board, 0, sizeof(double)*900*2); board[0][i][j] = 1; for(int c=0;c&<50;++c){ run(c%2); } for(int x=0;x&<30;++x){ for(int y=0;y&<30;++y){ pnotstay[x][y] *= 1-board[0][x][y]; } } } } double res=0; for(int i=0;i&<30;++i){ for(int j=0;j&<30;++j){ res += pnotstay[i][j]; } } printf("%.6lf ", res); }

---

空方格數的期望等於每個方格是空的概率之和。

每個格子是空的概率可以每個跳蚤50步之後的概率分布得出。

完畢

---

湊個熱鬧，來個python版的，耗時130毫秒：

In [1]: def px2y(x, y, n):
...: if abs(x-y) != n and (x//n != y//n or abs(x-y) != 1):
...: return 0.0
...: elif x in [0, n-1, n*(n-1), n*n-1]:
...: return 1/2
...: elif x &< n or x &> n*(n-1) or x%n in [0, n-1]:
...: return 1/3
...: else:
...: return 1/4
...:
...: def solve(n, m):
...: tran = np.array([[px2y(x, y, n) for y in range(n*n)] for x in range(n*n)])
...: return (1 - np.linalg.matrix_power(tran, m)).prod(axis=0).sum()
...:

In [2]: %time print(solve(30, 50))
330.721154014
Wall time: 133 ms


---

其實這是道 水題,答案是900/e...

--------------------------------正經的分界線---------------------------------------

這是一個Markov Chain Process

空方格數的期望等於每個方格是空的概率之和。

這個轉移矩陣是個對角矩陣,但是邏輯上比較繞.

(*矩陣大小n[Times]n,時間time,精度acc*)
move[n_, time_, acc_: 9] := Module[{L}, L = Table[0, {n^2}, {n^2}];
L[[2, 1]] = 1/2; L[[n + 1, 1]] = 1/2; L[[n - 1, n]] = 1/2;
L[[2 n, n]] = 1/2;
L[[n^2 - 2 n + 1, n^2 - n + 1]] = 1/2;
L[[n^2 - n + 2, n^2 - n + 1]] = 1/2; L[[n^2 - n, n^2]] = 1/2;
L[[n^2 - 1, n^2]] = 1/2;
Do[L[[i - 1, i]] = L[[i + 1, i]] = L[[i + n, i]] = 1/3, {i, 2,
n - 1}];
Do[L[[i - n, i]] = L[[i + 1, i]] = L[[i + n, i]] = 1/3, {i, n + 1,
n^2 - 2 n + 1, n}];
Do[L[[i - n, i]] = L[[i - 1, i]] = L[[i + n, i]] = 1/3, {i, 2 n,
n^2 - n, n}];
Do[L[[i - n, i]] = L[[i - 1, i]] = L[[i + 1, i]] = 1/3, {i,
n^2 - n + 2, n^2 - 1}];
Do[L[[n (i - 2) + j, n (i - 1) + j]] =
L[[n (i - 1) + j - 1, n (i - 1) + j]] =
L[[n i + j, n (i - 1) + j]] =
L[[n (i - 1) + j + 1, n (i - 1) + j]] = 1/4, {i, 2, n - 1}, {j,
2, n - 1}];
SetPrecision[Total[Times @@@ (1 - MatrixPower[L // N, time])], acc]]

(*代入題設*)
move[30, 50] // Timing


輸出很快

{0.296875, 330.721154}

哈哈,看來Mathematica矩陣運算也是不輸於MatLab的.

具體理論可以看這裡markov process

實際計算過程就一行,主要是這個矩陣L很不好說清楚.

Total[Times @@@ (1 - MatrixPower[L // N, time])]


---

推薦閱讀：