當(dāng)前位置：首頁(yè) > 编程资源 > 编程问答 >内容正文

编程问答

GitModel|Task04|随机模拟

發(fā)布時(shí)間：2024/3/7 编程问答 37 豆豆

生活随笔收集整理的這篇文章主要介紹了 GitModel|Task04|随机模拟小編覺得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.

動(dòng)手學(xué)隨機(jī)模擬
- 1. 如何通過隨機(jī)模擬估計(jì)看漲期權(quán)的報(bào)酬分布
- - 1.1 金融知識(shí)：股票與看漲期權(quán)
  - 1.2 問題分析與確定建模思路
  - 1.3 如何模擬股票價(jià)格：布朗運(yùn)動(dòng)
- 1.4 如何更真實(shí)地模擬股票價(jià)格：幾何布朗運(yùn)動(dòng)
- - 1.5 估計(jì)看漲期權(quán)的報(bào)酬分布
  - 1.6 隨機(jī)模擬的幾個(gè)關(guān)鍵點(diǎn)
- 2.隨機(jī)模擬的關(guān)鍵：隨機(jī)數(shù)生成
- - 2.1 隨機(jī)數(shù)生成的基礎(chǔ)：均勻分布的隨機(jī)數(shù)
  - - 構(gòu)造滿足某個(gè)離散分布的隨機(jī)數(shù)
  - 2.2 分布函數(shù)生成隨機(jī)數(shù)(逆變換)
  - 2.3 拒絕采樣生成隨機(jī)數(shù)
  - 2.4 MCMC（蒙特卡洛馬爾可夫）生成隨機(jī)數(shù)

動(dòng)手學(xué)隨機(jī)模擬

1. 如何通過隨機(jī)模擬估計(jì)看漲期權(quán)的報(bào)酬分布

關(guān)鍵字:隨機(jī)模擬，看漲期權(quán)，報(bào)酬分布

1.1 金融知識(shí)：股票與看漲期權(quán)

權(quán)利金
行權(quán)價(jià)格
行權(quán)
買進(jìn)看漲期權(quán)
期權(quán)的標(biāo)的資產(chǎn)
期權(quán)買方享有權(quán)利而不用承擔(dān)義務(wù)

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 # 用來正常顯示負(fù)號(hào) price_arr = np.linspace(0, 40, 10000) # 定義價(jià)格波動(dòng)，從0～100波動(dòng) 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()

購(gòu)買了看漲期權(quán)，收益上不封頂，虧損最多就是合約金（權(quán)利金）
一般來說，權(quán)利金的數(shù)額是遠(yuǎn)遠(yuǎn)低于標(biāo)的資產(chǎn)（股票）本身的資產(chǎn)價(jià)格。換句話來說，期權(quán)是有杠桿的，現(xiàn)在你用2元的本金撬動(dòng)了幾十元的收益，這是個(gè)十分關(guān)鍵的地方

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

看漲期權(quán)能得到報(bào)酬的關(guān)鍵點(diǎn)就是標(biāo)的資產(chǎn)的價(jià)格，對(duì)于看漲期權(quán)來說，標(biāo)的資產(chǎn)的價(jià)格越高，我們獲得報(bào)酬相應(yīng)的也會(huì)越高。

問題的關(guān)鍵是:

（1）在某個(gè)時(shí)刻如何估計(jì)股票的價(jià)格分布
（2）在某個(gè)時(shí)刻結(jié)合行權(quán)價(jià)格與股票的價(jià)格分布，模擬看漲期權(quán)的報(bào)酬分布

股票價(jià)格是一個(gè)隨機(jī)過程，而看漲期權(quán)的價(jià)格（也是衍生品的價(jià)格）是隨機(jī)過程的函數(shù)
估計(jì)股票價(jià)格=>使用隨機(jī)過程近似股票價(jià)格這個(gè)隨機(jī)過程

