机器算法验证 - 在 MCMC 中管理高自相关 - 吾爱随笔录

在 MCMC 中管理高自相关

机器算法验证自相关多层次分析马尔可夫链蒙特卡罗锯齿

2022-03-08 08:38:42

我正在为使用 R 和 JAGS 的元分析构建一个相当复杂的分层贝叶斯模型。简化一点，模型的两个关键级别有其中是第个观察值的终点（在这种情况下，转基因与非转基因作物产量），是研究 j 的影响 γ是各种研究级变量的影响（该国的经济发展状况）研究完成，作物种类，研究方法等）由函数族和

y_{i j} = α_{j} + ϵ_{i}

$y_{ij} = \alpha_j + \epsilon_i$

α_{j} = \sum_{h} γ_{h (j)} + ϵ_{j}

$\alpha_j = \sum_h \gamma_{h(j)} + \epsilon_j$

y_{i j}

$y_{ij}$

i

$i$

j

$j$

α_{j}

$\alpha_j$

j

$j$

γ

$\gamma$

h

$h$

ϵ

$\epsilon$ s 是误差项。请注意， s 不是虚拟变量的系数。相反，对于不同的研究水平值，存在不同的变量。例如，发展中国家有，发达国家\gamma_{developing} 。

γ

$\gamma$

γ

$\gamma$

γ_{d e v e l o p i n g}

$\gamma_{developing}$

γ_{d e v e l o p e d}

$\gamma_{developed}$

我主要对估计的值感兴趣。这意味着从模型中删除研究级别的变量不是一个好的选择。 $\gamma$

几个研究级变量之间存在高度相关性，我认为这在我的 MCMC 链中产生了很大的自相关性。此诊断图说明了链轨迹（左）和产生的自相关（右）：

由于自相关，我从 4 个链中的每个 10,000 个样本中获得了 60-120 的有效样本大小。

我有两个问题，一个明显客观，另一个更主观。

除了细化、添加更多链和运行采样器更长时间之外，我可以使用哪些技术来管理这个自相关问题？我所说的“管理”是指“在合理的时间内产生相当好的估计”。在计算能力方面，我在 MacBook Pro 上运行这些模型。
这种自相关程度有多严重？此处和John Kruschke 的博客上的讨论表明，如果我们运行模型足够长的时间，“块状自相关可能已经全部被平均化了”（Kruschke），所以这并不是什么大问题。

这是生成上述图的模型的 JAGS 代码，以防万一有人有兴趣深入了解细节：

model {
for (i in 1:n) {
    # Study finding = study effect + noise
    # tau = precision (1/variance)
    # nu = normality parameter (higher = more Gaussian)
    y[i] ~ dt(alpha[study[i]], tau[study[i]], nu)
}

nu <- nu_minus_one + 1
nu_minus_one ~ dexp(1/lambda)
lambda <- 30

# Hyperparameters above study effect
for (j in 1:n_study) {
    # Study effect = country-type effect + noise
    alpha_hat[j] <- gamma_countr[countr[j]] + 
                    gamma_studytype[studytype[j]] +
                    gamma_jour[jourtype[j]] +
                    gamma_industry[industrytype[j]]
    alpha[j] ~ dnorm(alpha_hat[j], tau_alpha)
    # Study-level variance
    tau[j] <- 1/sigmasq[j]
    sigmasq[j] ~ dunif(sigmasq_hat[j], sigmasq_hat[j] + pow(sigma_bound, 2))
    sigmasq_hat[j] <- eta_countr[countr[j]] + 
                        eta_studytype[studytype[j]] + 
                        eta_jour[jourtype[j]] +
                        eta_industry[industrytype[j]]
    sigma_hat[j] <- sqrt(sigmasq_hat[j])
}
tau_alpha <- 1/pow(sigma_alpha, 2)
sigma_alpha ~ dunif(0, sigma_alpha_bound)

# Priors for country-type effects
# Developing = 1, developed = 2
for (k in 1:2) {
    gamma_countr[k] ~ dnorm(gamma_prior_exp, tau_countr[k])
    tau_countr[k] <- 1/pow(sigma_countr[k], 2)
    sigma_countr[k] ~ dunif(0, gamma_sigma_bound)
    eta_countr[k] ~ dunif(0, eta_bound)
}

# Priors for study-type effects
# Farmer survey = 1, field trial = 2
for (k in 1:2) {
    gamma_studytype[k] ~ dnorm(gamma_prior_exp, tau_studytype[k])
    tau_studytype[k] <- 1/pow(sigma_studytype[k], 2)
    sigma_studytype[k] ~ dunif(0, gamma_sigma_bound)
    eta_studytype[k] ~ dunif(0, eta_bound)
}

# Priors for journal effects
# Note journal published = 1, journal published = 2
for (k in 1:2) {
    gamma_jour[k] ~ dnorm(gamma_prior_exp, tau_jourtype[k])
    tau_jourtype[k] <- 1/pow(sigma_jourtype[k], 2)
    sigma_jourtype[k] ~ dunif(0, gamma_sigma_bound)
    eta_jour[k] ~ dunif(0, eta_bound)
}

# Priors for industry funding effects
for (k in 1:2) {
    gamma_industry[k] ~ dnorm(gamma_prior_exp, tau_industrytype[k])
    tau_industrytype[k] <- 1/pow(sigma_industrytype[k], 2)
    sigma_industrytype[k] ~ dunif(0, gamma_sigma_bound)
    eta_industry[k] ~ dunif(0, eta_bound)
}
}

