Обзор модели оценки параметров полиномиальной регрессии Стэна

У меня есть следующая модель полиномиальной регрессии:

Версия изображения

Версия LaTeX

$ Y_i | \ mu_i, \ sigma ^ 2 \ sim \ text {Normal} (\ mu_i, \ sigma ^ 2), i = 1, \ dots, n \ \ text {независимый} $

$ \ mu_i = \ alpha + \ 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 \ бар {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 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. Это признак того, что я допустил ошибку? Как я уже сказал, я просмотрел весь свой код и математику, поэтому, насколько я могу судить, все кажется в порядке.