# 獲取工商銀行股票價(jià)格數(shù)據(jù) import baostock as bs #登錄系統(tǒng) 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) #### 打印結(jié)果集 #### 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 # 查看工商銀行的價(jià)格走勢(shì)圖 plt.figure(figsize=(8,6)) plt.plot(gsyh['date'], gsyh['close'], lw=2, alpha=0.8) plt.xlabel("date") plt.ylabel("Price") plt.title("工商銀行價(jià)格走勢(shì)圖") plt.show()

我們發(fā)現(xiàn)曲線有如下特點(diǎn)：

（1）股票價(jià)格走勢(shì)具有強(qiáng)烈的隨機(jī)性；
（2）股票價(jià)格的曲線是一條"不平滑"的曲線，與我們前面學(xué)習(xí)的可微分的曲線不太一樣。

可以嘗試使用物理中的布朗運(yùn)動(dòng)去近似股票價(jià)格的無規(guī)律波動(dòng),數(shù)學(xué)上對(duì)布朗運(yùn)動(dòng)的描述發(fā)展相比于物理上的布朗運(yùn)動(dòng)要慢一些

1.3 如何模擬股票價(jià)格：布朗運(yùn)動(dòng)

設(shè)有一個(gè)粒子在直線上隨機(jī)游動(dòng), 該粒子在每個(gè)單位時(shí)間內(nèi)等可能地向左或向右移動(dòng)一個(gè)單位長(zhǎng)度，觀察這個(gè)例子的軌跡圖

# 獲取隨機(jī)游動(dòng)的位置 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("隨機(jī)游動(dòng)的位置") plt.show()

我們將粒子的隨機(jī)游動(dòng)與股價(jià)走勢(shì)放在一起對(duì)比一下：

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("隨機(jī)游動(dòng)的位置") plt.subplot(1,2,2) plt.plot(gsyh['date'], gsyh['close'], lw=2, alpha=0.8) plt.xlabel("date") plt.ylabel("Price") plt.title("工商銀行價(jià)格走勢(shì)圖") plt.show()

通過對(duì)比一維布朗運(yùn)動(dòng)走勢(shì)和股票價(jià)格的走勢(shì)可以發(fā)現(xiàn)，這兩者在形狀上非常相似

布朗運(yùn)動(dòng)詳細(xì)定義
隨機(jī)過程 ${X(t),t?0}\{X(t), t \geqslant 0\}$ 如果滿足:
(1) $X (0) = 0$ ;
(2) ${X(t),t?0}\{X(t), t \geqslant 0\}$ 有平穩(wěn)獨(dú)立增量;
(3) 對(duì)每個(gè) $t > 0, X (t)$ 服從正態(tài)分布 $N(0,σ2t)N\left(0, \sigma^{2} t\right)$ ,
則稱 ${X(t),t?0}\{X(t), t \geqslant 0\}$ 為 Brown 運(yùn)動(dòng), 也稱為 Wiener 過程, 常記為 ${B(t),t?0}\{B(t), t \geqslant 0\}$ 或 ${W(t),t?0}\{W(t), t \geqslant 0\}$ 。
如果 $σ=1\sigma=1$ , 我們稱為標(biāo)準(zhǔn) Brown 運(yùn)動(dòng); 如果 $σ≠1\sigma \neq 1$ , 則可考慮 ${X(t)/σ,t?0}\{X(t) / \sigma, t \geqslant 0\}$ , 它是標(biāo)準(zhǔn) Brown 運(yùn)動(dòng)。

我們可以驗(yàn)證一維布朗運(yùn)動(dòng)（隨機(jī)游走）的定義（3），即：均值函數(shù)和方差函數(shù)分別為：0與 $σ2t\sigma^2 t$ 以及每個(gè)時(shí)刻的位置分布為正態(tài)分布。

plt.figure(figsize=(16, 6)) plt.subplot(1,2,1) ## 模擬10000次隨機(jī)漫步 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("隨機(jī)游動(dòng)的位置") 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()

