内容總結
- 學習datawhale的gitmodel教程。小郭為了鎖定價格波動風險，簽訂合約即買進看漲期權：提前給榴蓮超市2塊權利金，現在榴蓮30元一塊（期權的标的資産），下個月能用20元買到一塊榴蓮（即到期行使合約權利）。
- 期權是有杠杆的。
- 看漲期權能得到報酬的關鍵點就是标的資産的價格，對于看漲期權來說，标的資産的價格越高，我們獲得報酬相應的也會越高。
- 幾何布朗運動：
- 為了解決使用布朗運動描述股票價格走勢會出現負數的問題，我們嘗試給布朗運動加上一個僅和時間有關的漂移項, 以及一個尺度參數，這樣就可以得到帶漂移項的布朗運動：。但是因為與的取值随着時間
- 雖然股票價格不能是負數，但是股票價格的自然對數可以是負數。是以，我們假設使用表示股票價格（注意現在不是了），而是股票價格的自然對數，即：。然後對進行微分。
- 假設 為股票的價格, 則 為股價在無窮小的時間間隔内的變化量, 而 。化簡上述的式子得到帶随機性的微分方程，滿足該方程的股價即滿足一個幾何布朗運動（最後還要歐拉離散化得到離散時間模型）：
- 通過帶漂移項的布朗運動可以知道，的期望為，即t時刻的價格期望為。假如，那麼t時刻的價格期望也會大于0。是以如果我們堅信股市長期來看是一個向上的牛市（雖然很緩慢），那麼我們就應該接受短期波動而堅持長期持股，就像巴菲特的價值投資一樣。
- 很多随機模拟場景中需要構造自己想要的分布的随機數（采樣算法，從某個分布中采樣，采樣的樣本就是所謂的随機數），如：拒絕采樣、MCMC采樣等等。
文章目錄
- 内容總結
- 一、如何通過随機模拟估計看漲期權的報酬分布
- 1.1 金融知識：股票與看漲期權
- 1.2 問題分析與确定模組化思路
- 1.3 如何模拟股票價格：布朗運動
- 1.4 如何更真實地模拟股票價格：幾何布朗運動
- 1.5 估計看漲期權的報酬分布
- 1.6 随機模拟的幾個關鍵點
- 二、随機模拟的關鍵：随機數生成
- 2.1 随機數生成的基礎：均勻分布的随機數
- 2.2 分布函數生成随機數
- 2.3 拒絕采樣生成随機數
- 2.4 MCMC（蒙特卡洛馬爾可夫）生成随機數
- 附：時間安排
- Reference
一、如何通過随機模拟估計看漲期權的報酬分布
标題的三個關鍵詞：随機模拟、看漲期權、報酬分布。
這一個研究随機現象的問題，因為它想要我們研究“報酬分布”，先了解看漲期權，再定義看漲期權的報酬。
1.1 金融知識：股票與看漲期權
【榴蓮栗子】
假如你是一位榴蓮愛好者，衆所周知，榴蓮的價格波動是十分明顯的。現在榴蓮的價格是30元一塊，為了能在8月吃到美味的榴蓮，你跟GitModel超市簽訂了協定，協定規定：現在你必須給GitModel超市2元，那麼在8月的任何一天，你都能用20元買到一塊榴蓮。你認為，未來的榴蓮價格肯定會高于元，是以，你高高興興地簽了協定。是以，在8月份會出現一下幾種情況：
- 榴蓮價格在22元以上，你将以元的價格買到高于22元價值的榴蓮，你将會賺元。
- 榴蓮價格在20元以下，你可以選擇不行使購買的權利，你會從别的超市以的價格買來榴蓮，您将虧損2元。
- 榴蓮價格在20元至22元之間，你為了能吃到榴蓮，你也可以行權以20元的價格買到榴蓮。
隻要榴蓮的市價比總花費元高，那你就賺到了。而如果榴蓮的市價比總花費元低，或者等于總的花費元，那你可以選擇不行權，最多就虧損2元的合約費。
【看漲期權的作用】
由于你現在無法判斷榴蓮未來價格的漲跌，而期權的出現就可以幫你鎖定價格波動風險，如果榴蓮的價格一直上漲，你獲利的空間就會越來越大；但是如果榴蓮的價格一直下跌，因為期權買方（購買看漲期權的人）享有權利而不用承擔義務（相反，期權賣方（跟你簽訂看漲期權的人）承擔義務但不享有權利），是以你最大的損失就是2元的合約費。
【幾個專業術語】
- 簽訂合約，即買進看漲期權
- 期權的标的資産：榴蓮
- 權利金：2元的合約費
- 行權：到期行使合約權利
- 行權價格：20元
我們可以來對上面的看漲期權的例子的收益曲線進行繪制：
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_cost
my_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元的本金撬動了幾十元的收益：
- 假如現在房子是100萬，你用100萬買下房子，一段時間過後房子漲到了150萬，你賺了萬，收益率為：
- 假如你購買了一張看漲期權，價值為2萬，規定未來無論該房子價格如何，都可以使用100萬購買該房子。一段時間後，房子漲到了150萬，你想行使權利，但是你沒有100萬的交易費用怎麼辦？這時候你可以問銀行借100萬（假如不算利息），用100萬買了150萬的房子，再把房子以150萬價格賣出去，這時候你把得來的150萬中的100萬還給銀行，剩下50萬的盈利。這時候，我們計算下收益率：，也就是24倍，不是24%。是以，你是用2萬的本金撬動了48萬的盈利，但是你就算怎麼虧損，最多也就是虧掉2萬的權利金，但是虧損率是100%，這就是加了杠杆。
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() 
以上價格走勢圖的特點：
- 股票價格走勢具有強烈的随機性；
- 股票價格的曲線是一條"不平滑"的曲線，不可微分。
為了近似股票價格的性質，我們想到了實體世界中的布朗運動，什麼是實體中的布朗運動呢？1827年英國植物學家羅伯特 • 布朗（Robert Brown）在使用顯微鏡觀察水中花粉微粒運動時發現了微粒的無規則運動，這種無規則運動稱為：布朗運動。
布朗運動的無規律運動的性質很類似股票價格的無規則波動，是以可以嘗試使用實體中的布朗運動去近似股票價格的無規律波動。而數學上對布朗運動的描述發展相比于實體上的布朗運動要慢一些。布朗運動的定義和描述是由諾伯特 • 維納（Norbert Wiener）在 1918 年提出，是以布朗運動（随機過程）（Brownian motion）又稱為維納過程（Wiener process）。
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_ls
n_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() 通過對比一維布朗運動（随機遊動）走勢和股票價格的走勢可以發現，這兩者在形狀上非常相似，或許我們可以将股票價格走勢與一維布朗運動走勢聯系起來。
最早使用布朗運動分析股票價格走勢的人是路易斯 • 巴舍利耶，他在1900年就在博士論文中寫到如何使用布朗運動描述股票價格走勢和期權的價格。
【Bronw運動的定義】随機過程 如果滿足:
(1) ;
(2) 有平穩獨立增量;
(3) 對每個 服從正态分布 ,
則稱 為 Brown 運動, 也稱為 Wiener 過程, 常記為 或
如果 , 我們稱為标準 Brown 運動;
如果 , 則可考慮 , 它是标準 Brown 運動。
我們可以驗證一維布朗運動（随機遊走）的定義（3），即：均值函數和方差函數分别為：0與以及每個時刻的位置分布為正态分布。
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， 方差函數是一個關于時間的直線函數，index為60的位置分布也符合正态分布，以上都是布朗運動的表現。
通過實驗（非數學推導）的形式更加直覺地來說明為什麼使用布朗運動描述股票價格的走勢是合理的。首先，我們先來繪制20條随機遊動的樣本軌迹：
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() 從上圖可以看到：
- 布朗運動的軌迹線會頻繁穿越的水準線：圖中的12條樣本軌迹線大多數會頻繁穿越0基準線，隻有少部分曲線會不穿越0基準線（如：亮藍色的軌迹線），但是随着時間的推移，這條亮藍色的軌迹線也會穿越0基準線。我們把0基準線看成是股票的開盤價，頻繁穿越0基準線意味着股價很大機率會在開盤價上下波動，而不是始終維持在開盤價的上方或者下方。這一點的性質在很多股票投資的實踐經驗都有印證，例如：現在有一種引誘投資者投資的套路：誘多（如果大盤平穩的話，股價從開盤價一路飙升而收盤價又回到開盤價的套路，這種套路可以很輕易地把當日追漲者的資金套住）。
- 在任意時刻, 樣本軌迹的位置不會偏離正負一個标準差太遠：我們在布朗運動的定義中的（3）可以知道，對每個服從正态分布，一個标準差即，上圖中的紅色曲線可以包含絕大多數的樣本軌迹。這點可以說明，随着交易時間的推移，時刻的股票價格始終不會偏離開盤價（0基準線）
【基礎金融知識補充】
- 大盤：是指滬市的“上證綜合指數”和深市的“深證成分股指數”的股票。大盤指數是運用統計學中的指數方法編制而成的，反映股市總體價格或某類股價變動和走勢的名額。
- 上證綜合指數：以上海證券交易所 挂牌上市的全部股票（包括A股和B股）為樣本，以發行量 為權數（包括流通股本和非流通股本），以權重平均法 計算，以1990年12月19日為基日，基日指數 定為100點的股價指數。
- 深圳成分股指數：從深圳證券 交易所挂牌上市的 所有股票中抽取具有市場代表性的40家上市公司的股票為樣本，以流通股本為權數，以權重平 均法計算，以1994 年7月20日為基日，基日指數定為1000點的股價指數。 [1
上面的内容，我們主要闡述了布朗運動描述股價走勢的合理性，但是其不合理性也很明顯。其中，最大的一個缺點就是布朗運動可以是負數，因為布朗運動的樣本軌迹會頻繁穿越0基準線。然而，股票價格卻肯定不會是負數，這是布朗運動的理論與現實的沖突。為了解決這個問題，我們可以改進布朗運動描述股票價格走勢，使用幾何布朗運動描述股票價格。
1.4 如何更真實地模拟股票價格：幾何布朗運動
（1）加上一個僅和時間 有關的漂移項 ：為了解決使用布朗運動描述股票價格走勢會出現負數的問題，給布朗運動加上一個僅和時間 有關的漂移項 , 以及一個尺度參數 ，這樣就可以得到帶漂移項的布朗運動：。
即使是帶漂移項和尺度參數的布朗運動 , 它也不是描述股票價格走勢的最優選擇，因為 與 的取值随着時間
（2）使用對數性質：雖然股票價格不能是負數，但是股票價格的自然對數可以是負數。是以，我們假設使用表示股票價格（注意現在不是了），而是股票價格的自然對數，即：。是以，帶漂移項的布朗運動就變成了：
（3）微分求導：對進行微分(簡單的對數函數求導)，即：
為什麼要對進行微分呢？觀察上面的微分式子發現：假設 為股票的價格, 則 為股價在無窮小的時間間隔内的變化量, 而 ！化簡上述的式子得到：
式子：是一個微分方程，這個微分方程是帶有随機性的微分方程，簡稱"随機微分方程"。是以，滿足上述随機微分方程的股價 。
（4）歐拉離散化：現在我們就能用幾何布朗運動描述股票價格走勢了，但是這個随機微分方程怎麼描述股價走勢呢，既然連續的麻煩，我們可以使用歐拉離散化可以得到離散時間模型：
離散時間模型的參數含義：
- : 證券在t
接下來，我們模拟證券初始價格為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() 以上的案例給定了離散時間模型的參數值，這是一種隻會存在于習題中的情況，因為當我們想模拟一隻具體的股票時，這些參數資訊不會一下子彈到我們腦邊，我們需要根據實際資料估計其中的參數值。估計方法如下（不給證明了）：
# 使用工商銀行股票價格資料和幾何布朗運動模拟股價
## 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()
'''
login success!
login respond error_code:0
login respond error_msg:success
query_history_k_data_plus respond error_code:0
query_history_k_data_plus respond error_msg:success
==============開始模拟================
初始價格： 4.63
估計的mu = -0.07920499347383164
估計的sigma = 0.13628354313803137
''' # 第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() 上面完成了使用幾何布朗運動模拟股價走勢的任務。
最後，我們使用一個小案例結束本小節：我們通過帶漂移項的布朗運動可以知道， 的期望為 ，即t時刻的價格期望為。這是個十分有趣的結論，因為假如，那麼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 随機模拟的幾個關鍵點
我們回顧以下"看漲期權的報酬分布"這個案例的步驟，無論是布朗運動還是幾何布朗運動，都需要大量随機數的參與，如：布朗運動中随機遊動的案例的控制向左還是向右走的0-1分布的随機數、幾何布朗運動中控制随機項的服從正态分布的随機數。是以，我們可以将随機模拟分解為以下步驟：
- （1）确定事件的過程（仿真）；
- （2）在事件中添加合适的分布的随機數（随機模拟）；（關鍵）
- （3）得到模拟的結論；
那麼，在實際的程式設計中很多随機數都能通過python第三方包如：
numpy
、
scipy
等給出，這些第三方包給出的随機數分布往往是常用的，如：均勻分布、正态分布、指數分布等等。
但是在很多随機模拟的案例中，往往需要構造自己想要的分布的随機數，這類構造随機數的方法也叫采樣算法（因為要從某個分布中采樣，采樣的樣本就是所謂的随機數），如：拒絕采樣、MCMC采樣等等。
二、随機模拟的關鍵：随機數生成
解決上個任務中總結出來的規律：随機模拟的關鍵點在于給出滿足某個分布的随機數。構造随機數的方法非常多，這類方法也叫做采樣算法。
- 計算機生成随機數的基礎算法：均勻分布采樣
- 分布函數采樣
- 拒絕采樣
- MCMC采樣
2.1 随機數生成的基礎：均勻分布的随機數
一般來說，均勻分布是最簡單但同時是最重要的分布形式，如何對均勻分布采樣直接影響着其他分布的采樣。在計算機中生成[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) % m
uniform_ls.append(uniform_sample)
x0 = uniform_sample
return uniform_ls
n_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)
# 在numpy中，我們隻需要np.random.uniform(low, high, size)就可以構造size個low～high的均勻分布的随機數
uniform_ls = np.random.uniform(0, 1, 10)
print(uniform_ls)
'''
第一次實驗：
a = 11, c = 0, m = 8, x0 = 1:
[3, 1, 3, 1, 3, 1, 3, 1, 3, 1]
============================
第二次實驗：
a = 25214903917, c = 11, m = 281474976710656, x0 = 1:
[25214903928, 206026503483683, 245470556921330, 105707381795861, 223576932655868, 102497929776471, 87262199322646, 266094224901481, 44061996164032, 147838658590923]
[0.71364912 0.59925289 0.10660613 0.38554986 0.05717685 0.48277086
0.85735629 0.59281279 0.62027365 0.40008649]
''' 上面的線性同餘法是構造均勻分布的随機數,如果我們需要構造滿足某個離散分布的随機數，如：分别滿足機率這個離散分布的随機數。其實這個問題的解決方法并不難，需要用到均勻分布的随機數：
- 生成一個的均勻分布的随機數；
- 如果在，那麼， 如果在，那麼，如果在，那麼，如果在，那麼。
- 重複以上步驟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_ls
size = 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, 0.2982), (4, 0.1008), (3, 0.4992), (1, 0.1018)] 2.2 分布函數生成随機數
剛剛的線性同餘法隻是構造均勻分布的随機數，那能否像構造離散分布的随機數一樣，通過使用均勻分布的随機數構造别的分布的随機數呢？答案是肯定的，這個方法叫逆變換法。
逆變換法的理論依據是：假設 是一個連續型随機變量，它的累積分布函數為 ； 現在，定義随機變量 。那麼， 服從 上的均勻分布，即
換句話說，任意分布的分布函數值服從上的均勻分布，那麼我們可以通過這個理論通過均勻分布構造其他分布的随機數。
接下來，我們嘗試使用逆變換法通過均勻分布的随機數構造指數分布的随機數：
- 指數分布的機率密度函數為 ，它的累積分布函數為
- 産生 均勻分布的随機數 （當作是指數分布的分布函數值）
- 令（将滿足均勻分布的分布函數進行逆變換），
# 逆變換法構造指數分布的随機數
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_samples
lmd = 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 拒絕采樣生成随機數
逆變換法生成均勻分布外的采樣十分好用，但是這裡需要知道采樣分布的分布函數的形式，因為如果不知道采樣分布的分布函數，那麼将無法使用逆變換獲得随機數。實際上，有沒有分布的分布函數的表達不出來的呢？正态分布的分布函數就無法顯式表達，隻能用積分表達式。同時，逆變換法需要求分布函數和逆變換，方法太複雜。拒絕采樣可以輕松解決這些問題，拒絕采樣開始采用蒙特卡洛的思想構造随機數，我們先來看看下圖：
- 是我們需要采樣的分布，有時候十分複雜，這裡的，其中
- 是我們的構造的分布，這裡我們假設為[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() 我們從上圖可以看到，我們想要抽樣的分布是藍色的曲線，但是由于太複雜，我們無法直接采樣，是以我們借助輔助的分布，這個輔助的分布滿足一個特點：乘一個常數後的曲線能完全包住我們想要采樣的分布。拒絕采樣的過程如下：
- 通過輔助分布采樣，輔助分布可以選擇常見的均勻分布或者正态分布；
- 計算（圖中紅色的點）；
- 在區間範圍的均勻分布随機采樣，即在紅色的直線範圍按均勻分布采樣；
- 判斷的大小，如果（即位于下方）,則接受是來自分布的随機數；反之，如果如果（即位于上方）,則拒絕是來自分布的随機數。
- 重複以上過程，知道采樣個數為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 = 4
px_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_ls
size = 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（蒙特卡洛馬爾可夫）生成随機數
在前面随機過程的學習中，我們在馬爾可夫鍊中又一個重要的結論：隻要轉移機率矩陣合适，不管初始分布如何都會收斂到同一個平穩分布。那麼，我們會想到：能不能讓平穩後的平穩分布剛好就是我們需要采樣的分布呢？這樣做有什麼好處呢？事實上，如果我們能找到一個轉移機率矩陣，這個矩陣得到的平穩分布剛好就是我們想要采樣的分布，那麼我們可以讓平穩後的狀态就是我們所需要的采樣樣本，因為這些狀态都是由平穩分布得到的。
總結一下，我們目前具備的已知條件是平穩分布（我們需要采樣的分布），想要求一個轉移機率矩陣。直接求肯定是行不通的，怎麼可能由資訊量少的平穩分布推出資訊量多的轉移機率矩陣呢？那怎麼辦呢？下面，我們來解決這個問題，假設平穩分布為，轉移機率矩陣為。是以，如果是随便找一個轉移機率矩陣，。如果想讓平穩條件成立，可以引入，則：
最簡單的與是：
這樣, 我們就得到了我們想要采樣的分布 對應的馬爾科夫鍊狀态轉移矩陣 , 滿足:
是以， 我們有一般稱之為接受率，取值在
最後，我們總結下MCMC采樣的過程：
- 輸入任意的轉移機率矩陣，假設次轉移後能夠到達平穩，需要采個樣本。是以采樣算法從次後才真正進行，目的是放棄平穩前的采樣值，因為隻有平穩後的分布才是我們的目标采樣分布，如果采用了平穩前的樣本，這些樣本的分布将不是我們的目标分布。
- 随機産生一個初始值
-
forto:
- 從條件機率分布中采樣得到樣本
- 從均勻分布采樣
- 如果，則接受轉移，即
- 否則不接受轉移，即
- 樣本集（注意不能取前個樣本）即為我們需要的平穩分布對應的随機樣本。
現在，我們使用MCMC采樣來對，其中進行随機采樣：
- 假設為$\mu\sigma$的正态分布
- 初始樣本
# 使用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 0
def get_mcmc_normal(size):
mu = 0
sigma = 1
x_0 = 0 # 初始的随機狀态x0
n1 = 5000
n2 = size
p_sample_ls = []
t = 0
while len(p_sample_ls) < n2:
t += 1
x_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_1
if t >= n1:
p_sample_ls.append(x_0)
return p_sample_ls
p_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() 現在，我們已經幾乎能建立任何的分布的随機數，但是基于MCMC算法延伸的其他采樣算法還在繼續流行中，如：M-H采樣、吉布斯采樣等等。
附：時間安排
| Task | 内容 | 時間 |
| Task01 | 假設檢驗1：方法論與一進制數值檢驗 | 8-13 —— 8-18周四 |
| Task02 | 假設檢驗2：多元數值向量檢驗 | 8-19 —— 8-20周六 |
| Task03 | 假設檢驗3：分類資料檢驗 | 8-21 —— 8-22周一 |
| Task04 | 應用随機過程與仿真系統 | 8-23 —— 8-25周四 |
| Task05 | 金融量化分析與随機模拟 |