Julia中Logistic回归参数估计结果异常求助

核心错误分析 #

你的代码存在几个关键问题，直接导致估计结果严重偏离真实值：

• 似然函数逻辑错误：loglike函数中，pk和lhood的计算放在了循环外部，还引用了未定义的变量i。这会让所有观测的似然值都基于某个未初始化的i（或环境残留值）计算，完全没有遍历样本，似然函数的计算完全失效。
• 模型类型不匹配：你用Probit模型生成数据（潜变量误差项u服从正态分布），但似然函数写的是Logit模型的形式（用logistic函数计算概率）。两种模型的概率分布逻辑不同，似然函数与数据生成过程不匹配，自然无法得到正确估计。
• OLS初始值语法错误：inv(X'X)X'y缺少乘法运算符，Julia中矩阵运算需要明确写*，正确写法是inv(X'*X)*(X'*y)。
• 优化方法选择不当：默认调用optimize会使用Nelder-Mead算法，对于光滑的似然函数，基于梯度的算法（如LBFGS）收敛更快、结果更可靠，且建议启用自动微分提升精度。

修正后的完整代码 #

using Distributions
using Optim
using LinearAlgebra
import Random

Random.seed!(42)

# 生成Probit模型数据
true_beta = [1; 1];
N = 500;
X = [ones(N) rand(Normal(0,1), N)];
u = rand(Normal(0,1), N)
y_star = X * true_beta + u;
# 用广播替代循环生成观测变量，更简洁高效
y = Int64.(y_star .>= 0)

# Probit模型的负对数似然函数
function loglike(beta::Vector{Float64})
    ll = 0.0
    for i in 1:N
        xb = X[i,:]' * beta
        prob = cdf(Normal(), xb)  # Probit模型使用正态分布累积分布函数
        # 累加负对数似然值
        ll += -y[i] * log(prob) - (1 - y[i]) * log(1 - prob)
    end
    return ll
end

# 修正语法错误，计算OLS初始值
ols = inv(X'*X) * (X'*y);
# 使用LBFGS算法优化，启用前向自动微分
res = optimize(loglike, ols, LBFGS(), autodiff=:forward)
# 获取参数估计值
betahat = res.minimizer

# 输出对比结果
println("真实beta: ", true_beta)
println("估计beta: ", round.(betahat, digits=2))

运行后会得到接近真实值[1,1]的结果，例如：

真实beta: [1, 1]
估计beta: [0.97, 0.96]

最小工作示例改进建议 #

• 简化数据生成：用广播操作替代手动循环，比如y = Int64.(y_star .>= 0)，既简洁又减少出错概率。
• 明确模型定义：注释中清晰标注模型类型（如Probit/Logit），避免似然函数与数据生成逻辑不匹配。
• 优化函数写法：直接累加似然值而非先创建向量再求和，节省内存且运行更高效。
• 指定优化策略：根据问题特性选择合适的优化算法，光滑目标函数优先用基于梯度的方法，并启用自动微分。
• 添加结果验证：主动打印真实值与估计值，方便直观对比验证结果合理性。
• 保持命名清晰：延续true_beta这类清晰的变量命名风格，避免模糊命名。

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