目錄

• 動手學隨機模擬
• • 1. 如何通過隨機模擬估計看漲期權的報酬分布
  • • 1.1 金融知識：股票與看漲期權
    • 1.2 問題分析與確定建模思路
    • 1.3 如何模擬股票價格：布朗運動
  • 1.4 如何更真實地模擬股票價格：幾何布朗運動
  • • 1.5 估計看漲期權的報酬分布
    • 1.6 隨機模擬的幾個關鍵點
  • 2.隨機模擬的關鍵：隨機數生成
  • • 2.1 隨機數生成的基礎：均勻分布的隨機數
    • • 構造滿足某個離散分布的隨機數
    • 2.2 分布函數生成隨機數(逆變換)
    • 2.3 拒絕采樣生成隨機數
    • 2.4 MCMC（蒙特卡洛馬爾可夫）生成隨機數

動手學隨機模擬

1. 如何通過隨機模擬估計看漲期權的報酬分布

關鍵字:隨機模擬，看漲期權，報酬分布

1.1 金融知識：股票與看漲期權

• 權利金
• 行權價格
• 行權
• 買進看漲期權
• 期權的標的資產
  期權買方享有權利而不用承擔義務

import numpy as np import pandas as pd import matplotlib.pyplot as plt %matplotlib inline plt.rcParams['font.sans-serif']=['SimHei','Songti SC','STFangsong'] plt.rcParams['axes.unicode_minus'] = False # 用來正常顯示負號 price_arr = np.linspace(0, 40, 10000) # 定義價格波動，從0～100波動 def get_earnings(price: float, contract_sum: float, contract_cost: float) -> float:return price - (contract_sum + contract_cost) if price > contract_sum else -contract_costmy_earn = [get_earnings(price=price, contract_sum=20, contract_cost = 2) for price in price_arr] plt.figure(figsize=(8,6)) plt.plot(price_arr, my_earn, lw=4, alpha=0.75) plt.xlabel("Price") plt.ylabel("my_earns") plt.title("榴蓮合約的收益曲線") plt.show()

購買了看漲期權，收益上不封頂，虧損最多就是合約金（權利金）
一般來說，權利金的數額是遠遠低于標的資產（股票）本身的資產價格。換句話來說，期權是有杠桿的，現在你用2元的本金撬動了幾十元的收益，這是個十分關鍵的地方

1.2 問題分析與確定建模思路

看漲期權能得到報酬的關鍵點就是標的資產的價格，對于看漲期權來說，標的資產的價格越高，我們獲得報酬相應的也會越高。

問題的關鍵是:

• （1）在某個時刻如何估計股票的價格分布
• （2）在某個時刻結合行權價格與股票的價格分布，模擬看漲期權的報酬分布

股票價格是一個隨機過程，而看漲期權的價格（也是衍生品的價格）是隨機過程的函數
估計股票價格=>使用隨機過程近似股票價格這個隨機過程

# 獲取工商銀行股票價格數據 import baostock as bs #登錄系統 lg=bs.login() #顯示登錄返回信息 print('login respond error_code:'+lg.error_code) print('login respond error_msg:'+lg.error_msg) stock_code = 'sh.601398' rs = bs.query_history_k_data_plus(stock_code,"date,close",start_date='2016-01-04', end_date='2021-12-31',frequency="d", adjustflag="3") print('query_history_k_data_plus respond error_code:'+rs.error_code) print('query_history_k_data_plus respond error_msg:'+rs.error_msg) #### 打印結果集 #### data_list = [] while (rs.error_code == '0') & rs.next():# 獲取一條記錄，將記錄合并在一起data_list.append(rs.get_row_data())gsyh = pd.DataFrame(data_list, columns=rs.fields) gsyh['close'] = gsyh['close'].astype('float') gsyh['date'] = pd.to_datetime(gsyh['date'], format='%Y-%m-%d') gsyh # 查看工商銀行的價格走勢圖 plt.figure(figsize=(8,6)) plt.plot(gsyh['date'], gsyh['close'], lw=2, alpha=0.8) plt.xlabel("date") plt.ylabel("Price") plt.title("工商銀行價格走勢圖") plt.show()

我們發現曲線有如下特點：

