MCS：连续随机变量——Chi-Square分布

Chi-Square 分布

Chi-Square分布是统计分析上最常用的分布之一，用于检验变量的方差的变化，Chi-Square变量用${\mathcal{X}}^2$表示，自由度参数用$k$表示（正整数）通常$k<=100$,令$w={\mathcal{X}}^2$,概率密度函数为：

$$f(w)=\frac{w^{(k/2-1)}{e}^{-w/2}}{2^{k/2}\Gamma (k/2)}，w>=0$$

$\Gamma (k)$:Gamma函数:
$$\Gamma (k)={\int }_0^{\infty }{t}^{k-1}{e}^{-t}dt，k>0$$
如果$k$为正整数(Gamma = Erlang):
$$\Gamma (k)=(k-1)!$$

$w$(${\mathcal{X}}^2$)的均值和方差为：$E(w)=k，V(w)=2k$

$$E({\mathcal{X}}^2)=k$$
$$V({\mathcal{X}}^2)=2k$$

自由度为$k$的卡方变量与标准正态变量的关系：

$${\mathcal{X}}^2=z_1^2+z_2^2+...+z_k^2$$

当$k>30$，根据中心极限定理，chi-square的概率密度将呈现出正态分布，${\mathcal{X}}^2\sim N(k,2k)$,根据这个关系，$\alpha$-percent chi-square近似值：

$${\mathcal{X}}_{\alpha }^2\approx k+z_{\alpha }\sqrt{2k}$$

• $z$ : 标准正态变量， $P(z>z_{\alpha })=\alpha$, $P({\mathcal{X}}^2>{\mathcal{X}}_{\alpha }^2)\approx \alpha$
• $k$ : 均值，$2k$ : 方差

同Gamma分布的关系：

当Gamma分布的参数$\theta =2,k^{\prime }=k/2$时，他们具有相同的形状。

生成随机的Chi-Square变量

根据卡方变量跟标准正态变量之间的关系，可以通过K个标准的正态变量来生成自由度为k的卡方变量，计算步骤如下：

1. ${\mathcal{X}}^2=0$
2. For i in range(1, k+1).
3. Generate a standard normal variate z.
4. ${\mathcal{X}}^2={\mathcal{X}}^2+z^2$.
5. Next i.
6. Return ${\mathcal{X}}^2$.

根据Chi-Square变量跟Gamma变量的关系，也可以通过Gamma分布生成随机Chi-Square变量。

例：假设一个随机卡方变量自由度为$k=3$,模拟生成服从这样分布的随机变量，过程如下：

1. Generate three random standard normal variates : $z_1,z_2,z_3=0.47,-0.81,1.04$
2. ${\mathcal{X}}^2=0.47^2+-0.81^2+1.04^2=1.9586$
3. Return ${\mathcal{X}}^2=1.9586$

例：假设一个随机卡方变量的自由度为239，模拟生成一个这样的随机变量：

注：根据Chi-Square分布的特点，当自由度较高时(k > 30),Chi-Square分布将呈现出正态分布的形状。(N(k, 2k)),根据这个性质可以通过标准正态分布生成Chi-Square变量。

1. Generate a random standard normal $z\sim N(0,1)，z=1.34$
2. ${\mathcal{X}}^2\approx int(k+z\times \sqrt{2k}+0.5)=269$
3. Return ${\mathcal{X}}^2=269$

import numpy as np
import matplotlib.pyplot as plt

def generate_chi_square_var(k=1):
    if k <= 30:
        z = np.random.normal(0, 1, size=(k))**2
        x = z.sum()
    else:
        z = np.random.normal(0, 1)
        x = int(k + z * np.sqrt(2*k) + 0.5)
    return x

$k=1、10、30、50$

x_k1 = [generate_chi_square_var() for i in range(10000)]
x_k10 = [generate_chi_square_var(10) for i in range(10000)]
x_k30 = [generate_chi_square_var(30) for i in range(10000)]
x_k50 = [generate_chi_square_var(50) for i in range(10000)]
plt.figure(figsize=(12, 5))
plt.hist(x_k1, label="$k=1$")
plt.hist(x_k10, label="$k=10$", alpha=0.5)
plt.hist(x_k30, label="$k=30$", alpha=0.5)
_=plt.hist(x_k50, label="$k=50$", alpha=0.5)
plt.legend()
plt.grid()