Несколько параметров нашей модели не сходятся в WinBUGS

Question

Я установил бейесовскую модель, которая включает некоторые уравнения из литературы. Я хочу одновременно смоделировать десять параметров (P, R, Rear, K, R20, mu, a, b, tau и sigma) в соответствии с данными мониторинга. Модель синтаксически корректна, загружает данные и хорошо компилируется. Я обновил 1000000. Он может дать результаты каждого параметра.

Но три параметра (mu, a и b) нашей модели, похоже, всегда не сходятся в WinBUGS. Когда я только симулировал эти три параметра, они также не смогли сходиться.

Проблема в том, что я не знаю, как это настроить. Я пытался изменить распределение pior некоторых параметров много раз. Но бесполезно. Я также не знаю, как оценить модель, потому что результаты DIC показывают «неопределенные реальные результаты».

Я искал наш веб-сайт и обнаружил, что был также похожий вопрос.

Я получил следующий ответ: «WinBUGS уже несколько лет, и прошло 8 лет с момента его последнего обновления! Я думаю, вы должны забыть об этом, потому что есть несколько лучших альтернатив. Вы можете попробовать практически один и тот же код в Jags или Stan, в котором оба могут использоваться в R через rJags и RStan. Stan особенно важен, потому что он использует MCMC, который сходится во многих ситуациях, чего нет в WinBUGS. " Сетевой метаанализ WinBUGS Weibull

Однако я также новичок в Джагсе или Стене. Я не знаю, как изменить его на код STAN.

Буду признателен всем за предложения о том, как настроить это (структуру модели, уравнение, распределение пор, начальные значения и т. Д.), Чтобы он работал правильно. Кроме того, я стараюсь делать анализы с использованием R.

Заранее спасибо.

Вот моя модель

    model; 
{ 
for( i in 1 : N ) {
lambda[i] ~ dnorm(u[i],tau)  }

for( i in 1 : N ) {
u[i] <- P[i] - R[i] + Rear[i]
P[i] <- mu * IZ[i]  
IZ[i] <- I[i] * exp(e[i])  
e[i] <- (-1) * Kd[i]
Kd[i] <- a + b * log(Tur[i])

R[i] <- R20 * pow(1.047,T[i] - 20)
Rear[i] <- K * (Os[i] - O[i])
Os[i] <- 14.62 - 0.3671 * T[i] + (0.004497 * T[i]) * T[i] - 0.966 * Sa[i] + (0.00205 * Sa[i]) * T[i] + (2.739E-4 * Sa[i]) * Sa[i]

T[i] ~ dnorm( 0.0,1.0E-6)
I[i] ~ dnorm( 0.0,1.0E-6) 
Tur[i] ~ dnorm(10, 0.1)I(0,)  
Sa[i] ~ dnorm( 0.0,1.0E-6)
O[i] ~ dnorm( 0.0,1.0E-6)  
  }


tau ~ dgamma(0.001,0.001)
R20 ~ dnorm( 0.0,1.0E-6)
K ~ dnorm( 0.0,1.0E-6)
mu ~ dnorm( 0.0,1.0E-6)
b ~ dnorm( 0.0,1.0E-6)
a ~ dnorm( 0.0,1.0E-6)
sigma <- 1 / sqrt(tau) 
} 

Вот мои данные

list(lambda=c(0.3, 0, 0.03, 0.12,   0.13,   0.12,   0.03,   0.27,   0.29,   0.02,   0.2,    0.25,   0.26,   0.15,   0.16,   0.74,   0.3,    0.4,    0.4,    0.28,   0.15,   -0.15,  -0.07,  0.02,   -0.13,  -0.3,   -0.26,  -0.36,  -0.26,  -0.28,  -0.26,  -0.32,  -0.18,  -0.29,  -0.27,  -0.09,  -0.32,  -0.21,  -0.18,  -0.16,  -0.23,  -0.18,  -0.16,  -0.13,  -0.18,  -0.07,  -0.15,  -0.11,  -0.03,  0.01,   0.03,   0.12,   -0.07,  0.12,   0.08,   0.18,   0.24,   0.3,    0.08,   0.35,   0.27,   0.29,   0.02,   0.2,    0.25,   0.26,   0.24,   0.2,    0.2,    -0.09,  -0.2,   -0.16,  -0.16,  -0.08,  -0.06,  -0.17,  -0.32,  -0.14,  -0.18,  -0.32,  -0.16,  -0.28,  -0.04,  -0.27,  -0.15,  -0.06,  -0.15,  -0.11,  -0.15,  -0.11,  -0.03,  -0.25,  -0.02,  -0.06,  -0.16,  -0.08,  0.08),
  T=c(27.45,27.56,27.33,27.39,27.61,27.61,28,28.4,28.78,29,29.5,29.94,30,30.22,30.56,31,31.4,31.8,31.9,30.82,30.15,30.25,30.25,29.86,29.94,29.76,29.52,29.26,29.02,28.8,28.59,28.33,28.09,27.85,27.63,27.46,27.31,27.17,27.04,26.91,26.77,26.63,26.52,26.41,26.29,26.13,26.01,25.88,25.79,25.79,25.77,25.75,25.76,25.86,25.8,25.85,26.01,26.15,26.27,26.49,26.74,27,27.56,27.87,27.83,28.02,28.36,28.67,29.11,29.45,29.55,29.35,28.67,28.52,28.05,27.98,27.91,27.55,27.26,26.99,26.72,26.46,26.21,25.98,25.71,25.44,25.23,25.01,24.79,24.6,24.41,24.23,24.05,23.88,23.73,23.59,23.48),
I=c(24,52,141,138,179,193,148,306,210,423,306,258,235,385,157,382,598,517,427,603,279,196,235,124,132,87,42,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,20,24,52,141,138,179,193,148,306,210,423,306,258,235,385,157,382,598,517,427,603,279,196,235,124,132,87,42,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,20,59,90,109,112,86,143,206,248,244,260,384,361,460,563,398,204,510,614,608,588,523,430,332,246,161,73,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,22,90),
 Sa=c(2.67,2.67,2.68,2.68,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.68,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.68,2.68,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.68,2.67,2.67,2.67,2.68,2.67,2.67,2.67,2.68,2.67,2.67,2.67,2.68,2.67,2.68,2.67,2.68,2.67,2.68,2.71,2.68,2.6,2.68,2.68,2.68,2.68,2.69,2.69,2.69,2.69,2.68,2.68,2.69,2.68,2.7,2.69,2.69,2.695,2.697,2.699,2.701,2.703,2.69,2.7,2.7,2.69,2.7,2.7,2.69,2.7,2.69,2.7,2.7,2.7,2.7,2.7,2.7,2.7,2.7,2.69,2.69,2.7,2.7),
 O=c(6.23,6.53,6.53,6.56,6.68,6.81,6.93,6.96,7.23,7.52,7.54,7.74,7.99,8.25,8.4,8.56,9.3,9.6,10,10.4,10.68,10.83,10.68,10.61,10.63,10.5,10.2,9.94,9.58,9.32,9.04,8.78,8.46,8.28,7.99,7.72,7.63,7.31,7.1,6.92,6.76,6.53,6.35,6.19,6.06,5.88,5.81,5.66,5.55,5.52,5.53,5.56,5.68,5.61,5.73,5.81,5.99,6.23,6.53,6.61,6.96,7.23,7.52,7.54,7.74,7.99,8.25,8.49,8.69,8.89,8.8,8.6,8.44,8.28,8.2,8.14,7.97,7.65,7.51,7.33,7.01,6.85,6.57,6.53,6.26,6.11,6.05,5.9,5.79,5.64,5.53,5.5,5.25,5.23,5.17,5.01,4.93),   

N=97)