隨機(jī)游動(dòng)的均值函數(shù)為0， 方差函數(shù)是一個(gè)關(guān)于時(shí)間 $t$ 的直線函數(shù)，index為60的位置分布也符合正態(tài)分布，以上都是布朗運(yùn)動(dòng)的表現(xiàn)。

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）布朗運(yùn)動(dòng)的軌跡線會(huì)頻繁穿越 $p os = 0$ 的水平線
我們把0基準(zhǔn)線看成是股票的開盤價(jià)，頻繁穿越0基準(zhǔn)線意味著股價(jià)很大概率會(huì)在開盤價(jià)上下波動(dòng)，而不是始終維持在開盤價(jià)的上方或者下方。這一點(diǎn)的性質(zhì)在很多股票投資的實(shí)踐經(jīng)驗(yàn)都有印證
（2）在任意時(shí)刻 $t$ , 樣本軌跡的位置 $B (t)$ 不會(huì)偏離正負(fù)一個(gè)標(biāo)準(zhǔn)差太遠(yuǎn)：我們?cè)诓祭蔬\(yùn)動(dòng)的定義中的（3）可以知道，對(duì)每個(gè) $t > 0, X (t)$ 服從正態(tài)分布 $N(0,σ2t)N\left(0, \sigma^{2} t\right)$ ，一個(gè)標(biāo)準(zhǔn)差即 $σt\sigma \sqrt{t}$ ，上圖中的紅色曲線可以包含絕大多數(shù)的樣本軌跡。這點(diǎn)可以說明，隨著交易時(shí)間 $t$ 的推移， $t$ 時(shí)刻的股票價(jià)格始終不會(huì)偏離開盤價(jià)（0基準(zhǔn)線） $±t×\pm \sqrt{t} \times$ 價(jià)格波動(dòng)的標(biāo)準(zhǔn)差，就算偏離也不會(huì)太遠(yuǎn)（除了小概率事件）。
缺點(diǎn)
布朗運(yùn)動(dòng)可以是負(fù)數(shù)，因?yàn)椴祭蔬\(yùn)動(dòng)的樣本軌跡會(huì)頻繁穿越0基準(zhǔn)線。然而，股票價(jià)格卻肯定不會(huì)是負(fù)數(shù)，這是布朗運(yùn)動(dòng)的理論與現(xiàn)實(shí)的矛盾。為了解決這個(gè)問題，我們可以改進(jìn)布朗運(yùn)動(dòng)描述股票價(jià)格走勢(shì)，使用幾何布朗運(yùn)動(dòng)描述股票價(jià)格。

1.4 如何更真實(shí)地模擬股票價(jià)格：幾何布朗運(yùn)動(dòng)

為了解決使用布朗運(yùn)動(dòng)描述股票價(jià)格走勢(shì)會(huì)出現(xiàn)負(fù)數(shù)的問題，我們嘗試給布朗運(yùn)動(dòng)加上一個(gè)僅和時(shí)間 $t$ 有關(guān)的漂移項(xiàng) $μt\mu t$ , 以及一個(gè)尺度參數(shù) $σ\sigma$ ，這樣就可以得到帶漂移項(xiàng)的布朗運(yùn)動(dòng)：

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

我們假設(shè)使用 $S (t)$ 表示股票價(jià)格而 $X (t)$ 是股票價(jià)格的自然對(duì)數(shù)，即： $X (t) = l n (S (t))$ 。因此，帶漂移項(xiàng)的布朗運(yùn)動(dòng)就變成了：
$\mu t+\sigma B(t)$

我們對(duì) $\mu t+\sigma B(t)$ 進(jìn)行微分(簡(jiǎn)單的對(duì)數(shù)函數(shù)求導(dǎo))，即：
$dS(t)S(t)=μdt+σdB(t)\frac{d S(t)}{S(t)}=\mu d t+\sigma d B(t)$
$d S (t) / S (t)$ 就是這段間隔內(nèi)的收益率
化簡(jiǎn)上述的式子 $dS(t)S(t)=μdt+σdB(t)\frac{d S(t)}{S(t)}=\mu d t+\sigma d B(t)$ 得到：
$S(t)=\mu S(t) d t+\sigma S(t) d B(t)$

