使用Optim包实现威布尔分布极大似然估计的代码报错排查

嘿，我来帮你搞定这个问题！你遇到的DomainError其实是两个小问题叠加导致的，咱们一步步拆解清楚：

错误根源 #

1. 参数约束缺失：威布尔分布的尺度参数λ和形状参数k必须是严格正数，但你用的默认Nelder-Mead优化算法是无约束的——迭代过程中它会尝试各种参数值，包括非正的λ或k。这时候(x/lambda)会变成负数（毕竟x是从威布尔分布生成的正数），当k-1不是整数时，负数的非整数次幂会产生复数，直接触发DomainError。
2. PDF公式写错了：你定义的pdf_weibull里，指数部分是exp((-x/lambda)^k)，正确的威布尔PDF应该是exp(-(x/lambda)^k)——少了个括号导致符号完全错了，这会让似然函数计算彻底偏离正确结果！

解决方案：两种靠谱的修复方式 #

方案1：用带边界约束的优化算法（推荐） #

Optim里的L-BFGS-B算法支持设置参数的上下界，我们直接把λ和k的下界设为极小的正数（比如1e-6），避免算法探索非法参数值：

using Optim, Distributions

# 生成模拟数据
rng = MersenneTwister(1234);
d = Weibull(1,1)
x = rand(d,1000)

# 修正后的威布尔PDF函数
function pdf_weibull(x_val, lambda, k)
    k/lambda * (x_val/lambda)^(k-1) * exp(-(x_val/lambda)^k)
end

# 负对数似然函数
function neg_log_likelihood(x::Vector, lambda, k)
    total = 0.0
    for val in x
        total += log(pdf_weibull(val, lambda, k))
    end
    -total  # 极大似然等价于最小化负对数似然
end

# 包装成Optim需要的单参数函数
obj_func(vars) = neg_log_likelihood(x, vars[1], vars[2])

# 使用Fminbox+LBFGS设置参数边界：λ和k都≥1e-6
result = optimize(obj_func, [1e-6, 1e-6], [Inf, Inf], [1.0, 1.0], Fminbox(LBFGS()))

# 查看估计结果
println("估计的尺度参数λ: ", result.minimizer[1])
println("估计的形状参数k: ", result.minimizer[2])

方案2：对数变换参数（无约束优化） #

如果不想用约束算法，我们可以对λ和k取对数，把正参数空间映射到整个实数空间——这样不管算法怎么迭代，转换后的参数都是合法的正数：

using Optim, Distributions

# 生成模拟数据
rng = MersenneTwister(1234);
d = Weibull(1,1)
x = rand(d,1000)

# 修正后的威布尔PDF函数
function pdf_weibull(x_val, lambda, k)
    k/lambda * (x_val/lambda)^(k-1) * exp(-(x_val/lambda)^k)
end

# 基于对数参数的负对数似然
function neg_log_likelihood_log(x::Vector, log_lambda, log_k)
    lambda = exp(log_lambda)
    k = exp(log_k)
    total = 0.0
    for val in x
        total += log(pdf_weibull(val, lambda, k))
    end
    -total
end

# 包装函数
obj_func_log(vars) = neg_log_likelihood_log(x, vars[1], vars[2])

# 初始值对应λ=1、k=1的对数（即0.0）
result = optimize(obj_func_log, [0.0, 0.0])

# 转换回原始参数
lambda_est = exp(result.minimizer[1])
k_est = exp(result.minimizer[2])

println("估计的尺度参数λ: ", lambda_est)
println("估计的形状参数k: ", k_est)

验证结果 #

两种方法跑出来的估计值应该很接近真实值λ=1、k=1（毕竟模拟数据是用这两个参数生成的），你可以自己测试一下~

内容的提问来源于stack exchange，提问作者Pedro Rafael