• （1）股票價格走勢具有強烈的隨機性；
• （2）股票價格的曲線是一條"不平滑"的曲線，與我們前面學習的可微分的曲線不太一樣。

可以嘗試使用物理中的布朗運動去近似股票價格的無規律波動,數學上對布朗運動的描述發展相比于物理上的布朗運動要慢一些

1.3 如何模擬股票價格：布朗運動

設有一個粒子在直線上隨機游動, 該粒子在每個單位時間內等可能地向左或向右移動一個單位長度，觀察這個例子的軌跡圖

# 獲取隨機游動的位置 def get_particle_pos(n_times):position_ls = [0]for i in range(n_times-1):left_or_right = np.random.rand()if left_or_right >= 0.5:position_ls.append(1.0)else:position_ls.append(-1.0)position_ls = np.cumsum(position_ls)return position_lsn_times = 10000 position_arr = get_particle_pos(n_times = n_times) plt.figure(figsize=(8, 6)) plt.plot(np.arange(n_times), position_arr, lw=2, alpha=0.8) plt.xlabel("step") plt.ylabel("position") plt.title("隨機游動的位置") plt.show()

我們將粒子的隨機游動與股價走勢放在一起對比一下：

plt.figure(figsize=(16, 6)) plt.subplot(1,2,1) plt.plot(np.arange(n_times), position_arr, lw=2, alpha=0.8) plt.xlabel("step") plt.ylabel("position") plt.title("隨機游動的位置") plt.subplot(1,2,2) plt.plot(gsyh['date'], gsyh['close'], lw=2, alpha=0.8) plt.xlabel("date") plt.ylabel("Price") plt.title("工商銀行價格走勢圖") plt.show()

通過對比一維布朗運動走勢和股票價格的走勢可以發現，這兩者在形狀上非常相似

布朗運動詳細定義
隨機過程 $\{X(t),t?0\}$ 如果滿足:
(1) $X(0)=0$;
(2) $\{X(t),t?0\}$ 有平穩獨立增量;
(3) 對每個 $t>0,X(t)$ 服從正態分布 $N\left(0,{\sigma }^{2}t\right)$,
則稱 $\{X(t),t?0\}$ 為 Brown 運動, 也稱為 Wiener 過程, 常記為 $\{B(t),t?0\}$ 或 $\{W(t),t?0\}$ 。
如果 $\sigma =1$, 我們稱為標準 Brown 運動; 如果 $\sigma \ne 1$, 則可考慮 $\{X(t)/\sigma ,t?0\}$, 它是標準 Brown 運動。

我們可以驗證一維布朗運動（隨機游走）的定義（3），即：均值函數和方差函數分別為：0與${\sigma }^{2}t$以及每個時刻的位置分布為正態分布。

plt.figure(figsize=(16, 6)) plt.subplot(1,2,1) ## 模擬10000次隨機漫步 simulate_time = 10000 n_times = 100 result_ls = [] for i in range(simulate_time):position_arr = get_particle_pos(n_times=n_times)plt.plot(np.arange(n_times), position_arr, lw=2, alpha=0.8)result_ls.append(position_arr.tolist()) pos_mean = np.mean(result_ls, axis=0) pos_var = np.var(result_ls, axis=0) plt.plot(np.arange(n_times), pos_mean, lw=4, alpha=0.8, c='red', label='mean_func') plt.plot(np.arange(n_times), pos_var, lw=4, alpha=0.8, c='blue', label='var_func') plt.xlabel("step") plt.ylabel("position") plt.title("隨機游動的位置") plt.legend() plt.subplot(1,2,2) ## 選取index=60的位置分布 pos_60 = np.array(result_ls)[:, 60] plt.hist(pos_60, bins=100, density=1) plt.title("index=60的位置分布") plt.show()

隨機游動的均值函數為0， 方差函數是一個關于時間$t$的直線函數，index為60的位置分布也符合正態分布，以上都是布朗運動的表現。