1个回答

根据 user777 的建议，我的第一个问题的答案似乎是“使用 Stan”。在 Stan 中重写模型后，以下是轨迹（老化后 4 链 x 5000 次迭代）：
以及自相关图：

好多了！为了完整起见，这里是 Stan 代码：

data {                          // Data: Exogenously given information
// Data on totals
int n;                      // Number of distinct finding i
int n_study;                // Number of distinct studies j

// Finding-level data
vector[n] y;                // Endpoint for finding i
int study_n[n_study];       // # findings for study j

// Study-level data
int countr[n_study];        // Country type for study j
int studytype[n_study];     // Study type for study j
int jourtype[n_study];      // Was study j published in a journal?
int industrytype[n_study];  // Was study j funded by industry?

// Top-level constants set in R call
real sigma_alpha_bound;     // Upper bound for noise in alphas
real gamma_prior_exp;       // Prior expected value of gamma
real gamma_sigma_bound;     // Upper bound for noise in gammas
real eta_bound;             // Upper bound for etas
}

transformed data {
// Constants set here
int countr_levels;          // # levels for countr
int study_levels;           // # levels for studytype
int jour_levels;            // # levels for jourtype
int industry_levels;        // # levels for industrytype
countr_levels <- 2;
study_levels <- 2;
jour_levels <- 2;
industry_levels <- 2;
}

parameters {                    // Parameters:  Unobserved variables to be estimated
vector[n_study] alpha;      // Study-level mean
real<lower = 0, upper = sigma_alpha_bound> sigma_alpha;     // Noise in alphas

vector<lower = 0, upper = 100>[n_study] sigma;          // Study-level standard deviation

// Gammas:  contextual effects on study-level means
// Country-type effect and noise in its estimate
vector[countr_levels] gamma_countr;     
vector<lower = 0, upper = gamma_sigma_bound>[countr_levels] sigma_countr;
// Study-type effect and noise in its estimate
vector[study_levels] gamma_study;
vector<lower = 0, upper = gamma_sigma_bound>[study_levels] sigma_study;
vector[jour_levels] gamma_jour;
vector<lower = 0, upper = gamma_sigma_bound>[jour_levels] sigma_jour;
vector[industry_levels] gamma_industry;
vector<lower = 0, upper = gamma_sigma_bound>[industry_levels] sigma_industry;


// Etas:  contextual effects on study-level standard deviation
vector<lower = 0, upper = eta_bound>[countr_levels] eta_countr;
vector<lower = 0, upper = eta_bound>[study_levels] eta_study;
vector<lower = 0, upper = eta_bound>[jour_levels] eta_jour;
vector<lower = 0, upper = eta_bound>[industry_levels] eta_industry;
}

transformed parameters {
vector[n_study] alpha_hat;                  // Fitted alpha, based only on gammas
vector<lower = 0>[n_study] sigma_hat;       // Fitted sd, based only on sigmasq_hat

for (j in 1:n_study) {
    alpha_hat[j] <- gamma_countr[countr[j]] + gamma_study[studytype[j]] + 
                    gamma_jour[jourtype[j]] + gamma_industry[industrytype[j]];
    sigma_hat[j] <- sqrt(eta_countr[countr[j]]^2 + eta_study[studytype[j]]^2 +
                        eta_jour[jourtype[j]] + eta_industry[industrytype[j]]);
}
}

model {
// Technique for working w/ ragged data from Stan manual, page 135
int pos;
pos <- 1;
for (j in 1:n_study) {
    segment(y, pos, study_n[j]) ~ normal(alpha[j], sigma[j]);
    pos <- pos + study_n[j];
}

// Study-level mean = fitted alpha + Gaussian noise
alpha ~ normal(alpha_hat, sigma_alpha);

// Study-level variance = gamma distribution w/ mean sigma_hat
sigma ~ gamma(.1 * sigma_hat, .1);

// Priors for gammas
gamma_countr ~ normal(gamma_prior_exp, sigma_countr);
gamma_study ~ normal(gamma_prior_exp, sigma_study);
gamma_jour ~ normal(gamma_prior_exp, sigma_study);
gamma_industry ~ normal(gamma_prior_exp, sigma_study);
}

其它你可能感兴趣的问题

上一篇如何验证神经网络的实现是否正常工作？下一篇理解分数差分公式