У меня есть следующая модель полиномиальной регрессии:
Версия изображения
LaTeXВерсия
$ Y_i |\ mu_i, \ sigma ^ 2 \ sim \ text {Normal} (\ mu_i, \ sigma ^ 2), i = 1, \ dots, n \ \ text {independent} $
$ \ mu_i = \альфа + \ beta_1 x_ {i1} + \ beta_2 x_ {i2} + \ beta_3 x_ {i1} ^ 2 + \ beta_4 x_ {i2} ^ 2 + \ beta_5 x_ {i1} x_ {i2} $
$ \ alpha \ sim \ text {некоторые подходящие предыдущие} $
$ \ beta_1, \ dots, \ beta_5 \ sim \ text {некоторые подходящие априоры} $
$ \ sigma ^ 2\ sim \ text {некоторые подходящие предыдущие} $
Я хочу взять в качестве входных данных размер выборки и векторы наблюдений для $ y_i $, $ x_ {i1} $ и $ x_ {i2} $.Код для этого следующий:
data{
int<lower=1> n;
vector[n] x1;
vector[n] x2;
vector[n] y;
}
Я хочу стандартизировать (центрировать и масштабировать) две входные переменные, чтобы получить стандартизированные переменные регрессора x1_std и x2_std.Код для этого находится в блоке transformed data следующим образом:
transformed data{
real bar_x1;
real x1_sd;
vector[n] x1_std;
real bar_x2;
real x2_sd;
vector[n] x2_std;
real y_sd;
bar_x1 = mean(x1);
x1_sd = sd(x1);
x1_std = (x1 - bar_x1)/x1_sd; // centered and scaled
bar_x2 = mean(x2);
x2_sd = sd(x2);
x2_std = (x2 - bar_x2)/x2_sd; // centered and scaled
y_sd = sd(y);
}
Затем я хочу подогнать вышеприведенную модель полиномиальной регрессии, используя стандартизированные переменные регрессора, и вернуть оценки для параметров регрессии $ \ alpha $., $ \ beta_1 $ и $ \ dots, \ beta_5 $, по и по оригинальной и стандартизированной шкале.
Исходя из этого, если я не ошибаюсь, формулы преобразования из стандартизированных параметров в исходную шкалу выглядят следующим образом:
Версия изображения
LaTeX версия
$ \ alpha = \ tilde {\ alpha} - \dfrac {\ gamma_1} {s_1} \ bar {x} _1 - \ dfrac {\ gamma_2} {s_2} \ bar {x} _2 + \ dfrac {\ gamma_3} {s_1 ^ 2} \ bar {x} _1 ^ 2+ \ dfrac {\ gamma_4} {s_2 ^ 2} \ bar {x} _2 ^ 2 + \ dfrac {\ gamma_5} {s_1 s_2} \ bar {x} _1 \ bar {x} _2 $
$ \ beta_1 = \ left (\ dfrac {\ gamma_1} {s_1} - 2 \ dfrac {\ gamma_3} {s_1 ^ 2} \ bar {x} _1 - \ dfrac {\ gamma_5} {s_1 s_2} \ bar {x} _2 \ right) $
$ \ beta_2 = \ left (\ dfrac {\ gamma_2} {s_2} - 2 \ dfrac {\ gamma_4} {s_2 ^ 2} \ bar {x} _2 - \ dfrac {\ gamma_5} {s_1 s_2} \ bar {x} _1 \ right) $
$ \ beta_3 = \ dfrac {\ gamma_3} {s_1 ^ 2} $
$ \ beta_4 = \ dfrac {\ gamma_4} {s_2 ^ 2} $
$ \ beta_5 = \ dfrac {\ gamma_5} {s_1 s_2} $
Код, реализующий это, содержится в блоке generated quantities следующим образом:
alpha = alpha_std - beta1_std*bar_x1/x1_sd - beta2_std*bar_x2/x2_sd
+ (beta3_std*bar_x1^2)/x1_sd^2 + (beta4_std*bar_x2^2)/x2_sd^2
+ (beta5_std*bar_x2*bar_x1)/(x1_sd*x2_sd);
beta1 = beta1_std/x1_sd - 2*beta3_std*bar_x1/x1_sd^2
- beta5_std*bar_x2/(x1_sd*x2_sd);
beta2 = beta2_std/x2_sd - 2*beta4_std*bar_x2/x2_sd^2
- beta5_std*bar_x1/(x1_sd*x2_sd);
beta3 = beta3_std/x1_sd^2;
beta4 = beta4_std/x2_sd^2;
beta5 = beta5_std/(x1_sd*x2_sd);
Вся моя модель выглядит следующим образом:
data{
int<lower=1> n;
vector[n] x1;
vector[n] x2;
vector[n] y;
}
transformed data{
real bar_x1;
real x1_sd;
vector[n] x1_std;
real bar_x2;
real x2_sd;
vector[n] x2_std;
real y_sd;
bar_x1 = mean(x1);
x1_sd = sd(x1);
x1_std = (x1 - bar_x1)/x1_sd; // centered and scaled
bar_x2 = mean(x2);
x2_sd = sd(x2);
x2_std = (x2 - bar_x2)/x2_sd; // centered and scaled
y_sd = sd(y);
}
parameters{
real<lower=0> sigma;
real alpha_std;
real beta1_std;
real beta2_std;
real beta3_std;
real beta4_std;
real beta5_std;
}
transformed parameters {
real mu[n];
for(i in 1:n) {
mu[i] = alpha_std + beta1_std*x1_std[i]
+ beta2_std*x2_std[i] + beta3_std*x1_std[i]^2
+ beta4_std*x2_std[i]^2 + beta5_std*x1_std[i]*x2_std[i];
}
}
model{
alpha_std ~ normal(0, 10);
beta1_std ~ normal(0, 2.5);
beta2_std ~ normal(0, 2.5);
beta3_std ~ normal(0, 2.5);
beta4_std ~ normal(0, 2.5);
beta5_std ~ normal(0, 2.5);
sigma ~ exponential(1 / y_sd);
y ~ normal(mu, sigma);
}
generated quantities {
real alpha;
real beta1;
real beta2;
real beta3;
real beta4;
real beta5;
alpha = alpha_std - beta1_std*bar_x1/x1_sd - beta2_std*bar_x2/x2_sd
+ (beta3_std*bar_x1^2)/x1_sd^2 + (beta4_std*bar_x2^2)/x2_sd^2
+ (beta5_std*bar_x2*bar_x1)/(x1_sd*x2_sd);
beta1 = beta1_std/x1_sd - 2*beta3_std*bar_x1/x1_sd^2
- beta5_std*bar_x2/(x1_sd*x2_sd);
beta2 = beta2_std/x2_sd - 2*beta4_std*bar_x2/x2_sd^2
- beta5_std*bar_x1/(x1_sd*x2_sd);
beta3 = beta3_std/x1_sd^2;
beta4 = beta4_std/x2_sd^2;
beta5 = beta5_std/(x1_sd*x2_sd);
}
Я использую hills набор данных из пакета R's MASS:
library(MASS)
hills[18, 3] <- 18.65 # Fixing transcription error
x1 <- hills$dist
x2 <- hills$climb
y <- hills$time
n <- length(x1)
data.in <- list(x1 = x1, x2 = x2, y = y, n = n)
model.fit <- sampling(example, data.in)
А теперь явыводить стандартизированную (alpha_std, beta1_std, beta2_std, beta3_std, beta4_std, beta5_std) и оригинальную шкалу (alpha, beta1, beta2, beta3, beta4, beta5) параметры регрессии:
print(model.fit, pars = c("alpha_std", "alpha", "beta1_std", "beta2_std", "beta3_std", "beta4_std", "beta5_std", "beta1", "beta2", "beta3", "beta4", "beta5", "sigma"), probs = c(0.05, 0.5, 0.95), digits = 5)
Я был бы очень признателен, если бы люди могли посоветоватьправильно ли я поступил по этому поводу.Я также дважды и трижды проверил математику, поэтому я думаю это должно быть правильно.Несмотря на это, меня беспокоит то, что beta4 составляет 0,00000.Это признак того, что я допустил ошибку?Как я уже сказал, я просмотрел весь свой код и математику, поэтому, насколько я могу судить, все кажется нормальным ...