I’m building a hierarchical population model for cosmologically distributed objects (in redshift) where I do not know all of the distances. Thus, I need to sample from some redshift distribution for the objects with unknown distances.

The posterior looks like this:

propto_Pr(_psi_v.pdf (1.4 MB)

The term I need to sample is the dV/dOmega term.

It is obviously unnormalized because it gives a rate density. My naive approach is to design some functions that build this density and then sample from it’s log:

```
vector obj_rate(vector z, real rho, real zp, real p1, real p2){
/*
The rate in units on N/yr/Gpc^3
*/
vector[num_elements(z)] val;
// pow is not vectorized so
// we need to do a loop
for (i in 1:num_elements(z)){
val[i] = (1. + pow(z[i] / zp, p2));
}
// elementwise division
return rho * (1. + p1 * z) ./ val;
}
vector differential_comoving_volume(vector z){
/*
the differential comoving volume in units
of Gpc^3
*/
vector[num_elements(z)] td;
real dh;
dh = 4326.009494949495;
td = (comoving_transverse_distance(z)/3.086E24);
return (dh* td .* td ./ a_factor(z)) * 1E-9; // Gpc^3
}
vector C(vector z, real rho, real zp, real p1, real p2){
/*
The cosmological number
*/
return obj_rate(z, rho, zp, p1, p2) .* differential_comoving_volume(z) ./ (1+z);
}
real redshift_distribution_lpdf( vector z, real rho, real zp, real p1, real p2){
vector[num_elements(z)] prob;
real lprob;
prob = C(z,rho,zp,p1,p2);
lprob = sum(log(prob));
return lprob;
}
```

and then in my model, I sample with this statement:

```
parameters {
real L0_mu;
real<lower=log10(2.)> G0_mu;
real E0_mu;
real T0_mu;
real<lower=0> L0_sigma;
real<lower=0> G0_sigma;
real<lower=0> E0_sigma;
real<lower=0> T0_sigma;
vector<lower=0, upper=10>[N_unknown] z_unknown;
```

…

```
z_unknown ~ redshift_distribution(rho,zp,p1,p2);
```

Is it ok to sample from a distribution in this way? I’m worried about constants out from of the posterior, etc. Is my “log prob.” that I define for the redshift distribution really a valid way to think about the probability?

Cheers!