方差分布 #
我们上次讲到了标准误和正态分布
。在其中，我们只看了样本平均值 ($\bar{X}$) 的分布，没有看样本方差 ($S^2$) 的分布。这次我们来看一下样本方差的分布有什么规律。
依然，我们以问题为导向。
全球人类的智商分布符合正态分布，
$\mathcal{N}(100, 15^2)$。我们从全球$80$亿人中随机抽取$100$个人，这$100$人智商的方差大于等于$250$的概率是多少？
首先，我们再次回到智商的正态分布图：
import numpy as np
import matplotlib.pyplot as plt
# Set the seed for reproducibility
np.random.seed(42)
# Parameters for the normal distribution
mean = 100 # Mean of the distribution
std_dev = 15 # Standard deviation of the distribution
size = 100000 # Number of data points
# Generate the data
data = np.random.normal(mean, std_dev, size)
population_mu = np.mean(data)
population_var = np.var(data)
# Create a histogram of the data
plt.figure(figsize=(6, 4))
plt.hist(data, bins="auto", density=False,
alpha=0.6, color='g', edgecolor='black')
# Add titles and labels
plt.title(f"Histogram of IQ Distribution, sample Size = {size}, "\
f"μ = {population_mu:.2f}, σ² = {population_var:.2f}")
plt.xlabel('IQ')
plt.ylabel('# of People')
# Show the plot
plt.show()
假设我们从母体中抽取 $M$ 个大小为 $N$ 的样本。也就是说，从母体中，每次抽取 $N$ 个人，重复 $M$ 次。每次是一个样本，其方差计算公式为：
$$ S^2 = \frac{1}{N-1} \sum_{i=1}^{N} (X_i - \bar{X})^2 $$
其中：
$N$ 是样本大小
$\bar{X}$ 是样本平均值
$X_i$ 代表样本中第 $i$ 个人的智商
下面，我们设 $M = 2000$，然后取不同的 $N$，画出 $S^2$ 的分布图。
def obtain_sample_vars(data, sample_size, num_samples):
"""sample from data, each time to get a sample of sample_size
and calcuate the sample's variance.
Continue for num_samples times.
"""
sample_vars = np.zeros(num_samples)
for i in range(num_samples):
sample = np.random.choice(data, size=sample_size, replace=False)
sample_vars[i] = np.var(sample, ddof=1)
return sample_vars
num_samples = 2000
sample_sizes = [10, 100, 1000, 10000]
sample_vars = {size: obtain_sample_vars(
data, size, num_samples) for size in sample_sizes}
# Plotting
plt.figure(figsize=(12, 10))
# Create subplots
for i, size in enumerate(sample_sizes):
sub_data = sample_vars[size]
mu = np.mean(sub_data)
sigma = np.std(sub_data)
plt.subplot(2, 2, i + 1) # 2x2 grid, subplot index starts at 1
plt.hist(sub_data, bins="auto", density=False,
alpha=0.7, color='steelblue', edgecolor='black')
# Plot vertical line at the mean
plt.axvline(mu, color='red', linestyle='--', linewidth=2)
plt.title(rf'Sample Variance with $N$ = {size}, '
rf'$\mathbb{{E}}[S^2] = {mu:.2f}$, '
rf'$S (S^2) = {sigma:.2f}$')
plt.xlabel('$S^2$')
plt.ylabel('# of Samples')
# Adjust layout
plt.tight_layout()
# Show the plot
plt.show()
print(f"母体平均值 (μ) = {population_mu:.2f}, 母体方差 (σ²) = {population_var:.2f}")
for sample_size in sample_sizes:
sub_data = sample_vars[sample_size]
sample_mean = np.mean(sub_data)
sample_var = np.var(sub_data)
sample_var_estimates = 2*population_var**2/(sample_size - 1)
print(f"N = {sample_size}, S² 的平均值 = {sample_mean:.2f}, "\
f"S² 方差 = {sample_var:.2f}, " \
f"2*σ⁴/(N-1) = {sample_var_estimates:.2f}")
$$\mathrm{E}[S^2] = \sigma^2$$
也就是说当 $N \ge 2$ 时，$S^2$ 的平均值和母体方差是 ($\sigma^2$) 基本上是一样的。
$S^2$ 的方差基本上符合这样的规律：
$$Var(S^2) = \frac{2 \cdot \sigma^4}{N-1}$$
另外，我们也看到，当 $N$ 比较大的时候，$S^2$ 的分布越来越趋近于正态分布：
$$S^2 \sim \mathcal{N} \left(\sigma^2, \frac{2 \cdot \sigma^4}{N-1} \right)$$
但是当 $N$ 小的时候，明显无法用正态分布来描述 $S^2$ 的分布：
num_samples = 2000
sample_sizes = [3, 5, 10, 20]
sample_vars = {size: obtain_sample_vars(
data, size, num_samples) for size in sample_sizes}
# Plotting
plt.figure(figsize=(12, 10))
# Create subplots
for i, size in enumerate(sample_sizes):
sub_data = sample_vars[size]
mu = np.mean(sub_data)
sigma = np.std(sub_data)
plt.subplot(2, 2, i + 1) # 2x2 grid, subplot index starts at 1
plt.hist(sub_data, bins="auto", density=False,
alpha=0.7, color='steelblue', edgecolor='black')
# Plot vertical line at the mean
plt.axvline(mu, color='red', linestyle='--', linewidth=2)
plt.title(rf'Sample Variance with $N$ = {size}, '
rf'$\mathbb{{E}}[S^2] = {mu:.2f}$, '
rf'$S (S^2) = {sigma:.2f}$')
plt.xlabel('$S^2$')
plt.ylabel('# of Samples')
# Adjust layout
plt.tight_layout()
# Show the plot
plt.show()
该怎么解决这个问题？
卡方分布 #
我们现在跳出来，暂时不去想现在这个问题，而来了解一下卡方分布 。
假设我们从标准正态分布 $\mathcal N(0, 1)$ 中随机抽取 $k$ 个数，$Z_1, ... , Z_k$ 然后求他们的平方和
$$Q = \sum_{i=1}^{k} Z_{i}^{2}$$
我们说 $Q$ 服从自由度为 $k$ 的卡方分布，$\chi_{k}^{2}$。
# Parameters for the normal distribution
mean = 0 # Mean of the distribution
std_dev = 1 # Standard deviation of the distribution
size = 100000 # Number of data points
# Generate the data
data = np.random.normal(mean, std_dev, size)
population_mu = np.mean(data)
population_var = np.var(data)
# Create a histogram of the data
plt.figure(figsize=(6, 4))
plt.hist(data, bins="auto", density=False,
alpha=0.6, color='g', edgecolor='black')
# Add titles and labels
plt.title(f"Standard Normal Distribution Simulation")
plt.xlabel('X')
plt.ylabel('Frequency')
# Show the plot
plt.show()
接下来，我们抽样。我们分别看一下 $k = 1, 2, 3, 5, 10, 20$ 的情况。每种情况我们抽 $2000$ 个样本。
def obtain_q(data, sample_size, num_samples):
"""sample from data, each time to get a sample of sample_size (k)
and calcuate Q.
Continue for num_samples times.
"""
qs = np.zeros(num_samples)
for i in range(num_samples):
sample = np.random.choice(data, size=sample_size, replace=False)
qs[i] = sum(k**2 for k in sample)
return qs
num_samples = 2000
sample_sizes = [1, 2, 3, 5, 10, 20]
qs = {size: obtain_q(
data, size, num_samples) for size in sample_sizes}
from scipy.stats import chi2
# Plotting
plt.figure(figsize=(12, 10))
# Create subplots
for i, size in enumerate(sample_sizes):
sub_data = qs[size]
plt.subplot(3, 2, i + 1) # 2x2 grid, subplot index starts at 1
plt.hist(sub_data, bins="auto", density=True,
alpha=0.7, color='steelblue', edgecolor='black')
# Plot the theoretical chi-square distribution line
x = np.linspace(0, np.max(sub_data), 1000) # Define the x-axis range
plt.plot(x, chi2.pdf(x, df=size),
color="green", lw=2,
label=f'Chi-Square PDF (df={size})') # Plot PDF line
plt.title(rf'Q with $k$ = {size}')
plt.xlabel('$Q$')
plt.ylabel('Density')
plt.legend()
# Adjust layout
plt.tight_layout()
# Show the plot
plt.show()
科学家根据标准正态分布的特性推导出了
卡方分布的概率密度函数:
$$f(x; k) = \frac{1}{2^{k/2} \Gamma(k/2)} x^{(k/2) - 1} e^{-x/2}, \quad \text{for } x \geq 0$$
其中 $x$ 就是 $X$ 轴上的任意数字。$k$ 是自由度，也就是我们上面的例子中每次随机抽取几个数。
$\Gamma$ 函数的定义：
$$\Gamma(n) = \int_0^{\infty} t^{n-1} e^{-t} \, dt$$
上面我把精确的卡方分布密度函数也画出来了。
上面是从 $\mathcal N(0, 1)$ 随机抽样，那如果换成更普遍的情况呢？
假设我们从正态分布 $\mathcal N(\mu, \sigma^2)$ 中随机抽取 $k$ 个数，$X_1, ... , X_k$。
$\sum_{X_i}^{2}$ 确实不是卡方分布，但是
$$Z_i = \frac{X_i - \mu}{\sigma}$$
$$Q = \sum_{i=1}^{k} Z_{i}^{2}$$
$Q$ 依旧符合卡方分布。
继续假设，如果我们不知道母体的平均值 ($\mu$) 但是知道 $\sigma$，也就是说我们只能用样本平均值 $\bar{X}$ 来获取 $Z$ 的近似值:
$$T_i = \frac{X_i - \bar{X}}{\sigma}$$
统计学家证明
$$Q = \sum_{i=1}^{k} T_{i}^{2} = \sum_{i=1}^{k} \frac{(X_i - \bar{X})^2}{\sigma^2}$$
的分布依旧是一个卡方分布，但是自由度要减去一：
$$Q \sim \chi_{k-1}^{2}$$
因为
$$S^2 = \frac{\sum_{i=1}^{k} (X_i - \bar{X})^2}{k-1}$$
所以
$$Q = \frac{S^2 (k-1)}{\sigma^2}$$
所以，我们最后的结论是：
$$\frac{S^2 (k-1)}{\sigma^2} \sim \chi_{k-1}^{2}$$
也就是说 $\frac{S^2 (k-1)}{\sigma^2}$ 服从自由度为 $k-1$ 的卡方分布。其中
$k$ 是样本量
$S^2$ 是样本方差
$\sigma$ 是母体标准差
解决问题 #
现在我们回到最初的问题：
全球人类的智商分布符合正态分布，
$\mathcal{N}(100, 15^2)$。我们从全球$80$亿人中随机抽取$100$个人，这$100$人智商的方差大于等于$250$的概率是多少？
我们看到 $k=100$。$S^2 \ge 250 $ 的概率也就是 在一个自由度为 $99$ 的卡方分布中，$x \ge \frac{S^2 (k-1)}{\sigma^2} = \frac{250 \times 99}{15^2}$ 的概率。
df = 99
threshold = (250*df)/(15**2)
prob = 1 - chi2.cdf(threshold, df)
x = np.linspace(0, 300, 1000)
y = chi2.pdf(x, df)
plt.figure(figsize = (10, 6))
plt.plot(x, y, label=f'Chi-Square PDF (df = {df})', color = "steelblue")
plt.fill_between(x, y,
# where=np.logical_or(x >= threshold, x <= 80),
where=(x >= threshold),
color='orange',
alpha = 0.5,
label=f'P(x >= {threshold:.2f})')
plt.title('Chi-Square Distribution with Highlighted Region')
plt.xlabel('$\chi^2$')
plt.ylabel('Probability Density')
plt.axvline(threshold, color="red",
linestyle='--',
label=f'Threshold = {threshold:.2f}')
plt.legend()
plt.grid(True)
plt.show()
prob
总结 #
我们从一个均值为 $\mu$、方差为 $\sigma^2$ 的正态分布 $\mathcal N(\mu, \sigma^2)$ 中随机取 $k$ 个独立同分布
(i.i.d) 的随机变量 $X_1, X_2, ... , X_k$。这些随机变量构成一个样本。该样本的方差为 $S^2$。取无数个样本。
对样本的方差进行标准化：
$$Q = \frac{(k-1) S^2}{\sigma^2}$$
$Q$ 服从自由度为 $k-1$ 的卡方分布：
$$Q \sim \chi_{k-1}^{2}$$
dfs = [1, 2, 3, 4, 6, 9]
x = np.linspace(0, 8, 1000) # Define the x-axis range
plt.figure(figsize=(10, 6))
# Colors for each chi-square distribution
colors = ['olive', 'limegreen', 'cyan', 'blue', 'purple', 'red']
for df, color in zip(dfs, colors):
plt.plot(x, chi2.pdf(x, df=df), color=color, lw=2,
label=rf'$k={df}$') # Plot PDF line
plt.xlabel(r'$\chi^2$', fontsize = 14)
plt.ylabel(r'Probability Density', fontsize = 14)
plt.title(r'$\chi^2_k$ Distributions for Different Degrees of Freedom', fontsize = 16)
plt.legend(title = 'Degrees of Freedome', fontsize = 12)
plt.grid(True)
plt.ylim(0, 0.5) # Restrict y-axis limits
plt.show()
最后一次修改于 2024-12-07