式子： $S(t)=\mu S(t) d t+\sigma S(t) d B(t)$ 是一個(gè)微分方程，這個(gè)微分方程是帶有隨機(jī)性的微分方程，簡(jiǎn)稱"隨機(jī)微分方程"。因此，滿足上述隨機(jī)微分方程 $S(t)=\mu S(t) d t+\sigma S(t) d B(t)$ 的股價(jià) $S (t)$ 是一個(gè)幾何布朗運(yùn)動(dòng)。

我們可以使用歐拉離散化可以得到離散時(shí)間模型：
$St=St?Δte(μ?12σ2)Δt+σεtΔtS_{t}=S_{t-\Delta t} e^{\left(\mu-\frac{1}{2} \sigma^{2}\right) \Delta t+\sigma \varepsilon_{t} \sqrt{\Delta t}}$
離散時(shí)間模型的參數(shù)含義：

$S_{t}$ : 證券在t時(shí)刻的價(jià)格;
$St?Δt\mathrm{S}_{\mathrm{t}-\Delta \mathrm{t}}$ : 證券在t $?Δt-\Delta \mathrm{t}$ 時(shí)刻的價(jià)格；
$μ\mu$ ：證券收益率的期望值;
$σ\sigma$ ：證券收益率的波動(dòng)率;
$εt\varepsilon_{t}$ : 服從正態(tài)分布（期望為 0 , 方差為 1 ）

接下來，我們模擬證券初始價(jià)格為10（日收益率均值漂移項(xiàng)為0.005，波動(dòng)率為0.25），模擬天數(shù)為300天，次數(shù)為5000次的幾何布朗運(yùn)動(dòng)價(jià)格。

# 我們模擬證券初始價(jià)格為10（日收益率均值漂移項(xiàng)為0.005，波動(dòng)率為0.25），模擬天數(shù)為300天，次數(shù)為5000次的幾何布朗運(yùn)動(dòng)價(jià)格 S0 = 10 # 初始價(jià)格 n_days = 300 # 模擬天數(shù) mu = 0.005 sigma = 0.25 delta_t = 1/n_days n_times = 5000 # 模擬次數(shù) price_arr = np.zeros((n_times, n_days)) # n_times*n_days個(gè)價(jià)格 price_arr[: , 0] = S0 # 初始價(jià)格 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條股價(jià)樣本軌跡") 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條股價(jià)樣本軌跡") plt.show()

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

以上的案例給定了離散時(shí)間模型的參數(shù)值，這是一種只會(huì)存在于習(xí)題中的情況，因?yàn)楫?dāng)我們想模擬一只具體的股票時(shí)，這些參數(shù)信息不會(huì)一下子彈到我們腦邊，我們需要根據(jù)實(shí)際數(shù)據(jù)估計(jì)其中的參數(shù)值。估計(jì)方法如下（不給證明了）：
$μ=E[ln?(St+ΔtSt)]/Δtσ2=var?[ln?(St+ΔtSt)]/Δt\begin{aligned} \mu &=E\left[\ln \left(\frac{S_{t+\Delta t}}{S_{t}}\right)\right] / \Delta t \\ \sigma^{2} &=\operatorname{var}\left[\ln \left(\frac{S_{t+\Delta t}}{S_{t}}\right)\right] / \Delta t \end{aligned}$

# 使用工商銀行股票價(jià)格數(shù)據(jù)和幾何布朗運(yùn)動(dòng)模擬股價(jià) ## 1. 獲取2021年的工商銀行股票價(jià)格數(shù)據(jù) # 登陸系統(tǒng) 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) #### 打印結(jié)果集 #### 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.估計(jì)幾何布朗運(yùn)動(dòng)中的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.使用幾何布朗運(yùn)動(dòng)模擬股價(jià)走勢(shì) print("==============開始模擬================") S0 = price_ls[-1] # 初始價(jià)格 print("初始價(jià)格：", S0) n_days = 300 # 模擬天數(shù) sigma = np.sqrt(sigma2) print("估計(jì)的mu = ", mu) print("估計(jì)的sigma = ", sigma) delta_t = 1/n_days n_times = 10000 # 模擬次數(shù) price_arr = np.zeros((n_times, n_days)) # n_times*n_days個(gè)價(jià)格 price_arr[: , 0] = S0 # 初始價(jià)格 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條股價(jià)樣本軌跡") 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條股價(jià)樣本軌跡") plt.show()

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

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()