plt.figure(figsize=(12, 8)) simulate_times = 12 n_times = 10000 for i in range(simulate_times):position_arr = get_particle_pos(n_times=n_times)plt.plot(np.arange(n_times), position_arr, lw=2, alpha=0.75) plt.axhline(y=0, lw=6) plt.plot(np.arange(n_times), np.sqrt(np.arange(n_times)), lw=6, c='red', alpha=0.8, label=r'$t=pos^2$') plt.plot(np.arange(n_times), -np.sqrt(np.arange(n_times)), lw=6, c='red', alpha=0.8) plt.legend() plt.xlabel("t") plt.ylabel("position") plt.show()

• （1）布朗運動的軌跡線會頻繁穿越$pos=0$的水平線
  我們把0基準線看成是股票的開盤價，頻繁穿越0基準線意味著股價很大概率會在開盤價上下波動，而不是始終維持在開盤價的上方或者下方。這一點的性質在很多股票投資的實踐經驗都有印證
• （2）在任意時刻 $t$, 樣本軌跡的位置 $B(t)$ 不會偏離正負一個標準差太遠：我們在布朗運動的定義中的（3）可以知道，對每個 $t>0,X(t)$ 服從正態分布 $N\left(0,{\sigma }^{2}t\right)$，一個標準差即$\sigma \sqrt{t}$，上圖中的紅色曲線可以包含絕大多數的樣本軌跡。這點可以說明，隨著交易時間$t$的推移，$t$時刻的股票價格始終不會偏離開盤價（0基準線）$\pm \sqrt{t}\times$ 價格波動的標準差，就算偏離也不會太遠（除了小概率事件）。
  缺點
  布朗運動可以是負數，因為布朗運動的樣本軌跡會頻繁穿越0基準線。然而，股票價格卻肯定不會是負數，這是布朗運動的理論與現實的矛盾。為了解決這個問題，我們可以改進布朗運動描述股票價格走勢，使用幾何布朗運動描述股票價格。

1.4 如何更真實地模擬股票價格：幾何布朗運動

為了解決使用布朗運動描述股票價格走勢會出現負數的問題，我們嘗試給布朗運動加上一個僅和時間 $t$ 有關的漂移項 $\mu t$, 以及一個尺度參數 $\sigma$，這樣就可以得到帶漂移項的布朗運動：

因為 $X(t)$ 與 $B(t)$ 的取值隨著時間 $t$ 的變化也可以是負數, 但是股票的價格顯然不能是負數，這種辦法并沒有解決實際問題。

我們假設使用$S(t)$表示股票價格而$X(t)$是股票價格的自然對數，即：$X(t)=ln(S(t))$。因此，帶漂移項的布朗運動就變成了：
$$X(t)=ln(S(t))=\mu t+\sigma B(t)$$

我們對$X(t)=ln(S(t))=\mu t+\sigma B(t)$進行微分(簡單的對數函數求導)，即：
$$\frac{dS(t)}{S(t)}=\mu dt+\sigma dB(t)$$
$dS(t)/S(t)$ 就 是這段間隔內的收益率
化簡上述的式子$\frac{dS(t)}{S(t)}=\mu dt+\sigma dB(t)$得到：
$$dS(t)=\mu S(t)dt+\sigma S(t)dB(t)$$

式子：$dS(t)=\mu S(t)dt+\sigma S(t)dB(t)$是一個微分方程，這個微分方程是帶有隨機性的微分方程，簡稱"隨機微分方程"。因此，滿足上述隨機微分方程$dS(t)=\mu S(t)dt+\sigma S(t)dB(t)$的股價 $S(t)$ 是一個幾何布朗運動。

我們可以使用歐拉離散化可以得到離散時間模型：
$$S_t=S_{t?\Delta t}{e}^{\left(\mu ?\frac{1}{2}{\sigma }^2\right)\Delta t+\sigma {\epsilon }_t\sqrt{\Delta t}}$$
離散時間模型的參數含義：

• $S_t$ : 證券在t時刻的價格;
• ${\mathrm{S}}_{\mathrm{t}?\Delta \mathrm{t}}$ : 證券在t $?\Delta \mathrm{t}$ 時刻的價格；
• $\mu$ ：證券收益率的期望值;
• $\sigma$ ：證券收益率的波動率;
• ${\epsilon }_t$ : 服從正態分布（期望為 0 , 方差為 1 ）

接下來，我們模擬證券初始價格為10（日收益率均值漂移項為0.005，波動率為0.25），模擬天數為300天，次數為5000次的幾何布朗運動價格。

