2

我是斯坦的新手，所以希望你能指出我正确的方向。我会根据我的情况来确保我们在同一页面上......

如果我有一组单变量法线，文档会告诉我：

y ~ normal(mu_vec, sigma);

提供与未矢量化版本相同的模型：

for (n in 1:N)
    y[n] ~ normal(mu_vec[n], sigma);

但是矢量化版本（很多？）更快。好的，很好，很有道理。

所以第一个问题是：是否有可能在单变量正常情况下利用这种向量化加速，其中样本的mu和都随向量中的位置而变化。sigma即如果mu_vec和sigma_vec都是向量（在前一种情况下sigma是标量），那么是这样的：

y ~ normal(mu_vec, sigma_vec);

相当于：

for (n in 1:N)
    y[n] ~ normal(mu_vec[n], sigma_vec[n]);

如果是这样，是否有可比的加速？

好的。那就是热身。真正的问题是如何最好地接近上述的多变量等价物。

在我的特殊情况下，我N观​​察到一些变量的双变量数据y，我将其存储在N x 2矩阵中。（对于数量级，N大约1000在我的用例中。）

我的信念是，每个观察的每个组成部分的平均值是，每个组成部分的标准0差是每个观察结果1（我很高兴将它们硬编码为这样）。然而，我的信念是，相关性（ ）作为另一个观察变量（存储在元素向量中）rho的（简单）函数，因观察而异。例如，我们可能会说，我们的目标是从我们的数据中学习。即第 th 观察的协方差矩阵将是：xNrho[n] = 2*inverse_logit(beta * x[n]) - 1n in 1:Nbetan

[1,      rho[n]]
[rho[n], 1     ]

我的问题是将它放在一个 STAN 模型中的最佳方法是什么，这样它就不会那么慢了？是否有multi_normal分布的矢量化版本，以便我可以将其指定为：

y ~ multi_normal(vector_of_mu_2_tuples, vector_of_sigma_matrices)

或者也许是其他类似的表述？或者我需要写：

for (n in 1:N)
    y[n] ~ multi_normal(vector_of_mu_2_tuples[n], vector_of_sigma_matrices[n])

在设置vector_of_sigma_matrices并vector_of_mu_2_tuples在较早的块之后？

提前感谢您的任何指导！

编辑以添加代码

使用python，我可以按照我的问题的精神生成数据，如下所示：

import numpy as np
import pandas as pd
import pystan as pys
import scipy as sp
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import seaborn as sns

def gen_normal_data(N, true_beta, true_mu, true_stdevs):
    N = N

    true_beta = true_beta
    true_mu = true_mu
    true_stdevs = true_stdevs

    drivers = np.random.randn(N)
    correls = 2.0 * sp.special.expit(drivers*true_beta)-1.0

    observations = []
    for i in range(N):
        covar = np.array([[true_stdevs[0]**2, true_stdevs[0] * true_stdevs[1] * correls[i]],
                          [true_stdevs[0] * true_stdevs[1] * correls[i], true_stdevs[1]**2]])
        observations.append(sp.stats.multivariate_normal.rvs(true_mu, covar, size=1).tolist())
    observations = np.array(observations)
    
    return {
        'N': N,
        'true_mu': true_mu,
        'true_stdev': true_stdevs,
        'y': observations,
        'd': drivers,
        'correls': correls
    }

然后使用以下方法实际生成数据：

normal_data = gen_normal_data(100, 1.5, np.array([1., 5.]), np.array([2., 5.]))

这是数据集的样子（左窗格中的y彩色散点图和右窗格中的 by...所以想法是越高，越接近 ，越低，越接近。所以期望左窗格上的红点为“从左下到右上”，蓝点为“从左上到右下”，实际上它们是：correlsdriversdriver1correldriver-1correl

fig, axes = plt.subplots(1, 2, figsize=(15, 6))
x = normal_data['y'][:, 0]
y = normal_data['y'][:, 1]
correls = normal_data['correls']
drivers = normal_data['d']

for ax, colordata, cmap in zip(axes, [correls, drivers], ['coolwarm', 'viridis']):
    color_extreme = max(abs(colordata.max()), abs(colordata.min()))
    sc = ax.scatter(x, y, c=colordata, lw=0, cmap=cmap, vmin=-color_extreme, vmax=color_extreme)
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='5%', pad=0.05)
    fig.colorbar(sc, cax=cax, orientation='vertical')
fig.tight_layout()

散点图

