20 随机抽样
20.1 Inverse CDF Transformation（逆CDF变化法）
假设x的概率密度函数为f(x)，累计分布函数CDF为F(x)。我们可以通过如下步骤得到非均匀分布的随机样本。
图4.1: …
20.2 accept-reject sampling method（接受-拒绝抽样）
基本思想：假设p(x)不可以直接抽样，那么找一个可以直接抽样的分布，称为建议分布（proposal distribution） 。假设q(x)是建议分布的概率密度函数，并且有 kq(x)，k>0总在p(x)的上方。如图：
图20.1: …
具体算法如下：
图20.2: …
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams['font.sans-serif'] = ['FangSong']
plt.rcParams['axes.unicode_minus'] = False
def p(x): # 假设f1为我们想要进行抽样的分布
return 0.3/np.sqrt(np.pi*2)*np.exp(-(x+2)**2/2) + 0.7/np.sqrt(np.pi*2)*np.exp(-(x-4)**2/2)
def q(x,a,b): #定义建议分布q(x)
return 1/(b-a)
#初始设置
n=int(1e+6)
k=6
u = np.random.uniform(0,1,n)
z = np.random.uniform(-10,10,n)
#数据准备
data = {
'u':u,
'z':z,
'pz':p(z),
"kqz":k*q(z,-10,10)
}
df=pd.DataFrame(data)
df = df[df.u < df.pz/df.kqz]
z = df.loc[:,'z']
x = np.arange(-10,10,0.1)
y = p(x)
plt.plot([-10, 10], [k*1/10,k*1/10], label='建议分布') # 画出建议分布的图像
plt.plot(x,y,"b-")
plt.hist(z, bins=200, density=True,label='抽样数据', color='green')## (array([8.17626069e-05, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
## 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 8.17626069e-05,
## 1.63525214e-04, 8.17626069e-05, 1.63525214e-04, 0.00000000e+00,
## 0.00000000e+00, 3.27050428e-04, 1.63525214e-04, 4.90575642e-04,
## 1.06291389e-03, 9.81151283e-04, 1.14467650e-03, 1.38996432e-03,
## 1.55348953e-03, 1.47172693e-03, 2.69816603e-03, 3.51579210e-03,
## 4.08813035e-03, 3.18874167e-03, 4.82399381e-03, 6.54100856e-03,
## 6.29572073e-03, 9.56622501e-03, 1.12014772e-02, 1.08744267e-02,
## 1.38178806e-02, 1.73336727e-02, 1.85601118e-02, 2.29752926e-02,
## 2.44470195e-02, 2.46923073e-02, 3.24597550e-02, 3.58937845e-02,
## 3.68749357e-02, 4.56235347e-02, 5.48627093e-02, 5.44538962e-02,
## 6.14854804e-02, 6.40201212e-02, 6.64729994e-02, 7.26051950e-02,
## 8.15990817e-02, 9.07564937e-02, 9.08382563e-02, 9.13288320e-02,
## 1.00731532e-01, 1.04329086e-01, 1.06373152e-01, 1.16593478e-01,
## 1.13895311e-01, 1.21662759e-01, 1.16593478e-01, 1.14549412e-01,
## 1.18392255e-01, 1.19455169e-01, 1.25423839e-01, 1.15694089e-01,
## 1.07844879e-01, 1.15530564e-01, 1.10461282e-01, 1.03347935e-01,
## 1.01303870e-01, 9.90145170e-02, 8.91212416e-02, 8.55236869e-02,
## 8.49513486e-02, 7.29322454e-02, 7.10517054e-02, 6.28754447e-02,
## 6.09949048e-02, 5.51079971e-02, 5.14286798e-02, 4.76675999e-02,
## 4.07995409e-02, 3.84284253e-02, 3.47491080e-02, 2.91892507e-02,
## 2.81263368e-02, 2.19123787e-02, 1.92142126e-02, 1.76607231e-02,
## 1.43902188e-02, 1.35725928e-02, 9.81151283e-03, 9.32093719e-03,
## 6.94982159e-03, 6.13219552e-03, 6.37748334e-03, 5.06928163e-03,
## 4.98751902e-03, 3.02521646e-03, 5.39633206e-03, 3.51579210e-03,
## 3.76107992e-03, 4.25165556e-03, 5.64161988e-03, 6.21395813e-03,
## 7.84921027e-03, 8.42154852e-03, 9.97503805e-03, 1.32455423e-02,
## 1.30002545e-02, 1.58619457e-02, 1.83148240e-02, 2.07677022e-02,
## 2.35476308e-02, 3.13150785e-02, 3.50761584e-02, 4.31706565e-02,
## 4.63593981e-02, 5.41268458e-02, 6.26301569e-02, 6.74541507e-02,
## 8.15173191e-02, 9.15741198e-02, 1.03184410e-01, 1.10788332e-01,
## 1.15939377e-01, 1.32864236e-01, 1.42348699e-01, 1.54694852e-01,
## 1.68839783e-01, 1.79468922e-01, 2.00481912e-01, 2.11683389e-01,
## 2.15689757e-01, 2.26237133e-01, 2.31306415e-01, 2.44143144e-01,
## 2.57715737e-01, 2.55589909e-01, 2.70879517e-01, 2.68917214e-01,
## 2.82980383e-01, 2.78483439e-01, 2.77093475e-01, 2.78892252e-01,
## 2.76766425e-01, 2.79628116e-01, 2.67200200e-01, 2.50684153e-01,
## 2.49703002e-01, 2.33023430e-01, 2.40218539e-01, 2.23048392e-01,
## 2.05796482e-01, 1.98764897e-01, 1.81267700e-01, 1.72519101e-01,
## 1.64506365e-01, 1.43902188e-01, 1.37933518e-01, 1.25260314e-01,
## 1.15530564e-01, 1.03020885e-01, 9.30458467e-02, 8.02908800e-02,
## 7.62845123e-02, 6.29572073e-02, 5.73155875e-02, 5.06110537e-02,
## 4.40700451e-02, 3.66296479e-02, 2.96798263e-02, 2.68998977e-02,
## 2.55099334e-02, 1.83148240e-02, 1.61889962e-02, 1.48807945e-02,
## 1.21826284e-02, 1.06291389e-02, 7.35863463e-03, 6.54100856e-03,
## 6.29572073e-03, 4.33341817e-03, 3.51579210e-03, 3.51579210e-03,
## 2.45287821e-03, 2.20759039e-03, 1.47172693e-03, 4.90575642e-04,
## 8.99388676e-04, 5.72338249e-04, 5.72338249e-04, 1.63525214e-04,
## 3.27050428e-04, 3.27050428e-04, 2.45287821e-04, 2.45287821e-04,
## 8.17626069e-05, 2.45287821e-04, 0.00000000e+00, 8.17626069e-05]), array([-6.44277106, -6.369366 , -6.29596095, -6.2225559 , -6.14915084,
## -6.07574579, -6.00234074, -5.92893568, -5.85553063, -5.78212558,
## -5.70872052, -5.63531547, -5.56191042, -5.48850536, -5.41510031,
## -5.34169526, -5.2682902 , -5.19488515, -5.1214801 , -5.04807504,
## -4.97466999, -4.90126494, -4.82785988, -4.75445483, -4.68104978,
## -4.60764472, -4.53423967, -4.46083461, -4.38742956, -4.31402451,
## -4.24061945, -4.1672144 , -4.09380935, -4.02040429, -3.94699924,
## -3.87359419, -3.80018913, -3.72678408, -3.65337903, -3.57997397,
## -3.50656892, -3.43316387, -3.35975881, -3.28635376, -3.21294871,
## -3.13954365, -3.0661386 , -2.99273355, -2.91932849, -2.84592344,
## -2.77251839, -2.69911333, -2.62570828, -2.55230323, -2.47889817,
## -2.40549312, -2.33208806, -2.25868301, -2.18527796, -2.1118729 ,
## -2.03846785, -1.9650628 , -1.89165774, -1.81825269, -1.74484764,
## -1.67144258, -1.59803753, -1.52463248, -1.45122742, -1.37782237,
## -1.30441732, -1.23101226, -1.15760721, -1.08420216, -1.0107971 ,
## -0.93739205, -0.863987 , -0.79058194, -0.71717689, -0.64377184,
## -0.57036678, -0.49696173, -0.42355668, -0.35015162, -0.27674657,
## -0.20334151, -0.12993646, -0.05653141, 0.01687365, 0.0902787 ,
## 0.16368375, 0.23708881, 0.31049386, 0.38389891, 0.45730397,
## 0.53070902, 0.60411407, 0.67751913, 0.75092418, 0.82432923,
## 0.89773429, 0.97113934, 1.04454439, 1.11794945, 1.1913545 ,
## 1.26475955, 1.33816461, 1.41156966, 1.48497471, 1.55837977,
## 1.63178482, 1.70518987, 1.77859493, 1.85199998, 1.92540504,
## 1.99881009, 2.07221514, 2.1456202 , 2.21902525, 2.2924303 ,
## 2.36583536, 2.43924041, 2.51264546, 2.58605052, 2.65945557,
## 2.73286062, 2.80626568, 2.87967073, 2.95307578, 3.02648084,
## 3.09988589, 3.17329094, 3.246696 , 3.32010105, 3.3935061 ,
## 3.46691116, 3.54031621, 3.61372126, 3.68712632, 3.76053137,
## 3.83393642, 3.90734148, 3.98074653, 4.05415159, 4.12755664,
## 4.20096169, 4.27436675, 4.3477718 , 4.42117685, 4.49458191,
## 4.56798696, 4.64139201, 4.71479707, 4.78820212, 4.86160717,
## 4.93501223, 5.00841728, 5.08182233, 5.15522739, 5.22863244,
## 5.30203749, 5.37544255, 5.4488476 , 5.52225265, 5.59565771,
## 5.66906276, 5.74246781, 5.81587287, 5.88927792, 5.96268297,
## 6.03608803, 6.10949308, 6.18289814, 6.25630319, 6.32970824,
## 6.4031133 , 6.47651835, 6.5499234 , 6.62332846, 6.69673351,
## 6.77013856, 6.84354362, 6.91694867, 6.99035372, 7.06375878,
## 7.13716383, 7.21056888, 7.28397394, 7.35737899, 7.43078404,
## 7.5041891 , 7.57759415, 7.6509992 , 7.72440426, 7.79780931,
## 7.87121436, 7.94461942, 8.01802447, 8.09142952, 8.16483458,
## 8.23823963]), <BarContainer object of 200 artists>)[1]
SIGNOROVITCH J E. WU E Q. YU A P. 等. Comparative effectiveness without head-to-head trials: a method for matching-adjusted indirect comparisons applied to psoriasis treatment with adalimumab or etanercept[J]. Pharmacoeconomics, 2010, 28: 935-945.