第一節：問題

2015年的時候有朋友曾經給我分享過一個在知乎上遇到的隨機遊走問題，現在我已經找不到原始鏈接了，但是題目還記得很清楚，問題如下：

假設我們有一個一維坐標軸，一個粒子一開始在坐標原點處。拋一枚均勻的硬幣，如果硬幣向上，粒子向右走一步，如果硬幣向下，粒子向左走兩步。在坐標為 1 處有一個吸收壁，一旦粒子達到吸收壁，那麼它就不能再移動。問，粒子被吸收壁吸收的概率。

這是一個有單側吸收壁的概率隨機遊走問題。這個問題有解析解，也可以用蒙特卡洛模擬來計算數值解。這個問題非常有趣，我將這兩種解法羅列如下。

第二節：精確解

這個問題可以用迭代的方法求解。先將問題拓展為，如果硬幣朝上，那麼粒子向右移動一步；如果硬幣朝下，那麼粒子向左移動 步, . 假設粒子被吸收的概率為 ，那麼 滿足這個方程：

如果 , 那麼方程的解為

注意到， 始終為概率方程的一個解。這是一個平庸解。於是原問題的解就是該粒子被吸收壁吸收的概率為 .

我們知道，當方程次數太高的時候，我們沒有求根公式。所以對於任意的 ，我們不能解得精確解，而只能數值求解。或許對於這個特殊形式的方程，能找到求根公式，我也沒有嘗試。下面用牛頓法求該方程對於任意高次冪情況下的數值解, C++程序如下：

#include <iostream>#include <fstream>#include <cmath>using namespace std;double power(double p, int n){ double result = 1; for (int i = 0; i < n; ++i) result *= p; return result;}double f(double p, int n){ return power(p, n+1) - 2*p + 1;}double derivative(double p, int n){ return (n+1)*power(p, n) - 2;}double newton(int n, double p0){ double p = p0; int max_iteration = 20; int counter = 0; double eps = 1.0e-14; cout.precision(16); while(true) { counter++; if (counter > max_iteration) break; p = p0 - f(p0, n)/derivative(p0, n); double error = fabs(p - p0); cout << "counter = " << counter << ", error = " << error << ", value = " << p << endl; if (error < eps) break; p0 = p; } cout << endl; return p;}int main(int argc, char **argv){ int max_step_length = 10; double p0 = 0.1; ofstream writer("exact.txt"); for (int i = 1; i <= max_step_length; ++i) { writer << i << " " << newton(i, p0) << endl; } writer.close(); return 0;}

第三節：蒙特卡洛模擬結果

除了上面的精確解，我們也可以用蒙特卡洛模擬求解。實驗中我們有一個天然的隨機數生成器，就是材質均勻的硬幣。蒙特卡洛模擬過程中，我們用C++自帶的隨機數生成器代替硬幣，然後按照實驗的要求讓粒子隨機遊走。重複多次，可以算出粒子被吸收的頻率，當實驗次數足夠多的時候，按照大數定理，我們就可以得到粒子被吸收的概率。C++程序如下：

#include <iostream>#include <fstream>#include <cstdlib>#include <vector>#include <cassert>#include <cmath>#include <utility>using namespace std;double uniform(){ return rand()/float(RAND_MAX);}double uniform(double a, double b){ return a + (b-a)*uniform();}class MonteCarlo{ private: const int right_step_length = 1; int left_step_length = 1; const int no_return_point = -100; const int steps_limit = 1000; public: MonteCarlo(){} MonteCarlo(int left_step_length) { this->left_step_length = left_step_length; } void set_left_step_length(int left_step_length) { this->left_step_length = left_step_length; } bool simulation() { int current_position = 0; int counter = 0; while(current_position > no_return_point) { counter++; if (counter > steps_limit) { return false; } double random_value = uniform(); if (random_value < 0.5) { current_position -= left_step_length; } else { current_position += right_step_length; } if (current_position == 1) return true; } return false; } double simulation(int simulation_times) { assert(simulation_times >= 1); int success_counter = 0; for (int i = 0; i < simulation_times; ++i) { if (simulation()) success_counter += 1; else continue; } return float(success_counter)/simulation_times; }};vector<double> monte_carlo_simulation(int test_times, int left_step_length){ vector<double> results; MonteCarlo *monteCarlo = new MonteCarlo(left_step_length); int simulation_times = 10000; for (int i = 0; i < test_times; ++i) { results.push_back(monteCarlo->simulation(simulation_times)); } delete monteCarlo; return results;}template <typename T>void print_vector(const vector<T> &my_vector){ for (int i = 0; i < my_vector.size(); ++i) { cout << my_vector[i] << endl; }}template <typename T>double mean(const vector<T> &my_vector){ int length = my_vector.size(); assert(length >= 1); double s = 0; for (int i = 0; i < length; ++i) { s += my_vector[i]; } return s/float(length);}template <typename T>double standard_deviation(const vector<T> &my_vector){ int length = my_vector.size(); assert(length > 1); double mu = mean(my_vector); double s = 0; for (int i = 0; i < length; ++i) { s += (my_vector[i] - mu)*(my_vector[i] - mu); } return sqrt(s);}void random_walk(){ int scan_number = 10; vector<int> left_step_lengths; for (int i = 0; i < scan_number; ++i) { left_step_lengths.push_back(i+1); } vector< pair<double, double> > results; MonteCarlo *monteCarlo = new MonteCarlo(); int test_times = 20; for (int i = 0; i < scan_number; ++i) { cout << "i = " << i+1 << ", total = " << scan_number << endl; int left_step_length = left_step_lengths[i]; monteCarlo->set_left_step_length(left_step_length); vector<double> simulation_result = monte_carlo_simulation(test_times, left_step_length); pair<double, double> result(mean(simulation_result), standard_deviation(simulation_result)); results.push_back(result); } delete monteCarlo; ofstream writer("monte_carlo_results.txt"); for (int i = 0; i < results.size(); ++i) { writer << left_step_lengths[i] << " " << results[i].first << " " << results[i].second << endl; } writer.close();}int main(int argc, char **argv){ random_walk(); return 0;}

第四節：對比結果

前面給出了原問題的兩種解法，這兩種解法可以給出精度範圍內一致的結果。如下圖所示：

從圖中可以看出，當向左的步長增大的時候，粒子被吸收壁吸收的概率的顯著減小，而且最終趨近於0.5. 這個結果很容易理解，因為當向左的步長足夠大的時候，如果第一次拋硬幣結果朝下，那麼一旦粒子向左走了一次，那麼它能重新返回坐標原點的概率幾乎是趨近於等零了。所以當向左的步長 足夠大的時候，粒子可以被吸收的概率幾乎上就是第一次拋硬幣結果向上的概率，也就是0.5.

注意到當向左的步長為一的時候，這個粒子總是可以被吸收壁吸收，這是因為在布朗運動中有一個經典的結果：當空間維數小於等於二的時候，做布朗運動的粒子是可以返回初始點的，而空間維數大於二的時候，做布朗運動的粒子卻回不到原點。用一句通俗的話解釋就是，一個醉漢總是可以回家，但是一隻喝醉的小鳥卻回不了家了。這就是 Polyas Recurrence Theorem.