-
Notifications
You must be signed in to change notification settings - Fork 1
/
stan_code.txt
258 lines (189 loc) · 10 KB
/
stan_code.txt
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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
functions {
// This didn't turn out to be important when setting epilson > 0, but I include it here.
real my_normal_cdf_log(real x, real mu, real sig) {
real pull;
pull <- (x - mu)/sig;
if (pull > -25) {
return normal_cdf_log(x, mu, sig);
} else {
return 2.5/pull^4 - pull^-2. - 0.5*pull^2. -0.9189385332046 - log(-pull);
}
}
}
data {
int<lower=0> n_sne; // Number of SNe
int<lower=0> n_samples; // Number of SN samples
int<lower=0> n_calib; // Number of distance uncertainty systematics (e.g., zeropoints)
int n_x1c_star; // Number of redshift nodes per sample
int sample_list[n_sne]; // Which sample each SN is in
real <lower=0> redshifts[n_sne]; // The redshift for each SN. Union2.1 used z_CMB, but this could be improved
matrix [n_sne, n_x1c_star] redshift_coeffs; // How each SN depends on each node.
vector[3] obs_mBx1c [n_sne]; // SALT2 fits
matrix[3,3] obs_mBx1c_cov [n_sne]; // Covariance of SALT2 fits
matrix[3, n_calib] d_mBx1c_d_calib [n_sne]; // Sensitivity of each SN to each distance systematic
vector [n_sne] phighmass; // Probability of each SN being located in a high-mass host galaxy
int nzadd; // Integrating comoving distance using Simpson's rule. Add these redshifts to make sure that we have good sampling in redshift.
real redshifts_sort_fill [2*(n_sne + nzadd) - 1]; // Each SN redshift is a point for Simpson's rule. Fill in with intermediate points.
int unsort_inds[n_sne + nzadd]; // For converting to SN index.
real outl_mBx1c_uncertainties [3]; // Outlier distribution mB, x_1, c uncertainties
real outl_frac_prior_lnmean; // Log(fraction outliers) prior
real outl_frac_prior_lnwidth; // Log(fraction outliers) prior fractional width
real mB_cuts[n_samples]; // Magnitude limits for each sample
}
parameters {
real <lower = -20, upper = -18.> MB;
real <lower = -0.2, upper = 0.3> alpha_angle_low;
real <lower = -1.4, upper = 1.4> beta_angle_blue;
real <lower = -0.2, upper = 0.3> alpha_angle_high;
real <lower = -1, upper = 1.4> beta_angle_red;
real <lower = -0.2, upper = 0.3> delta_0;
real <lower = 0, upper = 1> delta_h;
real <lower = 0, upper = 1> Om;
real <lower=-1.3, upper = -0.3> log10_sigma_int[n_samples];
simplex [3] mBx1c_int_variance;
cholesky_factor_corr[3] Lmat;
real true_x1[n_sne];
real <lower = -1, upper = 2> true_c[n_sne];
matrix <lower = -1.5, upper = 1.5> [n_samples, n_x1c_star] x1_star;
matrix <lower = -0.4, upper = 0.4> [n_samples, n_x1c_star] c_star;
matrix <lower=-0.995, upper = 0.995> [n_samples, n_x1c_star] delta_c;
vector <lower= -0.5, upper = 0.5> [n_samples] log10_R_x1;
vector <lower=-1.5, upper = -0.5> [n_samples] log10_R_c;
vector [n_calib] calibs;
real <lower = 0, upper = 1> outl_frac;
}
transformed parameters {
vector [3] model_mBx1c [n_sne];
matrix [3,3] model_mBx1c_cov [n_sne];
vector [n_sne] mean_mB_by_SN;
vector [n_sne] var_mB_by_SN;
matrix [3,3] rho_int_mat;
matrix [3,3] mBx1c_int_variance_mat [n_samples];
vector [3] sig_int_vector;
real alpha;
real dalpha;
real beta;
real dbeta;
real alpha_eff;
real beta_eff;
real delta_eff;
vector [n_samples] R_x1;
vector [n_samples] R_c;
vector [n_sne] x1_star_by_SN;
vector [n_sne] R_x1_by_SN;
vector [n_sne] mu_c_by_SN;
vector [n_sne] sigma_c_by_SN;
vector [n_sne] alpha_c_by_SN;
real delta_c_eff;
real mean_c_eff;
real std_c_eff;
real Hinv_sort_fill [2*(n_sne + nzadd) - 1];
real r_com_sort[n_sne + nzadd];
real model_mu[n_sne];
vector [n_sne] outl_loglike;
vector [n_sne] PointPosteriors;
// -------------Begin numerical integration-----------------
for (i in 1: 2*(n_sne + nzadd) - 1) {
Hinv_sort_fill[i] <- 1./sqrt( Om*pow(1. + redshifts_sort_fill[i], 3) + (1. - Om) );
}
r_com_sort[1] <- 0.; // Redshift = 0 should be first element!
for (i in 2:(n_sne + nzadd)) {
r_com_sort[i] <- r_com_sort[i - 1] + (Hinv_sort_fill[2*i - 3] + 4.*Hinv_sort_fill[2*i - 2] + Hinv_sort_fill[2*i - 1])*(redshifts_sort_fill[2*i - 1] - redshifts_sort_fill[2*i - 3])/6.;
}
for (i in 1:n_sne) {
model_mu[i] <- 5.*log10((1. + redshifts[i])*r_com_sort[unsort_inds[i] + 1]) + 43.1586133146;
}
// -------------End numerical integration---------------
R_x1 <- exp(log(10.) * log10_R_x1);
R_c <- exp(log(10.) * log10_R_c);
// Convert from angles to slopes
alpha <- 0.5*(tan(alpha_angle_low) + tan(alpha_angle_high));
dalpha <- tan(alpha_angle_high) - tan(alpha_angle_low);
beta <- 0.5*(tan(beta_angle_blue) + tan(beta_angle_red));
dbeta <- tan(beta_angle_red) - tan(beta_angle_blue);
// Construct unexplained ("intrinsic") covariance matrix
rho_int_mat <- Lmat * Lmat';
model_mBx1c_cov <- obs_mBx1c_cov;
for (i in 1:n_samples) {
sig_int_vector[1] <- sqrt(mBx1c_int_variance[1])*pow(10, log10_sigma_int[i]); // This vector is in dispersion, not variance
sig_int_vector[2] <- sqrt(mBx1c_int_variance[2])*pow(10, log10_sigma_int[i])/0.13;
sig_int_vector[3] <- sqrt(mBx1c_int_variance[3])*pow(10, log10_sigma_int[i])/(-3.);
mBx1c_int_variance_mat[i] <- rho_int_mat .* (sig_int_vector * sig_int_vector');
}
for (i in 1:n_sne) {
// Broken-linear color and shape correction.
// If statements are ugly, but they don't seem to be any slower than using a hyperbolic interpolant.
if (true_x1[i] > 0) {
alpha_eff <- alpha + dalpha/2;
} else {
alpha_eff <- alpha - dalpha/2;
}
if (true_c[i] > 0) {
beta_eff <- beta + dbeta/2;
} else {
beta_eff <- beta - dbeta/2;
}
// Redshift-depedent host-mass correction
delta_eff <- delta_0*(( 1.9*(1 - delta_h)/(1 + 0.9*exp(0.95*log(10.)*redshifts[i])) ) + delta_h);
// Building the model of the observations
model_mBx1c[i][1] <- -(alpha_eff*true_x1[i] - beta_eff*true_c[i] - MB - model_mu[i] + delta_eff*phighmass[i]);
model_mBx1c[i][2] <- true_x1[i];
model_mBx1c[i][3] <- true_c[i];
model_mBx1c[i] <- model_mBx1c[i] + d_mBx1c_d_calib[i] * calibs;
// Building the covariance model
model_mBx1c_cov[i] <- obs_mBx1c_cov[i] + mBx1c_int_variance_mat[sample_list[i]];
// Builing the shape and color population model
x1_star_by_SN[i] <- dot_product(x1_star[sample_list[i]], redshift_coeffs[i]);
R_x1_by_SN[i] <- R_x1[sample_list[i]];
mean_c_eff <- dot_product(c_star[sample_list[i]], redshift_coeffs[i]);;
std_c_eff <- R_c[sample_list[i]];
delta_c_eff <- dot_product(delta_c[sample_list[i]], redshift_coeffs[i]);
// Remember: mean, std, and delta are properties of the distribution
// mu, sigma, alpha are parameters in the skew normal
mu_c_by_SN[i] <- mean_c_eff - 1.414*std_c_eff*delta_c_eff/sqrt(3.14159 - 2*delta_c_eff^2);
sigma_c_by_SN[i] <- std_c_eff/sqrt(1. - 0.6366*delta_c_eff^2);
alpha_c_by_SN[i] <- delta_c_eff/sqrt(1 - delta_c_eff^2);
// Building the selection efficiency model
mean_mB_by_SN[i] <- MB + model_mu[i] - alpha*x1_star_by_SN[i] - dalpha*R_x1_by_SN[i]/2.50663
+ beta*mean_c_eff + dbeta*R_c[sample_list[i]]/2.50663
- delta_eff*0.5; //No calibration in this term, as magnitude limits have same calibration uncertainty as distances!
var_mB_by_SN[i] <- model_mBx1c_cov[i][1,1] + 0.25
+ 0.797885*beta*dbeta*mean_c_eff*R_c[sample_list[i]]
+ (beta*R_c[sample_list[i]])^2 + 0.0908451*(dbeta*R_c[sample_list[i]])^2
+ 0.797885*alpha*dalpha*x1_star_by_SN[i]*R_x1_by_SN[i]
+ (alpha*R_x1_by_SN[i])^2 + 0.0908451*(dalpha*R_x1_by_SN[i])^2
+ 0.25*delta_eff^2;
// Building the relative log likelihood for the SN being an outlier
outl_loglike[i] <- log(1 - outl_frac)
+ multi_normal_log(obs_mBx1c[i], model_mBx1c[i], model_mBx1c_cov[i])
+ my_normal_cdf_log(mB_cuts[sample_list[i]] + d_mBx1c_d_calib[i][1] * calibs, obs_mBx1c[i][1], 0.5)
- log(normal_cdf(mB_cuts[sample_list[i]], mean_mB_by_SN[i], sqrt(var_mB_by_SN[i])) + 0.01) //No calibration in this term, see above!
- (log(outl_frac)
+ normal_log(obs_mBx1c[i][1], model_mBx1c[i][1], outl_mBx1c_uncertainties[1])
+ normal_log(obs_mBx1c[i][2], model_mBx1c[i][2], outl_mBx1c_uncertainties[2])
+ normal_log(obs_mBx1c[i][3], model_mBx1c[i][3], outl_mBx1c_uncertainties[3])
);
// Finally, the likelihood for each SN
PointPosteriors[i] <- log_sum_exp(log(1 - outl_frac)
+ multi_normal_log(obs_mBx1c[i], model_mBx1c[i], model_mBx1c_cov[i])
+ my_normal_cdf_log(mB_cuts[sample_list[i]] + d_mBx1c_d_calib[i][1] * calibs, obs_mBx1c[i][1], 0.5)
- log(normal_cdf(mB_cuts[sample_list[i]], mean_mB_by_SN[i], sqrt(var_mB_by_SN[i])) + 0.01) //No calibration in this term, see above!
, log(outl_frac)
+ normal_log(obs_mBx1c[i][1], model_mBx1c[i][1], outl_mBx1c_uncertainties[1])
+ normal_log(obs_mBx1c[i][2], model_mBx1c[i][2], outl_mBx1c_uncertainties[2])
+ normal_log(obs_mBx1c[i][3], model_mBx1c[i][3], outl_mBx1c_uncertainties[3])
);
}
}
model {
increment_log_prob(sum(PointPosteriors));
true_x1 ~ normal(x1_star_by_SN, R_x1_by_SN);
true_c ~ skew_normal(mu_c_by_SN, sigma_c_by_SN, alpha_c_by_SN);
for (i in 1:n_samples) {
x1_star[i] ~ cauchy(0, 1);
c_star[i] ~ cauchy(0, 0.3);
}
Lmat ~ lkj_corr_cholesky(1.0);
calibs ~ normal(0, 1);
outl_frac ~ lognormal(outl_frac_prior_lnmean, outl_frac_prior_lnwidth);
}