小節(jié)：我們通過帶漂移項(xiàng)的布朗運(yùn)動(dòng) $X(t)=μt+σB(t)X(t)=\mu t+\sigma B(t)$ 可以知道， $X (t)$ 的期望為 $E[X(t)]=μtE[X(t)]=\mu t$ ，即t時(shí)刻的價(jià)格期望為 $μt\mu t$ 。這是個(gè)十分有趣的結(jié)論，因?yàn)榧偃?span id="ze8trgl8bvbq" class="katex--inline"> $μ>0\mu > 0$ ，那么t時(shí)刻的價(jià)格期望也會(huì)大于0。因此如果我們堅(jiān)信股市長(zhǎng)期來看是一個(gè)向上的牛市（雖然很緩慢），那么我們就應(yīng)該接受短期波動(dòng)而堅(jiān)持長(zhǎng)期持股，就像巴菲特的價(jià)值投資一樣。

1.5 估計(jì)看漲期權(quán)的報(bào)酬分布

現(xiàn)在我們假設(shè)行權(quán)價(jià)格是4.7元，我們來分析報(bào)酬在第300天的時(shí)候的報(bào)酬分布情況（忽略權(quán)利金、遠(yuǎn)期利率和股息）：
當(dāng)股票價(jià)格大于行權(quán)價(jià)格，那么我將獲利；
當(dāng)股票價(jià)格小于行權(quán)價(jià)格，那么我將虧損；

# 繪制報(bào)酬分布圖 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] #剔除負(fù)數(shù) plt.figure(figsize=(8,6)) plt.hist(profit_sum, bins=40, density=True) plt.xlabel("Profit") plt.show()

1.6 隨機(jī)模擬的幾個(gè)關(guān)鍵點(diǎn)

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

（1）確定事件的過程（仿真）；
（2）在事件中添加合適的分布的隨機(jī)數(shù)（隨機(jī)模擬）；（關(guān)鍵）
（3）得到模擬的結(jié)論；
采樣算法：（拒絕采樣、MCMC采樣）
（因?yàn)橐獜哪硞€(gè)分布中采樣，采樣的樣本就是所謂的隨機(jī)數(shù)）

2.隨機(jī)模擬的關(guān)鍵：隨機(jī)數(shù)生成

計(jì)算機(jī)生成隨機(jī)數(shù)的基礎(chǔ)算法：均勻分布采樣
分布函數(shù)采樣
拒絕采樣
MCMC采樣

2.1 隨機(jī)數(shù)生成的基礎(chǔ)：均勻分布的隨機(jī)數(shù)

偽隨機(jī)數(shù)生成的方法有很多, 比較簡(jiǎn)單的一種方式比如:
$xn+1=(axn+c)modm\mathrm{x}_{\mathrm{n}+1}=(\mathrm{ax_n} +\mathrm{c}) \bmod \mathrm{m}$

其中, $X$ 為偽隨機(jī)序列,

$m, m > 0$ , 模數(shù), 也是生成序列的最大周期
$a, 0 < a < m$ , 乘數(shù)
$\leq c<m$ , 增量
$X0,0≤X0<mX_{0}, 0 \leq X_{0}<m$ , 種子點(diǎn) (起始值)

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

這種方法叫線性同余法，使用線性同余法時(shí)需要特別小心參數(shù)設(shè)置，否則生成的隨機(jī)數(shù)將會(huì)十分糟糕。下面通過一組實(shí)驗(yàn)看看不同的參數(shù)設(shè)置下生成的隨機(jī)數(shù)序列：

