做一晚水泥工歌曲网站,网站建设对策,wordpress地址和站点地址,网站开发常用图标系列文章目录
统计计算一|非线性方程的求解 统计计算二|EM算法#xff08;Expectation-Maximization Algorithm#xff0c;期望最大化算法#xff09; 统计计算三|Cases for EM 文章目录 系列文章目录一、基本概念#xff08;一#xff09;估算 π \pi π#xff08;二Expectation-Maximization Algorithm期望最大化算法 统计计算三|Cases for EM 文章目录 系列文章目录一、基本概念一估算 π \pi π二求积分三使用步骤 二、变换采样一基本概念二相关定理三示例 三、逆变换采样一相关定理二逆变换法三采样步骤四示例 四、接受拒绝采样一基本概念二工作原理三采样步骤四示例 五、重要性采样一基本概念二蒙特卡洛估计三示例 一、基本概念
蒙特卡洛方法为了解决某确定性问题把它变成一个概率模型的求解问题然后产生符合模型的大量随机数对产生的随机数进行分析从而求解问题的方法又称为随机模拟方法。
随机数设 X X X 是具有分布函数 F ( x ) F(x) F(x) 的随机变量从分布 F ( x ) F(x) F(x) 中随机抽样得到的序列 { x i , i 1 , 2 , . . . } \{x_i, i 1, 2, ...\} {xi,i1,2,...} 称为该分布的随机数序列 x i x_i xi 称为分布 F ( x ) F(x) F(x) 的随机数。
一估算 π \pi π
向正方形 D { ( x , y ) : x ∈ [ 0 , 1 ] , y ∈ [ 0 , 1 ] } D\{(x,y):x\in[0,1],y\in[0,1]\} D{(x,y):x∈[0,1],y∈[0,1]}内随机等可能投点落入四分之一圆 C { ( x , y ) : x 2 y 2 ≤ 1 , x 0 , y 0 } C\{(x,y):x^2y^2\leq 1,x0,y0\} C{(x,y):x2y2≤1,x0,y0}的概率为面积之比 p π 4 p\frac{\pi}{4} p4π。如果独立重复地投了 n n n个点落入 C C C中的点的个数为 ξ \xi ξ则有 ξ n ≈ π 4 , π ≈ π ^ 4 ξ n \frac{\xi}{n}\approx \frac{\pi}{4},\ \pi\approx \hat{\pi}\frac{4\xi}{n} nξ≈4π, π≈π^n4ξ 由于 ξ \xi ξ服从 B i n o m i a l ( n , π 4 ) Binomial(n,\frac{\pi}{4}) Binomial(n,4π)分布有 V a r ( π ^ ) π ( 4 − π ) 16 n Var(\hat{\pi})\frac{\pi(4-\pi)}{16n} Var(π^)16nπ(4−π) 由中心极限定理 π ^ \hat{\pi} π^近似服从 N ( π , π ( 4 − π ) 16 n ) N(\pi,\frac{\pi(4-\pi)}{16n}) N(π,16nπ(4−π))分布所以随机模拟误差的幅度大约在 ± 2 π ( 4 − π ) 16 n \pm 2\sqrt{\frac{\pi(4-\pi)}{16n}} ±216nπ(4−π) 随机模拟误差95%以上落入此区间
二求积分
将积分转化为期望来计算对于 Q ∫ a b h ( x ) d x Q\int_a^bh(x)dx Q∫abh(x)dx取 U ∼ U ( a , b ) U\sim U(a,b) U∼U(a,b)有 Q ( b − a ) ∫ a b h ( u ) 1 b − a d u ( b − a ) E [ h ( U ) ] Q(b-a)\int_a^bh(u)\frac{1}{b-a}du(b-a)E[h(U)] Q(b−a)∫abh(u)b−a1du(b−a)E[h(U)] 若取 { U i , i 1 , . . . , N } \{U_i,i1,...,N\} {Ui,i1,...,N}独立同 U ( a , b ) U(a,b) U(a,b)分布并设 Y i h ( U i ) , i 1 , 2 , . . . , N Y_ih(U_i),i1,2,...,N Yih(Ui),i1,2,...,N是 i i d iid iid随机变量列由强大数律 Y ˉ 1 N ∑ i 1 N h ( U i ) → E h ( U ) Q b − a , a . s . ( N → ∞ ) \bar{Y}\frac{1}{N}\sum_{i1}^Nh(U_i)\rightarrow Eh(U)\frac{Q}{b-a},\ a.s.(N\rightarrow ∞) YˉN1i1∑Nh(Ui)→Eh(U)b−aQ, a.s.(N→∞) 于是有 Q ^ ( b − a ) Y ˉ b − a N ∑ i 1 N h ( U i ) \hat{Q}(b-a)\bar{Y}\frac{b-a}{N}\sum_{i1}^Nh(U_i) Q^(b−a)YˉNb−ai1∑Nh(Ui)
由中心极限定理 N ( Q ^ − Q ) → d N ( 0 , ( b − a ) 2 V a r ( h ( U ) ) ) \sqrt{N}(\hat{Q}-Q)\xrightarrow{d} N(0,(b-a)^2Var(h(U))) N (Q^−Q)d N(0,(b−a)2Var(h(U))) V a r [ h ( U ) ] ∫ a b [ h ( u ) − E h ( U ) ] 2 1 b − a d u Var[h(U)]\int_a^b[h(u)-Eh(U)]^2\frac{1}{b-a}du Var[h(U)]∫ab[h(u)−Eh(U)]2b−a1du V a r [ ( h ( U ) ] Var[(h(U)] Var[(h(U)]可以用模拟样本 { Y i h ( U i ) } \{Y_ih(U_i)\} {Yih(Ui)}估计为 V a r ( h ( U ) ) ≈ 1 N ∑ i 1 N ( Y i − Y ˉ ) 2 Var(h(U))\approx\frac{1}{N}\sum_{i1}^N(Y_i-\bar{Y})^2 Var(h(U))≈N1i1∑N(Yi−Yˉ)2
三使用步骤
蒙特卡洛方法的理论基础是大数定律。样本数量越多则随机数的平均值就越接近期望也就是要计算的真实值。
将实际问题转化为求期望并定义要采样的随机变量计算机模拟采样过程处理产生的随机数得到期望
二、变换采样
一基本概念
如果随机变量 η η η 不容易抽样但是存在另一个容易抽样的随机变量 ξ ξ ξ 和随机变量 η η η 间具有一一对应关系即 η h ( ξ ) η h(ξ) ηh(ξ) 或 ξ h − 1 ( η ) ξ h^{−1}(η) ξh−1(η)同分布。那么可以先产生随机变量 ξ ξ ξ再由函数关系 h ( ⋅ ) h(·) h(⋅) 得到随机变量 η η η这种产生随机数的方法称为变换抽样法。
二相关定理
设随机变量 ξ \xi ξ具有概率密度函数 f ( x ) f(x) f(x)另有一函数 h ( ⋅ ) h(·) h(⋅)严格单调其反函数记为 h − 1 ( ⋅ ) h^{-1}(·) h−1(⋅)且导函数存在则 η h ( ξ ) \etah(\xi) ηh(ξ)是随机变量 ξ \xi ξ的函数其概率密度函数为 p ( z ) f ( h − 1 ( z ) ) ⋅ ∣ { h − 1 ( z ) } ′ ∣ p(z)f(h^{-1}(z))·|\{h^{-1}(z)\}| p(z)f(h−1(z))⋅∣{h−1(z)}′∣ 证明 三示例 用变换抽样法产生分布为 N ( µ , σ 2 ) N(µ, σ2) N(µ,σ2) 的随机数。 用变换抽样法产生分布为 G a m m a ( 1 / 2 , 3 ) Gamma(1/2, 3) Gamma(1/2,3) 的随机数。
三、逆变换采样
一相关定理
假设 X X X为一个连续随机变量其累计分布函数为 F X F_X FX此时可证明随机变量 Y F X ( X ) YF_X(X) YFX(X)服从区间 [ 0 , 1 ] [0,1] [0,1]上的均匀分。逆变换采样就是将上述过程反过来进行。
设连续型随机变量 η \eta η的分布函数 F ( x ) F(x) F(x)是连续且严格单调上升的分布函数其反函数存在且记为 F − 1 ( x ) F^{-1}(x) F−1(x)。则有
随机变量 F ( η ) F(\eta) F(η)服从 ( 0 , 1 ) (0,1) (0,1)上的均匀分布即 F ( η ) ∼ U ( 0 , 1 ) F(\eta)\sim U(0,1) F(η)∼U(0,1)对于随机变量 U ∼ U ( 0 , 1 ) U\sim U(0,1) U∼U(0,1) F − 1 ( U ) F^{-1}(U) F−1(U)的分布函数为 F ( x ) F(x) F(x) 证明 二逆变换法
逆变换法当随机变量 η \eta η的分布函数 F ( x ) F(x) F(x)的反函数存在且容易计算时可通过产生均匀分布的随机数来产生 η \eta η的随机数序列 { η i , i 1 , 2 , . . . } \{\eta_i,i1,2,...\} {ηi,i1,2,...}。这种产生非均匀分布随机数的方法称为逆变换法或反函数法。
三采样步骤
产生 U ( 0 , 1 ) U(0,1) U(0,1)的随机数序列 { u i , i 1 , 2 , . . . } \{u_i,i1,2,...\} {ui,i1,2,...} η \eta η的随机数序列为 η i F − 1 ( u i ) , i 1 , 2 , . . . \eta_iF^{-1}(u_i),i1,2,... ηiF−1(ui),i1,2,...
四示例
产生概率密度函数为 f(x) 的随机数其中 f ( x ) { x σ 2 e − x 2 2 σ 2 , x 0 0 , z ≤ 0 f(x) \begin{cases} \frac{x}{\sigma^2}e^{-\frac{x^2}{2\sigma^2}}, \text{$x0$} \\ 0, \text{$z\leq0$} \end{cases} f(x){σ2xe−2σ2x2,0,x0z≤0 产生分布函数为 F ( x ) F(x) F(x) 的随机数 η η η其中 F ( x ) x 2 x 2 , 0 ≤ x ≤ 1 F(x)\frac{x^2x}{2},0\leq x\leq 1 F(x)2x2x,0≤x≤1
四、接受拒绝采样
一基本概念
拒绝抽样是基于以下观察而提出的要在一维中抽样一个随机变量可以对二维笛卡尔图进行均匀随机抽样并将样本保留在其密度函数图形下的区域中。 想象将一个随机变量的密度函数绘制在一个大矩形板上并向其投掷飞镖。假设这些飞镖在整个板上均匀分布。现在移除所有落在曲线下方以外区域的飞镖。剩下的飞镖将在曲线下方的区域内均匀分布并且这些飞镖的 x 坐标将按照随机变量的密度分布。这是因为在曲线最高的地方也就是概率密度最大的地方飞镖着陆的空间最多。 拒绝抽样的一般形式假设板子的形状不一定是矩形而是根据某个提议分布的密度来确定该分布不一定归一化为 1。通常情况下将其视为某个已知的分布的倍数。提议分布中的每个点至少与想要抽样的分布一样高以便前者完全包围后者。否则想要抽样的曲线区域中的某些部分可能永远无法到达。
二工作原理
拒绝抽样的工作原理
从提议分布中在 x 轴上抽样一个点。在该 x 位置上画一条竖直线直到提议分布的概率密度函数的 y值。在这条线上从 0 到提议分布的概率密度函数的 y 值之间均匀抽样。如果抽样值大于该竖直线上所需分布的密度函数值则拒绝该 x 值并返回第 1 步否则该 x 值就是所需分布的一个样本。
拒绝抽样算法可以用于从任何曲线下方进行抽样无论函数是否积分为 1。事实上通过常数缩放函数对抽样的 x 位置没有影响。因此该算法可以用于从归一化常数未知的分布中进行抽样。
三采样步骤
提案分布 g g g为了从密度为 f f f的分布 X X X中获取样本利用了容易采样的密度函数为 g g g的分布 Y Y Y g g g就是提案分布
设 M M M为似然比 f ( x ) / g ( x ) f(x)/g(x) f(x)/g(x)的上界即一个常数满足 1 ≤ M ∞ 1\leq M∞ 1≤M∞。也就是说 M M M必须满足 f ( x ) ≤ M g ( x ) f(x)\leq Mg(x) f(x)≤Mg(x)对任意 x x x都成立因此 Y Y Y分布的支撑要包含 X X X的支撑
从分布 Y Y Y获取样本 y ∼ g y\sim g y∼g并从 U n i f ( 0 , 1 ) Unif(0,1) Unif(0,1)单位区间上的均匀分布获取样本 u u u检查是否 u f ( y ) / M g ( y ) uf(y)/Mg(y) uf(y)/Mg(y).( M ≥ 1 M\geq 1 M≥1) 成立则接受 y y y作为从 f f f中抽取的样本不成立则拒绝 y y y的值并重新获取样本
保留样本不大于值 y 的概率为
其中 P [ U ≤ f ( Y ) M g ( Y ) ] 1 M P[U\leq \frac{f(Y)}{Mg(Y)}]\frac{1}{M} P[U≤Mg(Y)f(Y)]M1为接受率接受率越大采样效率就越高。接受拒绝算法平均需要 M 次迭代才能获得样本并且M 越小越好即包络线越贴近目标分布越好可设 M max x f ( x ) g ( x ) M\max_x \frac{f(x)}{g(x)} Mxmaxg(x)f(x) 四示例
试用接受拒绝采样产生服从均值为 0方差为 1 的半正态分布的随机数 η η η该分布的概率密度函数为 p ( z ) { 2 / π e − z 2 2 , z ≥ 0 0 , z 0 p(z) \begin{cases} \sqrt{2/\pi}e^{-\frac{z^2}{2}}, \text{$z\geq0$} \\ 0, \text{$z0$} \end{cases} p(z){2/π e−2z2,0,z≥0z0 五、重要性采样
一基本概念
接受拒绝采样完美的解决了累积分布函数不可求时的采样问题。但是接受拒绝采样非常依赖于提议分布的选择如果提议分布选择的不好可能采样时间很长却获得很少满足分布的粒子。而重要性采样就解决了这一问题。
重要性采样是使用蒙特卡洛方法估算积分 (期望) 时提高对积分计算重要区域的抽样从而达到减少方差的目的。 μ E X ∼ f ( h ( X ) ) ∫ h ( x ) f ( x ) d x ∫ h ( x ) f ( x ) g ( x ) g ( x ) d x \muE_{X\sim f}(h(X))\int h(x)f(x)dx\int h(x)\frac{f(x)}{g(x)}g(x)dx μEX∼f(h(X))∫h(x)f(x)dx∫h(x)g(x)f(x)g(x)dx
或若 ∫ f ( x ) ≠ 1 \int f(x)\neq 1 ∫f(x)1, 也就是只知道分布成比例于某个函数差一个归一化常数则 μ E X ∝ f ( h ( X ) ) ∫ h ( x ) f ( x ) ∫ f ( x ) d x d x ∫ h ( x ) f ( x ) g ( x ) g ( x ) d x ∫ f ( x ) g ( x ) g ( x ) d x \muE_{X∝f}(h(X))\int h(x)\frac{f(x)}{\int f(x)dx}dx\frac{\int h(x)\frac{f(x)}{g(x)}g(x)dx}{\int \frac{f(x)}{g(x)}g(x)dx} μEX∝f(h(X))∫h(x)∫f(x)dxf(x)dx∫g(x)f(x)g(x)dx∫h(x)g(x)f(x)g(x)dx
二蒙特卡洛估计
上式建议用来估计 E h ( X ) Eh(X) Eh(X) 的一种 Monte Carlo 方法 从 g g g中抽取独立同分布的样本 X 1 , . . . , X n X_1,...,X_n X1,...,Xn并采用估计 μ ^ I S ∗ 1 n ∑ i h ( x i ) f ( x i ) g ( x i ) → E X ∼ g ( h ( x ) f ( x ) g ( x ) ) \hat{\mu}^*_{IS}\frac{1}{n}\sum_ih(x_i)\frac{f(x_i)}{g(x_i)}\rightarrow E_{X\sim g}(h(x)\frac{f(x)}{g(x)}) μ^IS∗n1i∑h(xi)g(xi)f(xi)→EX∼g(h(x)g(x)f(x)) 可以写成 w ∗ ( X i ) f ( X i ) / g ( X i ) w^*(X_i)f(X_i)/g(X_i) w∗(Xi)f(Xi)/g(Xi)是未标准化权重称为重要性比率 μ ^ I S ∗ 1 n ∑ i h ( X i ) w ∗ ( X i ) \hat{\mu}^*_{IS}\frac{1}{n}\sum_ih(X_i)w^*(X_i) μ^IS∗n1i∑h(Xi)w∗(Xi) 若是差一个比例常数的 f则 w ∗ ( X i ) w ∗ ( X i ) / ∑ i 1 n w ∗ ( X i ) w^*(X_i)w^*(X_i)/\sum_{i1}^nw^*(X_i) w∗(Xi)w∗(Xi)/∑i1nw∗(Xi)是标准化权重 μ ^ I S 1 n ∑ i h ( X i ) w ( X i ) \hat{\mu}_{IS}\frac{1}{n}\sum_ih(X_i)w(X_i) μ^ISn1i∑h(Xi)w(Xi) 估计值的方差是 V a r ( μ ^ ) 1 n V a r h ( X ) f ( X ) g ( X ) 1 n ∫ ( h ( x ) f ( x ) g ( x ) − μ 2 ) 2 g ( x ) d x 1 n { ∫ ( h 2 ( x ) f 2 ( x ) g ( x ) d x − μ 2 } \begin{aligned} Var(\hat{\mu})\frac{1}{n}Var\frac{h(X)f(X)}{g(X)} \\ \frac{1}{n}\int (\frac{h(x)f(x)}{g(x)}-\mu^2)^2g(x)dx\\ \frac{1}{n}\{ \int(\frac{h^2(x)f^2(x)}{g(x)}dx-\mu^2\}\\ \end{aligned} Var(μ^)n1Varg(X)h(X)f(X)n1∫(g(x)h(x)f(x)−μ2)2g(x)dxn1{∫(g(x)h2(x)f2(x)dx−μ2} 如果 h 2 ( x ) f 2 ( x ) / g 2 ( x ) μ 2 {h^2(x)f^2(x)}/{g^2(x)}\mu^2 h2(x)f2(x)/g2(x)μ2就有 V a r ( μ ^ ) 0 Var(\hat{\mu})0 Var(μ^)0即方差达到最小。此时 g ( x ) ∣ h ( x ) ∣ f ( x ) μ ∣ h ( x ) ∣ f ( x ) ∫ ∣ h ( x ) ∣ f ( x ) d x g(x)\frac{|h(x)|f(x)}{\mu}\frac{|h(x)|f(x)}{\int |h(x)|f(x)dx} g(x)μ∣h(x)∣f(x)∫∣h(x)∣f(x)dx∣h(x)∣f(x) 密度函数 g ( x ) g(x) g(x) 的最佳选择就是和被积函数 ∣ h ( x ) ∣ f ( x ) |h(x)|f(x) ∣h(x)∣f(x) 具有相同的形状。对积分值贡献越大的区域希望以较大的概率抽取到随机数。 在实际中 µ µ µ 是未知量因此无法选取 g ( x ) g(x) g(x)使得 V a r ( µ ^ ) 0 Var(\hat{µ}) 0 Var(µ^)0。通常情况下我们会选取一个形状接近 ∣ h ( x ) ∣ f ( x ) |h(x)|f(x) ∣h(x)∣f(x) 的函数作为 g ( x ) g(x) g(x)。 μ ∫ h ( x ) f ( x ) d x ∫ h ( x ) f ( x ) g ( x ) g ( x ) d x \mu\int h(x)f(x)dx\int h(x)\frac{f(x)}{g(x)}g(x)dx μ∫h(x)f(x)dx∫h(x)g(x)f(x)g(x)dx
三示例
例用重要抽样法估计 μ ∫ 0 1 x f ( x ) d x ∫ 0 1 e x d x \mu\int_0^1xf(x)dx\int_0^1e^xdx μ∫01xf(x)dx∫01exdx的估计值
参考 MCMC入门一蒙特卡罗方法与拒绝-接受采样
