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，因此没有充分证据表明我们的理论值是错误的~感兴趣的小伙伴可以调整样本容量，模拟次数和

值