为何Turing.jl中的分层贝叶斯Logit模型运行极慢？

问题背景 #

尝试构建贝叶斯逻辑回归模型分析个人支付次数与违约概率的关联，通过生成合成数据验证模型适配性。涉及变量：

• 二元变量payers：违约标记为1，否则为0
• 支付次数P
• 平均支付间隔avg_delay

合成数据生成代码 #

using Turing

is_bad_pay(x) = x > 70 ? 1 : 0

function simulate_payment_frequency(Pₙ, N)
    P = rand(DiscreteUniform(1, Pₙ), N)
    avg_delay = rand(LogNormal(2,2), N)
    payers = [is_bad_pay(x) for x in avg_delay]
    payers, P, avg_delay
end

模型定义代码 #

# fit simulated data
@model function freq_pay(prob_pay, number_payments, avg_delay)
    Num_payments = length(unique(number_payments))
    # hierarchical by quantile of number of payments
    αₛ ~ filldist(Normal(60, 10), Num_payments)
    βₛ ~ filldist(Normal(0, 1), Num_payments)
    v = @. logistic(αₛ[number_payments] + βₛ[number_payments]*(avg_delay))  
    # logistic regression 
    for i ∈ eachindex(v)
        prob_pay[i] ~ Bernoulli(v[i])
    end
end

问题现象 #

• 模拟仅含2次支付的样本时，模型拟合运行正常：

  synthetic_payers = simulate_payment_frequency(2, 100)
  s1_1 = sample(freq_pay(synthetic_payers[1], synthetic_payers[2], synthetic_payers[3]), NUTS(), 100)
  

• 模拟支付次数超过3次（如4次）的样本时，模型运行陷入停滞：

  synthetic_payers = simulate_payment_frequency(4, 100)
  s1_2 = sample(freq_pay(synthetic_payers[1], synthetic_payers[2], synthetic_payers[3]), NUTS(), 100)
  

核心问题与修正方案 #

1. 合成数据逻辑矛盾 #

你的合成数据中，payers仅由avg_delay决定，完全不依赖支付次数P。当支付次数类别增加，模型试图拟合P对违约概率的影响，但真实数据里无此关联，导致后验分布极度模糊，NUTS采样器无法找到有效路径，最终停滞。

2. 参数尺度不匹配 #

αₛ的先验为Normal(60,10)，而avg_delay由LogNormal(2,2)生成（均值约54.6），两者相加后数值过大，代入logistic函数会趋近于0或1，导致伯努利分布概率极端，梯度消失，采样效率骤降。

3. 分层样本量不足 #

当Pₙ=4、N=100时，每个支付次数类别的样本量仅约25个，分层模型样本量不足进一步加剧采样不稳定。

具体修正步骤 #

方案1：让合成数据包含支付次数的影响

修改数据生成逻辑，让payers与P关联，比如支付次数越多，违约概率越低：

function simulate_payment_frequency(Pₙ, N)
    P = rand(DiscreteUniform(1, Pₙ), N)
    # 支付次数越多，平均延迟的均值越低
    avg_delay = [rand(LogNormal(2 - 0.3*(p-1), 2)) for p in P]
    payers = [x > 70 ? 1 : 0 for x in avg_delay]
    payers, P, avg_delay
end

方案2：调整参数先验与变量尺度

缩小αₛ先验尺度，并对avg_delay标准化，避免数值极端：

@model function freq_pay(prob_pay, number_payments, avg_delay)
    Num_payments = length(unique(number_payments))
    αₛ ~ filldist(Normal(0, 5), Num_payments)
    βₛ ~ filldist(Normal(0, 1), Num_payments)
    # 标准化平均支付间隔，优化尺度匹配
    std_delay = (avg_delay .- mean(avg_delay)) ./ std(avg_delay)
    v = @. logistic(αₛ[number_payments] + βₛ[number_payments] * std_delay)  
    for i ∈ eachindex(v)
        prob_pay[i] ~ Bernoulli(v[i])
    end
end

方案3：增加样本量

提升模拟样本量，保证每个分层类别有足够数据：

synthetic_payers = simulate_payment_frequency(4, 500)
s1_2 = sample(freq_pay(synthetic_payers[1], synthetic_payers[2], synthetic_payers[3]), NUTS(), 200)

内容的提问来源于stack exchange，提问作者Emiliano Isaza Villamizar