# 我們模擬證券初始價格為10（日收益率均值漂移項為0.005，波動率為0.25），模擬天數為300天，次數為5000次的幾何布朗運動價格 S0 = 10 # 初始價格 n_days = 300 # 模擬天數 mu = 0.005 sigma = 0.25 delta_t = 1/n_days n_times = 5000 # 模擬次數 price_arr = np.zeros((n_times, n_days)) # n_times*n_days個價格 price_arr[: , 0] = S0 # 初始價格 for t in range(1, n_days):price_arr[:, t] = price_arr[:, t-1]*np.exp((mu-0.5*np.square(sigma))*delta_t+sigma*np.random.standard_normal(n_times)*np.sqrt(delta_t)) # 分別繪制一條軌跡和所有軌跡 plt.figure(figsize=(16, 6)) plt.subplot(1,2,1) plt.plot(np.arange(n_days), price_arr[0, :], lw=2, alpha=0.75) plt.xlabel("t") plt.ylabel("Price") plt.title("1條股價樣本軌跡") plt.subplot(1,2,2) for i in range(n_times):plt.plot(np.arange(n_days), price_arr[i, :], lw=2, alpha=0.75) plt.xlabel("t") plt.ylabel("Price") plt.title("5000條股價樣本軌跡") plt.show()

# 繪制某一天的價格分布 plt.figure(figsize=(8, 6)) plt.hist(price_arr[60, :], bins=50, density=True) plt.xlabel("price") plt.title("index=60的價格分布") plt.show()

以上的案例給定了離散時間模型的參數值，這是一種只會存在于習題中的情況，因為當我們想模擬一只具體的股票時，這些參數信息不會一下子彈到我們腦邊，我們需要根據實際數據估計其中的參數值。估計方法如下（不給證明了）：
$$\begin{array}{rl}\mu & =E\left[\ln \left(\frac{S_{t+\Delta t}}{S_t}\right)\right]/\Delta t\\ {\sigma }^2& =\mathrm{var}\left[\ln \left(\frac{S_{t+\Delta t}}{S_t}\right)\right]/\Delta t\end{array}$$

# 使用工商銀行股票價格數據和幾何布朗運動模擬股價 ## 1. 獲取2021年的工商銀行股票價格數據 # 登陸系統 lg = bs.login() # 顯示登陸返回信息 print('login respond error_code:'+lg.error_code) print('login respond error_msg:'+lg.error_msg) stock_code = 'sh.601398' rs = bs.query_history_k_data_plus(stock_code,"date,close",start_date='2021-01-01', end_date='2021-12-31',frequency="d", adjustflag="3") print('query_history_k_data_plus respond error_code:'+rs.error_code) print('query_history_k_data_plus respond error_msg:'+rs.error_msg) #### 打印結果集 #### data_list = [] while (rs.error_code == '0') & rs.next():# 獲取一條記錄，將記錄合并在一起data_list.append(rs.get_row_data()) gsyh_2021 = pd.DataFrame(data_list, columns=rs.fields) gsyh_2021['close'] = gsyh_2021['close'].astype('float') gsyh_2021['date'] = pd.to_datetime(gsyh_2021['date'], format='%Y-%m-%d')## 2.估計幾何布朗運動中的mu和sigma price_ls = gsyh_2021['close'].tolist() delta_t = 1/len(price_ls) mu = np.mean(np.log(price_ls[1:]) - np.log(price_ls[:-1])) / delta_t sigma2 = np.var(np.log(price_ls[1:]) - np.log(price_ls[:-1])) / delta_t ## 3.使用幾何布朗運動模擬股價走勢 print("==============開始模擬================") S0 = price_ls[-1] # 初始價格 print("初始價格：", S0) n_days = 300 # 模擬天數 sigma = np.sqrt(sigma2) print("估計的mu = ", mu) print("估計的sigma = ", sigma) delta_t = 1/n_days n_times = 10000 # 模擬次數 price_arr = np.zeros((n_times, n_days)) # n_times*n_days個價格 price_arr[: , 0] = S0 # 初始價格 for t in range(1, n_days):price_arr[:, t] = price_arr[:, t-1]*np.exp((mu-0.5*np.square(sigma))*delta_t+sigma*np.random.standard_normal(n_times)*np.sqrt(delta_t)) # 分別繪制一條軌跡和所有軌跡 plt.figure(figsize=(16, 6)) plt.subplot(1,2,1) plt.plot(np.arange(n_days), price_arr[0, :], lw=2, alpha=0.75) plt.xlabel("t") plt.ylabel("Price") plt.title("1條股價樣本軌跡") plt.subplot(1,2,2) for i in range(n_times):plt.plot(np.arange(n_days), price_arr[i, :], lw=2, alpha=0.75) plt.xlabel("t") plt.ylabel("Price") plt.title("10000條股價樣本軌跡") plt.show()