使用蛮力方法，我可以建立一个如下所示的 STAN 模型：

model_naked = pys.StanModel(
    model_name='naked',
    model_code="""
data {
    int<lower=0> N;
    vector[2] true_mu;
    vector[2] true_stdev;
    real d[N];
    vector[2] y[N];
}

parameters {
    real beta;
}

transformed parameters {
}

model {
    real rho[N];
    matrix[2, 2] cov[N];
    
    for (n in 1:N) {
        rho[n] = 2.0*inv_logit(beta * d[n]) - 1.0;
 
        cov[n, 1, 1] = true_stdev[1]^2;
        cov[n, 1, 2] = true_stdev[1] * true_stdev[2] * rho[n];
        cov[n, 2, 1] = true_stdev[1] * true_stdev[2] * rho[n];
        cov[n, 2, 2] = true_stdev[2]^2;
    }
    
    beta ~ normal(0, 10000);
    for (n in 1:N) {
        y[n] ~ multi_normal(true_mu, cov[n]);
    }
}
"""
)

这很适合：

fit_naked = model_naked.sampling(data=normal_data, iter=1000, chains=2)f = fit_naked.plot();
f.tight_layout()

跟踪测试版

但我希望有人能为“边缘化”方法指出正确的方向，在这种方法中，我们将双变量法线分解为一对可以使用相关性混合的独立法线。我需要这个的原因是在我的实际用例中，两个维度都是肥尾的。我很高兴将其建模为 student-t 分布，但问题是 STAN 只允许nu指定一个（不是每个维度一个），所以我认为我需要找到一种方法将 a 分解multi_student_t为一对独立student_t的，这样我就可以为每个维度分别设置自由度。

4

2 回答 2

2

单变量正态分布确实接受其任何或所有参数的向量，并且它比循环N观察以N使用标量参数调用它要快。

然而，加速只会是线性的，因为计算都是一样的，但如果你只调用一次，它只需要分配一次内存。总的墙时间受您必须执行的函数评估次数的影响更大，这取决于2^10 - 1每次 MCMC 迭代（默认情况下），但是您是否达到最大树深度取决于您尝试采样的后验分布的几何形状from，而这又取决于一切，包括您所依赖的数据。

双变量正态分布可以写成第一个变量的边际单变量正态分布和给定第一个变量的第二个变量的条件单变量正态分布的乘积。在 Stan 代码中，我们可以利用逐元素乘法和除法来编写它的对数密度，例如

target += normal_lpdf(first_variable | first_means, first_sigmas);
target += normal_lpdf(second_variable | second_means + 
  rhos .* first_sigmas ./ second_sigmas .* (first_variable - first_means),
  second_sigmas .* sqrt(1 - square(rhos)));

不幸的是，Stan 中更一般的多元正态分布没有输入协方差矩阵数组的实现。

于 2018-05-04T13:19:48.837 回答
2

这并不能完全回答您的问题，但是您可以通过删除一堆冗余计算并将比例转换为使用 tanh 而不是比例逆 logit 来提高程序的效率。我会摆脱缩放，只使用较小的 beta，但我留下了它，以便它应该得到相同的结果。

data {
  int<lower=0> N;
  vector[2] mu;
  vector[2] sigma;
  vector[N] d;
  vector[2] y[N];
}
parameters {
  real beta;
}
transformed data {
  real var1 = square(sigma[1]);
  real var2 = square(sigma[2]);
  real covar12 = sigma[1] * sigma[2];
  vector[N] d_div_2 = d * 0.5;
}
model {
  // note: tanh(u) = 2 * inv_logit(u / 2) - 1
  vector[N] rho = tanh(beta * d_div_2);
  matrix[2, 2] Sigma;
  Sigma[1, 1] = var1;
  Sigma[2, 2] = var2;
  // only reassign what's necessary with minimal recomputation
  for (n in 1:N) {
    Sigma[1, 2] = rho[n] * covar12;
    Sigma[2, 1] = Sigma[1, 2];
    y[n] ~ multi_normal(true_mu, Sigma);
  }
  // weakly informative priors fit more easily
  beta ~ normal(0, 8);
}

您还可以通过计算 Cholesky 分解作为函数rho和其他固定值并使用它来分解 - 它在多元正态中节省了求解器步骤。

另一个选择是直接写出 multi-student-t，而不是使用我们的内置实现。内置的可能不会快很多，因为整个操作很大程度上由矩阵求解主导。

于 2018-05-06T18:07:12.817 回答