高斯分布

正态分布,又称Gauss分布,其概率密度函数入下图所示

在这里插入图片描述
正态分布 N ( μ , σ ) N(\mu, \sigma) N(μ,σ)受到期望 μ \mu μ和方差 σ 2 \sigma^2 σ2的调控,其概率密度函数为

1 2 π σ 2 exp ⁡ [ − ( x − μ ) 2 2 σ 2 ] \frac{1}{\sqrt{2\pi\sigma^2}}\exp[-\frac{(x-\mu)^2}{2\sigma^2}] 2πσ2 1exp[2σ2(xμ)2]

μ = 0 \mu=0 μ=0 σ = 1 \sigma=1 σ=1时,为标准正态分布 N ( 0 , 1 ) N(0,1) N(0,1),对应概率分布函数为 Φ ( x ) = 1 2 π exp ⁡ [ − x 2 2 ] \Phi(x)=\frac{1}{\sqrt{2\pi}}\exp[-\frac{x^2}{2}] Φ(x)=2π 1exp[2x2]

逆高斯分布简介

在布朗运动中,高斯分布描述的是某一固定时刻距离的分布;而逆高斯分布则是达到固定距离所需时间的分布。

所以逆高斯分布的逆取自于物理意义,而非在表达式上有什么“逆”的地方,其PDF为

p ( x ) = λ 2 π x 3 exp ⁡ [ − λ ( x − μ ) 2 2 μ 2 x ] p(x)=\sqrt{\frac{\lambda}{2\pi x^3}}\exp[-\frac{\lambda(x-\mu)^2}{2\mu^2x}] p(x)=2πx3λ exp[2μ2xλ(xμ)2]

scipy.stats中提供了invgauss函数,但其默认 λ = 1 \lambda=1 λ=1。当 λ \lambda λ趋近于无穷大时,逆高斯分布趋向于高斯分布。

特别地,当 μ = λ = 1 \mu=\lambda=1 μ=λ=1时,逆高斯分布又被称为Wald分布,其概率密度函数表达式为

p ( x ) = 1 2 π x 3 exp ⁡ [ − ( x − 1 ) 2 2 x ] p(x)=\frac{1}{\sqrt{2\pi x^3}}\exp[-\frac{(x-1)^2}{2x}] p(x)=2πx3 1exp[2x(x1)2]

scipy.stats中也封装了wald函数,即wald

通过高斯分布构造逆高斯分布

通过正态分布和均匀分布的随机数,可以生成逆高斯分布的随机数,设 ν \nu ν服从标准正态分布;z服从标准均匀分布,则令

x = μ + μ 2 ν 2 2 λ − μ 2 λ 4 μ λ ν 2 + μ 2 ν 4 x=\mu+\frac{\mu^2\nu^2}{2\lambda}-\frac{\mu}{2\lambda}\sqrt{4\mu\lambda\nu^2+\mu^2\nu^4} x=μ+2λμ2ν22λμ4μλν2+μ2ν4

z ⩽ μ μ + x z\leqslant\frac{\mu}{\mu+x} zμ+xμ,则返回 x x x,否则返回 μ 2 x \frac{\mu^2}{x} xμ2

考虑到invgauss默认 λ = 1 \lambda=1 λ=1,所以 x x x的生成式改为 x = μ + μ 2 ν 2 2 − μ 2 4 μ ν 2 + μ 2 ν 4 x=\mu+\frac{\mu^2\nu^2}{2}-\frac{\mu}{2}\sqrt{4\mu\nu^2+\mu^2\nu^4} x=μ+2μ2ν22μ4μν2+μ2ν4

import numpy as np
from scipy.stats import norm, invgauss, uniform
import matplotlib.pyplot as plt

nu = norm.rvs(size=10000)
mu = 1      # 设mu=1
n2, m2 = nu**2, mu**2
x = mu + m2*n2/2-mu/2*np.sqrt(4*mu*n2+m2*nu**4)

z = uniform.rvs(size=10000)
flag = z<(mu/(mu+x))

rs = flag*x + (1-flag)*m2/x


rv = invgauss(mu)
st, ed = rv.interval(0.995)
xs = np.linspace(st, ed, 200)

def drawPDF(rs, xs, ys):
    plt.figure(figsize=(5,3))
    plt.hist(rs, density=True, bins=100, alpha=0.8)
    plt.plot(xs, ys)
    plt.tight_layout()
    plt.show()

drawPDF(rs, xs, rv.pdf(xs))

结果如下

Logo

腾讯云面向开发者汇聚海量精品云计算使用和开发经验,营造开放的云计算技术生态圈。

更多推荐