2020/11/11
为了便于计算，假设
之间相互独立，且
对
成立。令
由于指数分布是特殊的gamma分布，则由gamma分布的可加性知，
。
从而
的概率密度函数为
令
，则易得
的概率密度函数为
也可以通过定义求解
的分布函数，再求导得到其概率密度函数。从而
的期望为：
可以计算
的期望为：
从而，
的方差为：
令
表示样本均值，则样本均值的倒数为
，故样本均值的倒数的期望为
样本均值的倒数的方差为
整体思路就是，根据总体分布求样本和的分布，再求和的倒数的分布，计算出和的倒数的均值和方差，最后求样本均值的倒数的均值和方差。亦可用原概率密度函数直接对均值倒数求期望和方差，如有错误请指正~
期中考以后用蒙特卡洛模拟看看结果对不对
2020/12/3
回来更新啦~
一般情况：设样本的每个个体独立同分布服从于参数为
的指数分布，即
，也可写成
。则
。
令
，根据上述步骤，容易计算得到：
因此样本均值倒数的期望为：
方差为：
根据中心极限定理，通过蒙特卡洛模拟获取的期望和方差的样本数据，其样本均值的概率分布近似为正态分布，并且随着样本量(模拟次数)趋于无穷大，样本均值会收敛于期望值，因此我们可以通过样本均值估计理论值，并进行假设检验，验证理论值是否reasonable。原假设和备择假设分别为：
以及：
这里以
为例，模拟通过Python实现，代码如下：
import numpy as np
import random
from scipy.stats import norm
#生成服从指数分布的 x 总体
np.random.seed(0)
beta = 5
u = np.random.uniform(0, 1, 1000)
population = (-1/beta * np.log(1 - u)).tolist()
#初始化样本容量 n
n = 50
#通过 m 次抽样获得的目标值(样本均值倒数)的波动，用于检验方差
m = 100
#初始化模拟次数 M
M = 1000
#定义用于检验方差的函数，每调用一次这个函数，会返回一个方差估计值(基于 m 次抽取容量为 n 的样本数据)
def var_est(population, m, n):
data = [1/np.mean(random.sample(population, n)) for i in range(m)]
return np.var(data, ddof = 1)
#定义模拟函数，参数 target 用于选择要研究的统计量
def simulation(population, M = 1000, n = 50, target = 'mean'):
if target == 'mean':
theoretical_value = n * beta / (n-1)
data = [1/np.mean(random.sample(population, n)) for i in range(M)]
elif target == 'variance':
theoretical_value = n**2 * beta**2 / (n-1)**2 / (n-2)
data = [var_est(population, m, n) for i in range(M)]
return theoretical_value, data
#模拟两组数据，分别用于检验期望和方差
random.seed(1)
mean_theory, data_mean = simulation(population, M, n, target = 'mean')
variance_theory, data_var = simulation(population, M, n, target = 'variance')
#先看看期望和方差的理论值
mean_theory, variance_theory(5.1020408163265305, 0.5423087602387894)
#再看看样本数据的均值
np.mean(data_mean), np.mean(data_var)(5.080338272542856, 0.5403331158305243)
差异不大。
设定显著性水平为5%，分别计算两个假设检验问题的检验统计量和对应的p值：
test_mean = abs(np.mean(data_mean) - mean_theory) / np.std(data_mean, ddof = 1)
test_var = abs(np.mean(data_var) - variance_theory) / np.std(data_var, ddof = 1)
p_mean = 2 - 2 * norm.cdf(test_mean, loc = 0, scale = 1).round(3)
p_var = 2 - 2 * norm.cdf(test_var, loc = 0, scale = 1).round(3)
#输出 p 值
p_mean, p_var(0.976, 0.982)
两组假设检验的p值均大于0.05，因此没有充分证据表明我们的理论值是错误的~感兴趣的小伙伴可以调整样本容量，模拟次数和
值