可以從上述模擬中得到第300天的模擬價格分布，即：

import seaborn as sns palette = plt.get_cmap('tab20c') price_300 = price_arr[299, :] plt.figure(figsize=(8,6)) sns.histplot(price_300, kde = True, stat='probability', bins = 20, shrink = 1, color = palette.colors[0], edgecolor = palette.colors[-1]) plt.xlabel("Price") plt.show()

小節：我們通過帶漂移項的布朗運動$X(t)=\mu t+\sigma B(t)$可以知道，$X(t)$ 的期望為 $E[X(t)]=\mu t$，即t時刻的價格期望為$\mu t$。這是個十分有趣的結論，因為假如$\mu >0$，那么t時刻的價格期望也會大于0。因此如果我們堅信股市長期來看是一個向上的牛市（雖然很緩慢），那么我們就應該接受短期波動而堅持長期持股，就像巴菲特的價值投資一樣。

1.5 估計看漲期權的報酬分布

現在我們假設行權價格是4.7元，我們來分析報酬在第300天的時候的報酬分布情況（忽略權利金、遠期利率和股息）：
當股票價格大于行權價格，那么我將獲利；
當股票價格小于行權價格，那么我將虧損；

# 繪制報酬分布圖 price_300 = price_arr[299, :].tolist() contract_sum = 4.7 profit_sum = [price - contract_sum if price > contract_sum else 0 for price in price_300] #剔除負數 plt.figure(figsize=(8,6)) plt.hist(profit_sum, bins=40, density=True) plt.xlabel("Profit") plt.show()

1.6 隨機模擬的幾個關鍵點

我們可以將隨機模擬分解為以下步驟：

• （1）確定事件的過程（仿真）；
• （2）在事件中添加合適的分布的隨機數（隨機模擬）；（關鍵）
• （3）得到模擬的結論；
  采樣算法：（拒絕采樣、MCMC采樣）
  （因為要從某個分布中采樣，采樣的樣本就是所謂的隨機數）

2.隨機模擬的關鍵：隨機數生成

• 計算機生成隨機數的基礎算法：均勻分布采樣
• 分布函數采樣
• 拒絕采樣
• MCMC采樣

2.1 隨機數生成的基礎：均勻分布的隨機數

偽隨機數生成的方法有很多, 比較簡單的一種方式比如:
$${\mathrm{x}}_{\mathrm{n}+1}=(\mathrm{a}{\mathrm{x}}_{\mathrm{n}}+\mathrm{c})\mathrm{mod}\mathrm{m}$$

其中, $X$ 為偽隨機序列,

• $m,m>0$, 模數, 也是生成序列的最大周期
• $a,0<a<m$, 乘數
• $c,0\le c<m$, 增量
• $X_0,0\le X_0<m$, 種子點 (起始值)

${\mathrm{X}}_{\mathrm{n}}/n$ 近似服從 $(0,1)$ 上的均勻分布。

這種方法叫線性同余法，使用線性同余法時需要特別小心參數設置，否則生成的隨機數將會十分糟糕。下面通過一組實驗看看不同的參數設置下生成的隨機數序列：