# 線性同余法 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 # 實(shí)驗(yàn)次數(shù) ## 第一組實(shí)驗(yàn)： a, c, m, x0 = 11, 0, 8, 1 uniform_ls = get_uniform_sample(a, c, m, x0, n_times) print("第一次實(shí)驗(yàn)：") print("a = %d, c = %d, m = %d, x0 = %d:\n" % (a, c, m, x0)) print(uniform_ls) print("============================") ## 第二組實(shí)驗(yàn)： a, c, m, x0 = 25214903917, 11, 2**48, 1 # 這也是java.util.Random中的默認(rèn)參數(shù)設(shè)置 uniform_ls = get_uniform_sample(a, c, m, x0, n_times) print("第二次實(shí)驗(yàn)：") print("a = %d, c = %d, m = %d, x0 = %d:\n" % (a, c, m, x0)) print(uniform_ls)

構(gòu)造滿足某個(gè)離散分布的隨機(jī)數(shù)

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

生成一個(gè) $[0, 1]$ 的均勻分布的隨機(jī)數(shù) $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$ 。
重復(fù)以上步驟n_times次就能獲得n_times個(gè)符合離散分布的隨機(jī)數(shù)

# 生成10000個(gè)離散分布的隨機(jī)數(shù) from collections import Counter def get_discrete_sample(x_ls, p_ls, size):""":param x_ls: 類別序列:param p_ls: 類別的概率序列:param size: 生成隨機(jī)數(shù)個(gè)數(shù):return: 隨機(jī)數(shù)序列"""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 分布函數(shù)生成隨機(jī)數(shù)(逆變換)

通過使用均勻分布的隨機(jī)數(shù)構(gòu)造別的分布的隨機(jī)數(shù)呢？答案是肯定的，這個(gè)方法叫逆變換法。

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

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

指數(shù)分布的概率密度函數(shù)為 $f_{X}(x)=$ $λe?λx,x>0\lambda e^{-\lambda x}, x>0$ ，它的累積分布函數(shù)為
$FX(x)=1?e?λx,x≥0\mathrm{F}_{X}(x)=1-e^{-\lambda x}, \quad \mathrm{x} \geq 0$
產(chǎn)生 $(0, 1)$ 均勻分布的隨機(jī)數(shù) $V$ （當(dāng)作是指數(shù)分布的分布函數(shù)值）
令 $Y=?1λlog?(V)Y=-\frac{1}{\lambda} \log (V)$ （將滿足均勻分布的分布函數(shù)進(jìn)行逆變換）， $Y$ 即為服從指數(shù)分布的隨機(jī)數(shù)

逆變換法構(gòu)造指數(shù)分布的隨機(jī)數(shù) 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 拒絕采樣生成隨機(jī)數(shù)

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

$p (x)$ 是我們需要采樣的分布，有時(shí)候十分復(fù)雜，這里的 $e^{-(x - 0.3)^2} + 0.7e^{-\frac{(x - 2.0)^2}{0.3}}$ ，其中 $\le x \le 4$
$q (x)$ 是我們的構(gòu)造的分布，這里我們假設(shè)為[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)$ 太復(fù)雜，我們無法直接采樣，因此我們借助輔助的分布 $q (x)$ ，這個(gè)輔助的分布滿足一個(gè)特點(diǎn)： $q (x)$ 乘一個(gè)常數(shù) $M$ 后的 $Mq (x)$ 曲線能完全包住我們想要采樣的分布 $p (x)$ 。

拒絕采樣的過程如下：

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

# 通過拒絕采樣0.3*np.exp(-(x-0.3)**2)+0.7*np.exp(-(x-2.0)**2/0.3)分布的隨機(jī)數(shù) 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（蒙特卡洛馬爾可夫）生成隨機(jī)數(shù)

我們?cè)隈R爾可夫鏈中又一個(gè)重要的結(jié)論：只要轉(zhuǎn)移概率矩陣合適，不管初始分布如何都會(huì)收斂到同一個(gè)平穩(wěn)分布。

