Я новичок в Стене, так что надеюсь, что вы можете указать мне правильное направление. Я создам ситуацию, чтобы убедиться, что мы на одной странице ...
Если бы у меня была коллекция одномерных нормалей, документы говорят мне, что:
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]);
и если да, то есть ли сопоставимое ускорение?
Ok. Это разминка. Реальный вопрос заключается в том, как наилучшим образом приблизиться к многовариантному эквиваленту вышеизложенного.
В моем конкретном случае у меня есть N наблюдений двумерных данных для некоторой переменной y, которые я храню в матрице N x 2. (Для порядка, N составляет около 1000 в моем случае использования.)
Я считаю, что среднее значение каждого компонента каждого наблюдения равно 0, а стандартным значением каждого компонента каждого наблюдения является 1 (и я рад жестко закодировать их как таковые). Однако я полагаю, что корреляция (rho) варьируется от наблюдения к наблюдению как (простая) функция другой наблюдаемой переменной x (хранится в N элементном векторе). Например, мы можем сказать, что rho[n] = 2*inverse_logit(beta * x[n]) - 1 для n in 1:N, и наша цель - узнать о beta из наших данных. То есть ковариационная матрица для n-го наблюдения будет:
[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, окрашенная correls на левой панели и drivers на правой панели ... так что идея в том, что чем выше driver, тем ближе до 1 correl и чем ниже driver, тем ближе к -1 correl. Таким образом, можно было бы ожидать, что красные точки на левой панели будут «слева направо вверх», а синие точки быть «вверх-слева-вниз-вправо», и они действительно:
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()
Но я надеюсь, что кто-то может указать мне правильное направление для «маргинализованного» подхода, в котором мы разбиваем наш двумерный нормаль на пару независимых нормалей, которые можно смешать, используя корреляцию. Причина, по которой я нуждаюсь в этом, заключается в том, что в моем реальном случае использования оба аспекта имеют толстый хвост. Я счастлив смоделировать это как студенческий дистрибутив, но проблема в том, что STAN позволяет указывать только один nu (не один для каждого измерения), поэтому я думаю, что мне нужно будет найти способ декомпозировать multi_student_t в пару независимых student_t, так что я могу установить степени свободы отдельно для каждого измерения.