# 線性同余法 def get_uniform_sample(a, c, m, x0, n_times):uniform_ls = []for i in range(n_times):uniform_sample = (a*x0+c) % muniform_ls.append(uniform_sample)x0 = uniform_samplereturn uniform_lsn_times = 10 # 實驗次數 ## 第一組實驗： a, c, m, x0 = 11, 0, 8, 1 uniform_ls = get_uniform_sample(a, c, m, x0, n_times) print("第一次實驗：") print("a = %d, c = %d, m = %d, x0 = %d:\n" % (a, c, m, x0)) print(uniform_ls) print("============================") ## 第二組實驗： a, c, m, x0 = 25214903917, 11, 2**48, 1 # 這也是java.util.Random中的默認參數設置 uniform_ls = get_uniform_sample(a, c, m, x0, n_times) print("第二次實驗：") print("a = %d, c = %d, m = %d, x0 = %d:\n" % (a, c, m, x0)) print(uniform_ls)

構造滿足某個離散分布的隨機數

上面的線性同余法是構造均勻分布的隨機數,如果我們需要構造滿足某個離散分布的隨機數，如：$x=1,2,3,4$分別滿足概率$[0.1,0.3,0.5,0.1]$這個離散分布的隨機數。其實這個問題的解決方法并不難，需要用到均勻分布的隨機數：

• 生成一個$[0,1]$的均勻分布的隨機數$x$；
• 如果$x$在$[0,0.1)$，那么$x=1$， 如果$x$在$[0.1,0.4)$，那么$x=2$，如果$x$在$[0.4,0.9)$，那么$x=3$，如果$x$在$[0.9,1)$，那么$x=4$。
• 重復以上步驟n_times次就能獲得n_times個符合離散分布的隨機數

# 生成10000個離散分布的隨機數 from collections import Counter def get_discrete_sample(x_ls, p_ls, size):""":param x_ls: 類別序列:param p_ls: 類別的概率序列:param size: 生成隨機數個數:return: 隨機數序列"""p_ls = np.cumsum(p_ls) # [0.1, 0.4, 0.9, 1]x_sample_ls = []for i in range(size):p_uniform = np.random.uniform(0, 1)tmp_ls = [1 if p_uniform < p else 0 for p in p_ls]index = len(tmp_ls) - np.sum(tmp_ls)x_sample_ls.append(x_ls[index])return x_sample_lssize = 10000 x_ls = [1, 2, 3, 4] p_ls = [0.1, 0.3, 0.5, 0.1] x_sample_ls = get_discrete_sample(x_ls, p_ls, size) x_count = Counter(x_sample_ls) x_count = list(x_count.items()) x_ratio = [(x[0], x[1] / len(x_sample_ls)) for x in x_count] print(x_ratio)

2.2 分布函數生成隨機數(逆變換)

通過使用均勻分布的隨機數構造別的分布的隨機數呢？答案是肯定的，這個方法叫逆變換法。

逆變換法的理論依據是：假設 $\mathrm{X}$ 是一個連續型隨機變量，它的累積分布函數為 $F_X(x)$ ； 現在，定義隨機變量 $Y=F_X(X)$ 。那么， $Y$ 服從 $(0,1)$ 上的均勻分布，即
$$P(Y\le y)=y,0<y<1\text{.?}$$

換句話說，任意分布的分布函數值服從$(0,1)$上的均勻分布
接下來，我們嘗試使用逆變換法通過均勻分布的隨機數構造指數分布的隨機數：

• 指數分布的概率密度函數為 $f_X(x)=$ $\lambda e^{?\lambda x},x>0$ ，它的累積分布函數為
  $${\mathrm{F}}_X(x)=1?e^{?\lambda x},\mathrm{x}\ge 0$$

• 產生 $(0,1)$ 均勻分布的隨機數 $V$（當作是指數分布的分布函數值）

• 令 $Y=?\frac{1}{\lambda }\log (V)$（將滿足均勻分布的分布函數進行逆變換），$Y$ 即為服從指數分布的隨機數

逆變換法構造指數分布的隨機數 def get_exp_sample(lmd, size):uniform_samples = np.random.uniform(0, 1, size)exp_samples = []for uniform_sample in uniform_samples:exp_sample = -1.0/lmd*np.log(uniform_sample)exp_samples.append(exp_sample)return exp_sampleslmd = 1 size = 10000 exp_samples = get_exp_sample(lmd, size) plt.figure(figsize=(8, 6)) plt.hist(exp_samples, bins=50, density=True, alpha=0.75) plt.xlabel("X") plt.title("X~exp(1)") plt.show()