list(mu=0.5, a=1, b=2, K=0.5, R20=1, tau=1) 

Часть результатов

ode statistics 
node     mean    sd  MC error   2.5%    median  97.5%   start   sample
mu  713.5   620.3   17.57   0.07538 590.0   2162.0  4001    996000   

node     mean    sd  MC error   2.5%    median  97.5%   start   sample
a   8.924   2.585   0.08136 0.6363  9.707   11.45   4001    996000

временной ряд итерации

введите описание изображения здесь введите описание изображения здесь

Краткий обзор переменных, введенных в модель:

T : water temperature
I: solar irradiance at water surface
Sa: water salinity
O: disolved oxygen

(P, R, Rear, K, R20, mu, a, b, tau, and sigma)
P: primary prodution 
R: ecosystem respiration
Rear: molecular diffusivity of dissolved oxygen across the air–water interface
K: gas transfer velocity parmeter
R20: respiration at T= 20℃
mu: the slope of the photosynthesis–irradiance relationship at low light conditions
a and b: regression coefficients between turbidity and irradiance attenuation coefficient
tau: the precision?
sigma: the standard error?

Answer 1

Было бы написано что-то вроде этого на языке Стэн, я думаю

data {
  int<lower=1> N;
  vector[N] lambda;
  vector[N] T;
  vector[N] I;
  vector[N] Sa;
  vector[N] O;
}
parameters {
  vector<lower=0>[N] Tur;
  real<lower=0> sigma;
  real R20;
  real K;
  real mu;
  real a;
  real b;
}
model {
  vector[N] u;
  for (i in 1:N) {
    real Os = 14.62 - 0.3671 * T[i] + (0.004497 * T[i]) * T[i] - 0.966 * Sa[i] + (0.00205 * Sa[i]) * T[i] + (2.739E-4 * Sa[i]) * Sa[i];
    real Rear = K * (Os - O[i]);
    real R = R20 * 1.047^(T[i] - 20);
    real Kd = a + b * log(Tur[i]);
    real e = -1 * Kd;
    real IZ = I[i] * exp(e);
    real P = mu * IZ;
    u[i] = P - R + Rear;
  }
  target += normal_lpdf(lambda | u, sigma);

  target += normal_lpdf(Tur | 10, 3.162278); // do not need to truncate
  // other priors in the form of target += distribution_lpdf( | , );
}

Однако в вашей модели BUGS есть несколько запутанных вещей. Во-первых, у вас есть операторы выборки для T, I, Sa и O, даже если это явно данные, а операторы выборки содержат только константы. Константы не имеют отношения к предложениям параметров в Stan или тем, какие предложения принимаются, так что вы можете просто опустить все это. Если вам нужны случайные их отрисовки для прогнозирования в какой-то большей популяции, вы можете сделать это в Stan с блоком generated quantities в конце вашей программы Stan или позже в R.

Во-вторых, эти приоры, которые вы используете - хотя и распространены в моделях BUGS, - противоположны выводам. Для модели, подобной этой (и многих других моделей), вам придется выразить то, во что вы действительно верите об этих параметрах, вместо того, чтобы говорить пустые вещи, например, есть вероятность того, что mu будет между -1000 и 1000 Существует множество рекомендуемых приоров для Стэна. Кроме того, имейте в виду, что параметризация нормального распределения в Stan основана на ожидании и стандартном отклонении (а не на точности).

См. Также Приложение B пользователя Stan , руководство и FAQ . Если вы используете гораздо лучшие приоры, это может сработать, но обратите внимание на предупреждающие сообщения, которые вы, вероятно, получите при первой попытке (которые указывают на то, что заданное распределение, которое вы определили, не способствует выборке из).