26 概率统计应用
26.1 概率论例子
26.1.1 取帽子问题
26.1.1.1 问题介绍与推导
设\(n\)个人都戴了帽子参加聚会， 聚会时摘下帽子弄乱了， 散会时随机抽取一顶。 问：所有人都没有取到自己的帽子的概率是多少？
令\(A_i\)表示第\(i\)个人取对自己的帽子， 要求概率的事件为 \[ B = \bigcap_{i=1}^n A_i^c . \] 由概率论的Jordan公式, \[\begin{aligned} P(B) =& 1 - P(B^c) = 1 - P(\bigcup_{i=1}^n A_i) \\ =& 1 - \left[\sum_{i} P(A_i) - \sum_{i<j} P(A_i A_j) + \sum_{i<j<k} P(A_i A_j A_k) + \cdots + (-1)^{n-1} P(A_1 A_2 \cdots A_n) \right]. \end{aligned}\]
根据对称性， \(P(A_i)=P(A_1)=\frac{1}{n}\), \(P(A_i A_j) = P(A_1 A_2) = \frac{1}{n} \frac{1}{n-1}\), \(P(A_i A_j A_k) = P(A_1 A_2 A_3) = \frac{1}{n} \frac{1}{n-1} \frac{1}{n-2}\), 有\(m\)个\(A_i\)的交集的概率为\(\frac{(n-m)!}{n!}\)， 这样的交集共有\(\binom{n}{m} = \frac{n!}{m!(n-m)!}\)个， 所以 \[\begin{aligned} P(B) =& 1 - \sum_{m=1}^n (-1)^{m-1} \frac{n!}{m!(n-m)!} \frac{(n-m)!}{n!} \\ =& 1 - \sum_{m=1}^n (-1)^{m-1} \frac{1}{m!} \\ =& \sum_{m=0}^n (-1)^m \frac{1}{m!} . \end{aligned}\]
用计算机程序验证此公式。
26.1.1.2 理论公式的程序
function hat_prob_th(n)
local s, a
s = 1//1
a = 1//1
sn = 1
for k=1:n
a *= -1 // k
s += a
end
return s
end
hat_prob_th.(1:10) |> show
## Rational{Int64}[0//1, 1//2, 1//3, 3//8, 11//30,
## 53//144, 103//280, 2119//5760, 16687//45360, 16481//44800]26.1.1.3 用穷举法计算理论值
function hat_exaust(n::Int)
local Ai_vec, Ω, Ai
# 生成样本空间为1:n的所有排列
Ω = collect(permutations(1:n))
# 生成各个事件为所包含的基本事件的集合
Ai_vec = []
for i=1:n
# 对$A_i$
Ai = []
for ω ∈ Ω
if ω[i] == i
push!(Ai, ω)
end
end
push!(Ai_vec, Ai)
end # for i
noMatch = setdiff(Ω, union(Ai_vec...))
length(noMatch) // length(Ω)
end
hat_exaust.(1:10) |> show
## Rational{Int64}[0//1, 1//2, 1//3, 3//8, 11//30,
## 53//144, 103//280, 2119//5760, 16687//45360, 16481//44800]上面的程序中，
用\(1:n\)的所有\(n!\)个排列组成了样本空间\(\Omega\)，
对每个事件\(A_i\)，
都对每一个\(\omega \in \Omega\)查看是否\(\omega \in \Omega\)，
如果\(\omega \in \Omega\)就把\(\omega\)存入集合\(A_i\)中。
这样，
把所有\(A_i\)保存为一个集合的序列Ai_vec，
然后用union(Ai_vec...)计算所有这些\(A_i\)的并集，
用setdiff(Ω, union(Ai_vec...))计算所有\(A_i^c\)的交集，
用古典概型估计\(P(B)\)。
26.1.1.4 随机模拟方法
# 某种取法是否有人取对的函数
function any_match(perm, n)
return any(perm .== collect(1:n))
end
# 对n个人的问题模拟计算概率的函数
function hat_sim(n, N=1_000_000)
1 - sum(any_match(randperm(n), n) for _ in 1:N) / N
end
# 测试，取n=1:10
nn = 10
p_th = round.(float.(hat_prob_th.(1:nn)), digits=4)
p_sim = round.(hat_sim.(1:nn), digits=4)
[1:nn p_th p_sim]结果：
10×3 Matrix{Float64}:
1.0 0.0 0.0
2.0 0.5 0.4992
3.0 0.3333 0.333
4.0 0.375 0.3745
5.0 0.3667 0.3665
6.0 0.3681 0.3681
7.0 0.3679 0.3681
8.0 0.3679 0.3688
9.0 0.3679 0.3671
10.0 0.3679 0.3675注意到理论公式的值当\(n\)较大时 接近于\(e^{-1}\)的泰勒展开级数， 所以\(n\)较大时所有人都没有取对的概率近似等于\(e{^-1} = 0.3678794\)。
26.1.2 可重复分组问题
26.1.2.1 问题和推导
设有\(n\)个编号的球， 有放回随机抽取\(m\)次， 记录第\(i\)个球出现的次数为\(X_i\)， 考虑随机向量\(\boldsymbol X = (X_1, X_2, \dots, X_n)\)的分布。
这个模型也可以表述为有\(m\)个不可区分的球， 有\(n\)个编号的箱子， 将\(m\)个球依次随机地放入箱子， 问每个箱子中球的个数组合的分布。
这是多项分布， 所以 \[\begin{aligned} & P(X_1=k_1, X_2=k_2, \dots, X_n=k_n) \\ =& \binom{m}{k_1, k_2 \dots, k_n} (\frac{1}{n})^{k_1} (\frac{1}{n})^{k_2} \dots (\frac{1}{n})^{k_n} \\ =& \frac{m!}{k_1!k_2! \dots k_n!} n^{-m} . \end{aligned}\]
验证这个概率公式。
26.1.2.2 穷举所有取值
用穷举所有\((k_1, k_2, \dots, k_n)\)组合计算所有对应概率， 然后检查和是否等于1来检查这个公式有没有错误。
如何穷举\((k_1, k_2, \dots, k_n)\)组合？ 其中\(k_1 + k_2 + \dots + k_n = m\)。 先看有多少个这样的组合。 我们想象有\(n\)个箱子排成一排并编号\(1:n\)， 每次抽取的球放入对应号码的箱子中（因为是有放回，所以可以认为是放了同号码球的复制品）。 用\(n-1\)个1代表两个箱子之间的边界， 用\(m\)个0代表箱子中的球， 用\(n+m-1\)个0-1序列代表某种结果。 比如\(010010\)表示\(n=3\), \(m=4\)， \(X_1=1\), \(X_2=2\), \(X_3=1\)。 显然，总的可能个数是从\(n+m-1\)个二进制位中选出\(n-1\)个1的位置的组合数\(\binom{n+m-1}{n-1}\)。
计算所有结果数的函数， 和某一个结果的概率的函数如下：
using Combinatorics: multinomial, combinations
function n_xvalues(n, m)
binomial(n+m-1, n-1)
end
function prob(x::Vector, n, m)
@assert sum(x)==m
multinomial(x...) / n^m
end下面的函数将0-1序列表示的结果转换为多项分布向量结果， 并罗列多项分布所有可能结果：
# 从0-1数列表示转换成(x1, ..., xn)的表示
function bin2mult(bits::Vector)
n = sum(bits) + 1
x = zeros(Int, n)
kx = 1
for bit in bits
if bit == 0
x[kx] += 1
else
kx += 1
end
end
return x
end
# 罗列所有可能组合：
function make_xvalues(n,m)
# 找出0-1序列中n-1个1的位置，穷举所有可能
npos = Combinatorics.combinations(1:(n+m-1), n-1)
xs = Vector{Int}[]
bits = zeros(Int, n+m-1)
for pos in npos
bits[:] .= 0
bits[pos] .= 1
push!(xs, bin2mult(bits))
end
return xs
end下面的函数对某对\(n, m\)值计算分布概率和：
function prob_sum(n, m)
xs = make_xvalues(n, m)
probs = map(x -> prob(x, n, m), xs)
sum(probs)
end下面的程序对多个\(n, m\)组合检查概率和是否等于1：
for n=2:10, m=1:20
if !(prob_sum(n,m) ≈ 1.0)
@show n, m, prob_sum(n,m);
end
end结果发现当\(n \geq 9\), \(m \geq 20\)时发生错误。 这应该是累积误差造成的， 对\(n=10, m=20\)， 共有10015005种不同结果。
26.1.2.3 模拟验证
对某个\(n, m\)组合， 模拟\(N\)次， 然后比较所有可能结果的理论概率与模拟样本估计的频率。
using DataFrames
function sim(n, m, N=1_000_000)
# 随机向量所有可能取值组合的列表
xvalues = make_xvalues(n,m)
# 用字典来保存从随机向量取值组合到概率或频数的
prob_th = Dict(Tuple(x) => prob(x, n, m) for x in xvalues)
xfreqs = Dict(Tuple(x) => 0 for x in xvalues)
for i=1:N
# xvs用来保存某一次模拟的随机向量取值
xvs = zeros(Int, n)
# 随机有放回地从1:n中抽取m次
draws = [rand(1:n) for _ in 1:m]
for k=1:n
xvs[k] += count(d -> d==k, draws)
end
xfreqs[Tuple(xvs)] += 1
end
# 下面将随机向量每个取值组合、理论概率、模拟估计频率列表
xv_mat = reduce(vcat, x' for x in xvalues)
pr0 = zeros(size(xv_mat, 1))
pr1 = zeros(size(xv_mat, 1))
for i=1:size(xv_mat,1)
xvt = Tuple(xv_mat[i,:])
pr0[i] = prob_th[xvt]
pr1[i] = xfreqs[xvt] / N
end
xv_df = DataFrame(reduce(vcat, x' for x in xvalues),
:auto)
xv_df[!,:prob_th] = pr0
xv_df[!,:prob_sim] = pr1
xv_df
end
# 测试
Random.seed!(101)
sim(3,4)结果：
15×5 DataFrame
Row │ x1 x2 x3 prob_th prob_sim
│ Int64 Int64 Int64 Float64 Float64
─────┼──────────────────────────────────────────
1 │ 0 0 4 0.0123457 0.012519
2 │ 0 1 3 0.0493827 0.049381
3 │ 0 2 2 0.0740741 0.074013
4 │ 0 3 1 0.0493827 0.049132
5 │ 0 4 0 0.0123457 0.012562
6 │ 1 0 3 0.0493827 0.049522
7 │ 1 1 2 0.148148 0.148013
8 │ 1 2 1 0.148148 0.148371
9 │ 1 3 0 0.0493827 0.049498
10 │ 2 0 2 0.0740741 0.074574
11 │ 2 1 1 0.148148 0.148458
12 │ 2 2 0 0.0740741 0.07349
13 │ 3 0 1 0.0493827 0.049227
14 │ 3 1 0 0.0493827 0.049107
15 │ 4 0 0 0.0123457 0.012133 26.1.3 多人博弈连续胜出问题
例子来自Ross Introduction to Probablity Models 12th Ed P.439 例7.7。
设有\(n\)个人比赛， 每一轮有且仅有一个人获胜， 每个人在一轮中获胜的概率为\(p_1, \dots, p_n\)， 每当同一个人\(k\)次连续获胜时就结束一局比赛， 开始新的一局，\(k \geq 2\)。 问每个人赢得一局的概率是多少， 平均多长时间完成一局比赛。
用全期望公式和更新过程理论可以推导出每个人赢得一局的概率为 \[ w_i \propto \frac{p_i^k (1 - p_i)}{1 - p_i^k} . \] 平均每局的时间为 \[ E(T) = \left( \sum_{i=1}^n \frac{p_i^k (1 - p_i)}{1 - p_i^k} \right)^{-1} . \]
用随机模拟方法进行验证。
程序和结果如下：
# 多人连胜游戏。每次某人连胜k轮时结束一局比赛。
# 输入：
# pw: 每人的获胜概率
# k: 连胜多少轮算一局结束
# N: 模拟重复的轮数
using Random, Distributions
using Statistics, StatsBase
function simone(pw::Vector; k=2, N=1000_000)
n = length(pw)
rates = pw .^ k .* (1 .- pw) ./ (1 .- pw .^ k)
# 每人理论获胜概率
probs = rates ./ sum(rates)
# 平均每局时间的理论值
mean_time = 1 / sum(rates)
# 一次性生成N轮的输赢
xv = StatsBase.sample(1:n, StatsBase.Weights(pw),
N; replace=true)
winners = Int[]
times = Int[]
# 下面循环找到每一局，条件是最后的k个为相同胜者
istart = 1
i = k
while i <= N
if all(xv[i-k+1:i-1] .== xv[i])
## 找到k个连续，结束一局，开始新的一局
push!(winners, xv[i])
push!(times, i - istart + 1)
istart = i+1
i += k
else
i += 1
end
end # while i
# 获胜者频率估计
freqs = counts(winners, 1:n)
probs_est = freqs ./ sum(freqs)
# 平均每局时间估计
mean_time_est = mean(times)
println("理论概率与模拟估计概率：")
display(round.([probs probs_est], digits=6))
println("理论平均时间与模拟估计值：")
show(round.([mean_time, mean_time_est], digits=4))
return
end
function test(k=2, N=1000)
pw = [1/2, 1/3, 1/6]
simone(pw, k=k, N=N)
end
Random.seed!(101)
test(2, 1000_000)
## 理论概率与模拟估计概率：
## 3×2 Matrix{Float64}:
## 0.608696 0.608658
## 0.304348 0.303514
## 0.086957 0.087829
## 理论平均时间与模拟估计值：
## [3.6522, 3.6565]
test(3, 1000_000)
## 理论概率与模拟估计概率：
## 3×2 Matrix{Float64}:
## 0.707595 0.706294
## 0.254008 0.254546
## 0.038397 0.03916
## 理论平均时间与模拟估计值：
## [9.9063, 9.9466]26.2 统计学例子
26.2.1 威布尔参数估计
考虑威布尔分布的矩估计问题。 威布尔分布密度函数为 \[\begin{aligned} f(x; m, \eta) = \frac{m}{\eta} \left( \frac{x}{\eta} \right)^{m-1} \exp\left[-\left( \frac{x}{\eta} \right)^m \right], \ x > 0. \end{aligned}\] 期望和方差为 \[\begin{aligned} EX =& \Gamma(1 + \frac{1}{m}), \\ \text{Var}(X) =& \eta^2 \left[ \Gamma(1 + \frac{2}{m}) - (\Gamma(1 + \frac{1}{m}))^2 \right] . \end{aligned}\]
26.2.1.1 矩估计法
设样本均值、方差分别为\(\bar X\), \(S^2\)， 求解关于\((m,\eta)\)的非线性方程组 \[\begin{aligned} \Gamma(1 + \frac{1}{m}) =& \bar X, \\ \eta^2 \left[ \Gamma(1 + \frac{2}{m}) - (\Gamma(1 + \frac{1}{m}))^2 \right] =& S^2 . \end{aligned}\] 可以求得参数的矩估计\((\hat m, \hat\eta)\)。
这有一个危险， 就是这个非线性方程组可能解不唯一， 求解的迭代算法可能不收敛到方程组的零点。 举例说明：
using Random, Distributions, SpecialFunctions, Statistics
Random.seed!(101)
dist = Weibull(2.0, 4.0)
samp = rand(dist, 100)
xbar = mean(samp)
S2 = var(samp)
xbar, S2
## (3.8028830034916767, 3.567696384660564)用NLsolve包求解：
using NLsolve
function eqnf!(F, x)
m, η = x
F[1] = gamma(1 + 1.0/m) - xbar
F[2] = η^2 * (gamma(1 + 2/m) - gamma(1 + 1/m)^2) - S2
return
end
x0 = [1.0, 1.0]
res = nlsolve(eqnf!, x0, autodiff=:forward)其中eqnf!函数将计算得到的两个方程返回到输入的自变量F中。
x0为初值。
结果：
Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [1.0, 1.0]
* Zero: [0.3816718560241766, 0.14606373284767674]
* Inf-norm of residuals: 0.000000
* Iterations: 15
* Convergence: true
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 12
* Jacobian Calls (df/dx): 8结果收敛了， 方程组的结果的函数值都等于0。 但是， 估计的参数值为\(m=0.38\), \(\eta=0.15\)， 与真实值\((2,4)\)差距很大， 说明方程组解不唯一。
换一组初值， 使用模拟时的参数真实值作为初值：
x0 = [2.0, 4.0]
res = nlsolve(eqnf!, x0, autodiff=:forward)Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [2.0, 4.0]
* Zero: [0.5756992788212827, 0.7052533506778957]
* Inf-norm of residuals: 2.211473
* Iterations: 8
* Convergence: false
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: false
* Function Calls (f): 9
* Jacobian Calls (df/dx): 4结果没有收敛。 这说明威布尔分布不适合使用矩估计法。
26.2.1.2 最大似然估计法
设样本为\(x_1, x_2, \dots, x_n\)。 对数似然函数为 \[ \log L(m, \eta) = n \log(m) - n m \log(\eta) + (m-1) \sum_i \log x_i - \eta^{-m} \sum x_i^m . \] 为避免\(m>0\), \(\eta>0\)的约束， 令\(a = \log m\), \(b=\log\eta\)， 取目标函数为关于\(a, b\)的负对数似然函数： \[ h(a, b) = -n a + n m b - (m - 1) \sum_i \log x_i + \eta^{-m} \sum_i x_i^m . \] 用Optim包求解上述目标函数的最小值点。
using Optim, ForwardDiff
function mle_wei(x, x0=log.([1.0, 1.0]))
n = length(x)
sumlogx = sum(log.(x))
# 负对数似然函数
function objf(xv)
a, b = xv
m = exp(a)
eta = exp(b)
-n*a + n*m*b - (m-1) * sumlogx + sum(x .^ m) / eta^m
end
ores = Optim.optimize(objf, [0.0, 0.0], Newton(); autodiff = :forward)
#ores = Optim.optimize(objf, [0.0, 0.0])
display(ores)
a, b = Optim.minimizer(ores)
[exp(a), exp(b)]
end
mle_wei(samp) |> show * Status: success
* Candidate solution
Final objective value: 1.996358e+02
* Found with
Algorithm: Newton's Method
* Convergence measures
|x - x'| = 3.69e-10 ≰ 0.0e+00
|x - x'|/|x'| = 2.53e-10 ≰ 0.0e+00
|f(x) - f(x')| = 2.84e-14 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 1.42e-16 ≰ 0.0e+00
|g(x)| = 1.42e-13 ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 6
f(x) calls: 20
∇f(x) calls: 20
∇²f(x) calls: 6
[2.138207889557369, 4.299549737307604]估计的参数值为\((2.14, 4.30)\)。
26.2.1.3 最小二乘法
威布尔分布函数为 \[ F(x) = 1 - \exp\left[ - \left( \frac{x}{\eta} \right)^m \right], \ x > 0, \] 生存函数为 \[ S(x) = \exp\left[ - \left( \frac{x}{\eta} \right)^m \right], \ x > 0 . \] 设样本的次序统计量为\(t_1 \leq t_2 \leq \cdots \leq t_n\)， 估计\(S(t_i)\)为 \[ \hat S(t_i) = \frac{n - i + \frac{5}{8}}{n + \frac{1}{4}}, \ n=1,2,\dots,n, \] 注意到 \[ \log [-\log S(x)] = m \log x - m \log\eta, \] 记 \(y_i = \log(-\log(\hat S(t_i)))\), \(x_i = \log t_i\)， 考虑回归问题 \[ y_i = a + b x_i, \] 其中\(a = -m \log\eta\), \(b = m\)。 求解出\(a, b\)可得 \(m=b\), \(\eta = e^{-a/b}\)。
using GLM
using DataFrames
n = length(samp)
y = log.(-log.((n .- (1:n) .+ 5/8)/(n+1/4)))
x = log.(sort(samp))
df = DataFrame(x=x, y=y)
lmr = lm(@formula(y ~ x), df)
show(lmr)
a, b = coef(lmr)
m, eta = b, exp(-a/b)StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}},
GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{
Float64, Matrix{Float64}}}}, Matrix{Float64}}
y ~ 1 + x
Coefficients:
─────────────────────────────────────────────────────────────────────────
Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
─────────────────────────────────────────────────────────────────────────
(Intercept) -3.16163 0.0430582 -73.43 <1e-86 -3.24708 -3.07618
x 2.16524 0.032482 66.66 <1e-82 2.10078 2.2297
─────────────────────────────────────────────────────────────────────────
(2.165238989151171, 4.30671460249195)估计\(m = 2.17\), \(\eta = 4.31\)， 与最大似然估计结果相近。 最小二乘方法计算仅需要求解线性方程组， 计算量较小。
这个计算中的最小二乘回归是一个一元回归问题， 如果我们仅需要回归系数估计值， 可以不需要调用GLM包， 而是自己写一个小的一元回归系数计算函数：
function simple_reg(x, y)
n = length(x)
mx = mean(x)
my = mean(y)
b = sum((x .- mx) .* (y .- my)) / sum((x .- mx) .^2)
a = my - b * mx
a, b
end测试：
a, b = simple_reg(x, y)
m, eta = b, exp(-a/b)
## (2.1652389891511716, 4.30671460249195)26.2.2 三角形分布参数矩估计
考虑如下的一般三角形分布： \[ f(x; a, b, c) = \begin{cases} 2\frac{x-a}{(b-a)(c-a)}, & \ x \in [a, c), \\ 2\frac{b-x}{(b-a)(b-c)}, & \ x \in [c, b] . \end{cases} \]
可以计算得到前三阶中心矩为： \[\begin{aligned} m_1 =& \frac{1}{3}(a + b + c), \\ m_2 =& \frac{1}{6}(a^2 + b^2 + c^2 + ab + bc + ac), \\ m_3 =& \frac{1}{10}(a^3 + b^3 + c^3 + a^2 b + a^2 c + b^2 a + b^2 c + c^2 a + c^2 b + abc). \end{aligned}\]
用矩估计法估计\(a, b, c\)：
using Random, Distributions, Statistics
using NLsolve, ForwardDiff
Random.seed!(101)
# 生成模拟样本
a, b, c = 3, 5, 4
dist = TriangularDist(a,b,c)
n = 100
samp = rand(dist, n)
# k阶中心矩估计函数：
m_k(k,data) = 1/n*sum(data.^k)
mHats = [m_k(k,samp) for k in 1:3]
# 要求解的方程，函数值返回到F中
function eqnf!(F, x)
a, b, c = x
F[1] = 1/3*( a + b + c) - mHats[1]
F[2] = 1/6*( a^2 + b^2 + c^2 + a*b + b*c + c*a ) - mHats[2]
F[3] = 1/10*( a^3 + b^3 + c^3
+ a^2*(b+c) + b^2*(a+c) + c^2*(a+b) + a*b*c) - mHats[3]
return
end
nlOutput = nlsolve(eqnf!, [ minimum(samp), maximum(samp), mean(samp) ],
autodiff=:forward)
display(nlOutput)
sol = sort(nlOutput.zero)
aHat, cHat, bHat = sol
show((aHat, bHat, cHat))Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [3.1321401602888455, 4.847665660768514, 4.041969871032595]
* Zero: [3.0109895917758136, 5.051519606315913, 4.063400415006064]
* Inf-norm of residuals: 0.000000
* Iterations: 4
* Convergence: true
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 5
* Jacobian Calls (df/dx): 5
(3.0109895917758136, 5.051519606315913, 4.063400415006064)模拟的真实参数为\((a, b, c)= (3, 5, 4)\)。 估计结果\((a, b, c) = (3.01, 5.05, 4.06)\)。 注意一到三阶矩的公式关于\(a, b, c\)是对称的， 所以求得解后按\(a < c < b\)重排后输出估计结果。
26.2.3 伽马分布参数估计
设总体\(X\)服从\(\Gamma(\alpha, \lambda)\)分布，密度为 \[\begin{aligned} p(x; \alpha, \lambda) = \begin{cases} \frac{\lambda^\alpha x^{\alpha-1}}{\Gamma(\alpha)} e^{-\lambda x}, & x>0, \alpha>0, \lambda>0, \\ 0, & \text{其它} \end{cases} \end{aligned}\]
设有\(X\)的简单随机样本\(X_1, X_2, \dots, X_n\)，则对数似然函数为 \[\begin{aligned} l(\alpha, \lambda) =& n \alpha \log \lambda - n \log \Gamma(\alpha) + \alpha \sum \log X_i - \lambda \sum X_i - \sum \log X_i, \end{aligned}\] 给定\(\alpha>0\)，先关于\(\lambda\)求最大值。易见 \[\begin{aligned} \frac{\partial l}{\partial \lambda} =& \frac{n\alpha}{\lambda} - \sum X_i, \end{aligned}\] 解得稳定点\(\hat\lambda = \hat\lambda(\alpha) = \alpha/\bar X\), 其中\(\bar X = \frac{1}{n} \sum X_i\)。 又\(\frac{\partial^2 l}{\partial \lambda^2} = -\frac{n\alpha}{\lambda^2} < 0\), \(\forall \lambda>0\), 所以\(\hat\lambda(\alpha)\)是给定\(\alpha\)情形下\(l(\alpha,\lambda)\)关于自变量\(\lambda\)的最大值点。 令\(\tilde l(\alpha) = l(\alpha, \hat\lambda(\alpha))\), 则 \[\begin{aligned} \tilde l(\alpha) =& n \alpha \log\alpha - n \log \Gamma(\alpha) + \alpha \left[ \sum \log \frac{X_i}{\bar X} - n \right] - \sum \log X_i, \\ \tilde l'(\alpha) =& n \log\alpha - n \psi(\alpha) + \sum\log\frac{X_i}{\bar X}, \end{aligned}\] 其中\(\psi(\alpha) = \frac{d}{d\alpha} \log\Gamma(\alpha)\)称为digamma函数。 可以证明\(\tilde l'(\alpha)\)严格单调减， 方程\(\tilde l'(\alpha) = 0\)存在唯一解\(\hat\alpha\)， 且\(\hat\alpha = \mathop{\text{argmax}}_{\alpha>0} \tilde l(\alpha)\)， 从而\((\hat\alpha, \hat\lambda(\hat\alpha))\)为总体参数的最大似然估计。 原来二维的优化问题简化为一个一元函数\(\tilde l(\alpha)\)的优化问题， 可以用二分法或牛顿法求解方程\(\tilde l'(\alpha) = 0\)。
Julia模拟例子：
using Random, Distributions, SpecialFunctions, Statistics
Random.seed!(101)
α,λ=(2.0, 3.0)
n = 100
dist = Gamma(α, 1/λ)
samp = rand(dist, n)
xbar = mean(samp)
using Roots
c1 = sum(log.(samp ./ xbar))
function g(alpha)
n * log(alpha) - n*digamma(alpha) + c1
end
alphahat = Roots.find_zero(g, 1.0)
lambdahat = alphahat / xbar
alphahat, lambdahat
## (1.9778484307569368, 2.9858680740265013)26.2.4 二项分布参数三次方无偏估计
设总体服从二项分布\(\text{B}(m, \theta)\)， \(m \geq 3\)已知， \(0 < \theta < 1\)未知。 对样本\(X_1, \dots, X_n\)(\(n \geq 2\))， 求\(\theta^3\)的无偏估计。
先得到二项分布的矩母函数 \[ g(t) = (\theta e^t + 1 - \theta)^m , \] 利用矩母函数求出一、二、三阶矩， 有 \[\begin{aligned} m \theta =& E(X), \\ m(m-1) \theta^2 =& E(X^2) - E(X) , \\ \theta^3 =& \frac{1}{m(m-1)(m-2)}[E(X^3) - 3 E(X^2) + 2 E(X)] . \end{aligned}\]
用程序验证：
using Random, Distributions
using Statistics, StatsBase
function pcest(x, m)
n = length(x)
V3 = sum(x .^ 3) / n
V1 = sum(x) / n
V2 = sum(x .^ 2) / n
#1/(m*(m-1)*(m-2))*(v3 - (3*m-2)*xbar + 3*(m-1)*S2)
1/(m*(m-1)*(m-2))*(V3 - 3V2 + 2V1)
end
function mean_pce(theta; m = 10, n=15, N=100_000)
dist = Binomial(m, theta)
mean([pcest(rand(dist, n), m) for i=1:N])
end
Random.seed!(101)
[[theta^3, mean_pce(theta)] for theta=0.1:0.1:0.9]9-element Vector{Vector{Float64}}:
[0.0010000000000000002, 0.0009987555555555552]
[0.008000000000000002, 0.00795673888888889]
[0.027, 0.027034577777777777]
[0.06400000000000002, 0.06390299444444444]
[0.125, 0.12496742222222224]
[0.216, 0.215703]
[0.3429999999999999, 0.34297825]
[0.5120000000000001, 0.5120347444444445]
[0.7290000000000001, 0.7291053055555556]26.2.5 Satterthwaite两样本t检验统计量分布研究
统计学中一个常见问题是独立的两个正态总体的均值是否相等的检验。 常用的方法是假设方差相等， 然后使用独立两样本t检验法。 但实际中两个总体方差不一定相等。 Satterthwaite检验法使用如下的近似t检验统计量， 不要求方差相等： \[ T = \frac{\bar x - \bar y}{\sqrt{ \frac{S_x^2}{n_1} + \frac{S_y^2}{n_2} }} . \] 在\(\mu_x = \mu_y\)的零假设成立时， \(T\)近似服从\(t(v)\)分布， 其中 \[ v = \frac{ \left( \frac{S_x^2}{n_1} + \frac{S_y^2}{n_2}\right)^2 }{ \frac{S_x^2/n_1}{n_1 - 1} + \frac{S_y^2/n_2}{n_2 - 1} } . \] \(v\)不一定是整数值。 \(t\)分布可以推广到允许自由度非整数值。
这个近似分布的精确程度如何？ 我们针对不同的样本量， 不同的方差比值做一些模拟研究。
我们对给定的\(n_1, n_2\), \(\sigma_x^2\), \(\sigma_y\)平方， 模拟产生\(N\)组样本， 每组样本产生相应的T统计量\(t_i\)和相应的自由度\(v_i\)， \(i=1,2,\ldots,N\)。 用模拟方法研究\(T\)统计量是否服从\(t(v)\)分布。 这里的困难是， 每个不同的\(t_i\)服从不同的t分布， 所以无法直接做QQ图， 或者作直方图然后叠加理论分布密度图。 为此， 我们计算\(t_i\)的正态分数： 因为理论上\(t_i\)服从\(t(v_i)\)分布， 所以\(F_t(t_i, v_i)\)应服从U(0,1)均匀分布， 其中\(F_t(\cdot, n)\)是自由度为\(n\)的t分布的分布函数。 于是，\(\Phi^{-1}(F_t(t_i, v_i))\)应该服从标准正态分布， 其中\(\Phi(\cdot)\)是标准正态分布函数， \(\Phi^{-1}\)是标准正态分布分位数函数。 这样， 令\(z_i = \Phi^{-1}(F_t(t_i, v_i))\)， \(\{ z_i, i=1,2,\ldots, N \}\)应表现为来自标准正态分布的随机样本。 我们用正态QQ图验证。
using Statistics
# 计算Satterthwaite两样本t统计量和近似自由度
function tstat_satt(x, y)
n1 = length(x)
n2 = length(y)
sx2 = var(x)
sy2 = var(y)
v = (sx2/n1 + sy2/n2)^2 /
((sx2/n1)^2/(n1-1) + (sy2/n2)^2/(n2-1))
t = (mean(x) - mean(y)) / sqrt(sx2/n1 + sy2/n2)
return (t, v)
end
# 将t统计量按理论分布转换为正态分布
using Distributions
z_satt(t, v) = quantile(Normal(), cdf(TDist(v), t))
# 模拟
function sim_satt(n1=10, n2=10, sx2=1.0, sy2=1.0; N=10_000)
t_vec = zeros(N)
v_vec = zeros(N)
q_vec = zeros(N)
for i=1:N
x = randn(n1) .* sqrt(sx2)
y = randn(n2) .* sqrt(sy2)
t, v = tstat_satt(x, y)
t_vec[i] = t
v_vec[i] = v
end
z_vec = [z_satt(t_vec[i], v_vec[i]) for i=1:N]
# 作Z分数的正态QQ图
fig, ax, plt = qqnorm(z_vec, qqline=:identity)
display(fig)
fig
end相等样本量， 相等方差的情形：
using Random
Random.seed!(101)
using CairoMakie
CairoMakie.activate!()
sim_satt(10, 10, 1.0, 1.0);
相等样本量，2倍方差的情形：
sim_satt(10, 10, 2.0, 1.0);
不等样本量，2倍方差的情形：
sim_satt(20, 10, 2.0, 1.0);