2.3 拒絕采樣生成隨機數

逆變換法需要求分布函數和逆變換
拒絕采樣開始采用蒙特卡洛的思想構造隨機數，我們先來看看下圖：

• $p(x)$是我們需要采樣的分布，有時候十分復雜，這里的$p(x)=0.3e^{?(x?0.3{)}^2}+0.7e^{?\frac{(x?2.0{)}^2}{0.3}}$，其中$0\le x\le 4$
• $q(x)$是我們的構造的分布，這里我們假設為[0, 4]的均勻分布

# 拒絕采樣演示： M = 4 z0 = 2.5 x = np.linspace(0,4,1000) q_x = np.array([1.0/4]*1000) p_x = 0.3*np.exp(-(x-0.3)**2)+0.7*np.exp(-(x-2.0)**2/0.3) plt.figure(figsize=(8,6)) plt.plot(x, p_x, lw=4, alpha=0.8, label=r'$p(x) = 0.3 e^{-(x - 0.3)^2} + 0.7e^{-\frac{(x - 2.0)^2}{0.3}}$') plt.plot(x, M*q_x, lw=4, alpha=0.8, label=r'$Mq(x)$') plt.vlines(x=0, ymin=0, ymax=1, colors='gray',linestyles='--') plt.vlines(x=4, ymin=0, ymax=1, colors='gray',linestyles='--') plt.vlines(x=z0, ymin=0, ymax=1, colors='red',linestyles='--', lw=2) plt.scatter(z0, 0.3*np.exp(-(z0-0.3)**2)+0.7*np.exp(-(z0-2.0)**2/0.3), s=90, color='red') plt.text(z0+0.1, 0.3*np.exp(-(z0-0.3)**2)+0.7*np.exp(-(z0-2.0)**2/0.3), r'$u_0$') plt.legend() plt.xlabel("x") plt.ylabel("f(x)") plt.title("拒絕采樣演示圖") plt.show()

由于$p(x)$太復雜，我們無法直接采樣，因此我們借助輔助的分布$q(x)$，這個輔助的分布滿足一個特點：$q(x)$乘一個常數$M$后的$Mq(x)$曲線能完全包住我們想要采樣的分布$p(x)$。

拒絕采樣的過程如下：

• 通過輔助分布$q(x)$采樣$z_0$，輔助分布可以選擇常見的均勻分布或者正態分布；
• 計算$u_0=p(z_0)$（圖中紅色的點）；
• 在區間$(0,Mq(z_0)]$范圍的均勻分布隨機采樣$k$，即在紅色的直線范圍按均勻分布采樣$k$；
• 判斷$k$的大小，
  如果$k<u_0$（即$k$位于$u_0$下方）,則接受$z_0$是來自分布$p(x)$的隨機數；
  反之，如果$k\ge u_0$（即$k$位于$u_0$上方）,則拒絕$z_0$是來自分布$p(x)$的隨機數。
• 重復以上過程，知道采樣個數為size為止。

# 通過拒絕采樣0.3*np.exp(-(x-0.3)**2)+0.7*np.exp(-(x-2.0)**2/0.3)分布的隨機數 from scipy.stats import uniform def get_px(x):return 0.3*np.exp(-(x-0.3)**2)+0.7*np.exp(-(x-2.0)**2/0.3)def get_px_samples(size):M = 4px_sample_ls = []while len(px_sample_ls) <= size:z0 = np.random.uniform(0, 4)u0 = get_px(z0)k = np.random.uniform(0, M*uniform(0, 4).pdf(z0))if k < u0:px_sample_ls.append(z0)return px_sample_lssize = 10000 px_samples = get_px_samples(size) plt.figure(figsize=(8, 6)) plt.hist(px_samples, bins=50, density=True, alpha=0.75) plt.xlim(0,4) plt.xlabel("X") plt.title("X~p(x)") plt.show()

2.4 MCMC（蒙特卡洛馬爾可夫）生成隨機數

我們在馬爾可夫鏈中又一個重要的結論：只要轉移概率矩陣合適，不管初始分布如何都會收斂到同一個平穩分布。