如果我們能找到一個(gè)轉(zhuǎn)移概率矩陣，這個(gè)矩陣得到的平穩(wěn)分布剛好就是我們想要采樣的分布，那么我們可以讓平穩(wěn)后的狀態(tài)就是我們所需要的采樣樣本，因?yàn)檫@些狀態(tài)都是由平穩(wěn)分布得到的。

已知條件是平穩(wěn)分布（我們需要采樣的分布），想要求一個(gè)轉(zhuǎn)移概率矩陣

隨便找一個(gè)轉(zhuǎn)移概率矩陣 $Q$ ， $π(i)Q(i,j)≠π(j)Q(j,i)\pi(i) Q(i, j) \neq \pi(j) Q(j, i)$ 。如果想讓平穩(wěn)條件成立，可以引入 $α(i,j)\alpha(i, j)$ ，則：
$π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)\pi(i) Q(i, j) \alpha(i, j)=\pi(j) Q(j, i) \alpha(j, i)$
最簡(jiǎn)單的 $α(i,j)\alpha(i, j)$ 與 $α(j,i)\alpha(j, i)$ 是：
$α(i,j)=π(j)Q(j,i)α(j,i)=π(i)Q(i,j)\begin{aligned} &\alpha(i, j)=\pi(j) Q(j, i) \\ &\alpha(j, i)=\pi(i) Q(i, j) \end{aligned}$

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

最后，我們總結(jié)下MCMC采樣的過程：

輸入任意的轉(zhuǎn)移概率矩陣 $Q (i, j)$ ，假設(shè) $n_1$ 次轉(zhuǎn)移后能夠到達(dá)平穩(wěn)，需要采 $n_2$ 個(gè)樣本。因此采樣算法從 $n_1$ 次后才真正進(jìn)行，目的是放棄平穩(wěn)前的采樣值，因?yàn)橹挥衅椒€(wěn)后的分布才是我們的目標(biāo)采樣分布，如果采用了平穩(wěn)前的樣本，這些樣本的分布將不是我們的目標(biāo)分布。
隨機(jī)產(chǎn)生一個(gè)初始值 $x_0$
for $t = 0$ to $n_{1}+n_{2}-1$ :
- 從條件概率分布 $Q(x∣xt)Q\left(x \mid x_{t}\right)$ 中采樣得到樣本 $x_{*}$
- 從均勻分布采樣 $\sim uniform[0,1]$
- 如果 $u<α(xt,x?)=π(x?)Q(x?,xt)u<\alpha\left(x_{t}, x_{*}\right)=\pi\left(x_{*}\right) Q\left(x_{*}, x_{t}\right)$ ，則接受轉(zhuǎn)移 $xt→x?x_{t} \rightarrow x_{*}$ ，即 $x_{t+1}=x_{*}$
- 否則不接受轉(zhuǎn)移，即 $x_{t+1}=x_{t}$
樣本集 $(xn1,xn1+1,…,xn1+n2?1)\left(x_{n_{1}}, x_{n_{1}+1}, \ldots, x_{n_{1}+n_{2}-1}\right)$ （注意不能取前 $n_1$ 個(gè)樣本）即為我們需要的平穩(wěn)分布對(duì)應(yīng)的隨機(jī)樣本。

exp
現(xiàn)在，我們使用MCMC采樣來對(duì) $e^{-(x - 0.3)^2} + 0.7e^{-\frac{(x - 2.0)^2}{0.3}}$ ，其中 $\le x \le 4$ 進(jìn)行隨機(jī)采樣：

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

# 使用MCMC采樣來對(duì)正態(tài)分布進(jìn)行隨機(jī)采樣 # 定義任意分布Q的密度函數(shù) 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中獲取一個(gè)樣本 def get_q_sample(mu, sigma):return np.random.normal(mu, sigma) # 定義需要采樣的分布的密度函數(shù) 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 # 初始的隨機(jī)狀態(tài)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()

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。

gitmodel

上一篇：使用Java实现Comet风格的Web应
下一篇：基带信号、载波、带通信号