-
Notifications
You must be signed in to change notification settings - Fork 0
/
model2.stan
53 lines (51 loc) · 1.45 KB
/
model2.stan
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
// To add/fix:
// - GP for months
data {
int S; int T; int darkN;
real temp[T, S];
real lat[S]; real lon[S];
int month[T];
}
transformed data {
vector[S] zeroS;
for (i in 1:S) zeroS[i] <- 0;
}
parameters {
matrix[S, 12] baseline;
vector[12] trend;
vector<lower=0>[12] tau_month;
real trend_lat;
real trend_lon;
cholesky_factor_corr[S] LOmega;
vector<lower=0>[S] tau; // station (margin) sd's
vector[darkN] darkErr;
real<lower=0, upper=1> rho;
}
model {
matrix[S, S] LSigma[12];
matrix[T, S] TSerr;
real m; real decade; int i_dark;
i_dark <- 1;
for (t in 1:T) {
decade <- (t-T/2.0)/120.;
for (s in 1:S) {
if (temp[t, s]<5000) {
m <- baseline[s, month[t]]
+ trend[month[t]]/10*decade
+ trend_lat/100*lat[s]*decade
+ trend_lon/100*lon[s]*decade;
if (t>1) m <- m + rho*TSerr[t-1, s];
TSerr[t, s] <- temp[t, s] - m; }
else { TSerr[t, s] <- darkErr[i_dark]; i_dark <- i_dark+1; }
}}
for (i in 1:12) LSigma[i] <- diag_pre_multiply(tau*tau_month[i], LOmega);
for (i in 1:T) TSerr[i] ~ multi_normal_cholesky(zeroS, LSigma[month[i]]);
darkErr ~ normal(0, 30);
tau ~ lognormal(0.7, 1.0);
tau_month ~ lognormal(0.0, 1.0); // Fix one index?
LOmega ~ lkj_corr_cholesky(1.0);
for (i in 1:S) baseline[i] ~ normal(5, 50);
trend ~ normal(0, 30);
trend_lat ~ normal(0, 30);
trend_lon ~ normal(0, 30);
}