如果我們能找到一個轉移概率矩陣，這個矩陣得到的平穩分布剛好就是我們想要采樣的分布，那么我們可以讓平穩后的狀態就是我們所需要的采樣樣本，因為這些狀態都是由平穩分布得到的。

已知條件是平穩分布（我們需要采樣的分布），想要求一個轉移概率矩陣

隨便找一個轉移概率矩陣$Q$，$\pi (i)Q(i,j)\ne \pi (j)Q(j,i)$。如果想讓平穩條件成立，可以引入$\alpha (i,j)$，則：
$$\pi (i)Q(i,j)\alpha (i,j)=\pi (j)Q(j,i)\alpha (j,i)$$
最簡單的$\alpha (i,j)$與$\alpha (j,i)$是：
$$\begin{array}{rl}& \alpha (i,j)=\pi (j)Q(j,i)\\ & \alpha (j,i)=\pi (i)Q(i,j)\end{array}$$

$\alpha (i,j)$ 我們有一般稱之為接受率，$\alpha (i,j)$取值在 $[0,1]$ 之間, 可以理解為一個概率值（可以從均勻分布中采樣）。

最后，我們總結下MCMC采樣的過程：

• 輸入任意的轉移概率矩陣$Q(i,j)$，假設$n_1$次轉移后能夠到達平穩，需要采$n_2$個樣本。因此采樣算法從$n_1$次后才真正進行，目的是放棄平穩前的采樣值，因為只有平穩后的分布才是我們的目標采樣分布，如果采用了平穩前的樣本，這些樣本的分布將不是我們的目標分布。
• 隨機產生一個初始值$x_0$
• for $t=0$ to $n_1+n_2?1$ :
  - 從條件概率分布 $Q\left(x\mid x_t\right)$ 中采樣得到樣本 $x_{?}$
  - 從均勻分布采樣 $u～uniform[0,1]$
  - 如果 $u<\alpha \left(x_t,x_{?}\right)=\pi \left(x_{?}\right)Q\left(x_{?},x_t\right)$ ，則接受轉移 $x_t\to x_{?}$ ，即 $x_{t+1}=x_{?}$
  - 否則不接受轉移，即 $x_{t+1}=x_t$
• 樣本集 $\left(x_{n_1},x_{n_1+1},\dots ,x_{n_1+n_2?1}\right)$ （注意不能取前$n_1$個樣本）即為我們需要的平穩分布對應的隨機樣本。

exp
現在，我們使用MCMC采樣來對$p(x)=0.3e^{?(x?0.3{)}^2}+0.7e^{?\frac{(x?2.0{)}^2}{0.3}}$，其中$0\le x\le 4$進行隨機采樣：

• 假設$Q(i,j)$為$\mu $，$\sigma$的正態分布
• 初始樣本$x_0=0$

# 使用MCMC采樣來對正態分布進行隨機采樣 # 定義任意分布Q的密度函數 def Gussian_q(x,mu,sigma):return 1/(np.sqrt(2*np.pi)*sigma)*np.exp(-np.square((x-mu))/(2*np.square(sigma))) # 從任意分布Q中獲取一個樣本 def get_q_sample(mu, sigma):return np.random.normal(mu, sigma) # 定義需要采樣的分布的密度函數 def p(x):return 0.3*np.exp(-(x-0.3)**2) + 0.7* np.exp(-(x-2.0)**2/0.3) if (x>0) & (x<4) else 0def get_mcmc_normal(size):mu = 0sigma = 1x_0 = 0 # 初始的隨機狀態x0n1 = 5000n2 = sizep_sample_ls = []t = 0while len(p_sample_ls) < n2:t += 1x_1 = get_q_sample(mu, sigma)u = np.random.uniform(0, 1)if u < p(x_1)*Gussian_q(x_0, mu, sigma)/Gussian_q(x_1, mu, sigma):x_0 = x_1if t >= n1:p_sample_ls.append(x_0)return p_sample_lsp_sample_ls = get_mcmc_normal(size=10000) plt.figure(figsize=(8, 6)) plt.hist(p_sample_ls, bins=50, density=True, alpha=0.7) plt.xlim(0,4) plt.xlabel("x") plt.title("使用MCMC抽樣") plt.show()

總結

以上是生活随笔為你收集整理的GitModel|Task04|随机模拟的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。