Im trying to estimate some parameters (rc and tB0) for this model, but i dont know why RStan says that logprob is Nan, i did the same problem in emcee with no issues.
Could you help me to see the error? Thank you so much!
My program:
functions {
real tBB(real r, // Comovil Position
real tb0, // initial time
real rc // Void
) {
return (tb0)*exp(-(r/rc)^4.0);
}
real tBBprima(real r, // Comovil Position
real tb0, // initial time
real rc // Void
) {
real t0=13.7;
return (-4.0)*tBB(r,tb0,rc)*((r^3.0)/(rc^4.0));
}
real invmpom(real r, // Comovil Position
real tb0, // initial time
real rc // Void
) {
real t0 = 13.7;
real t1 = (2.0/9.0)^(-1.0/3.0);
real t2 = (t0 - tBB(r,tb0,rc))^(5.0/3.0);
real t3 = 3.0*(t0 - tBB(r,tb0,rc))+2.0*r*tBBprima(r, tb0, rc);
return (t1*t2)/t3;
}
real R(real r, // Comovil Position
real t, // Comovil Time
real tb0, // initial time
real rc // Void
) {
real t0=13.7;
return r*(((t-tBB(r,tb0,rc))/(t0-tBB(r,tb0,rc)))^(2.0/3.0));
}
real Nr(real r, // Comovil Position
real t, // Comovil Time
real tb0, // initial time
real rc // Void
) {
real t0=13.7;
real t1 = (8.0/3.0)*(r/rc)^4.0;
real t2 = tBB(r,tb0,rc)/(t0 - tBB(r,tb0,rc));
real t4 = tBB(r,tb0,rc)/(t - tBB(r,tb0,rc));
real t3 = 1.0/(1.0 - t1*t2);
return t1*t4*t3;
}
real dL(real z, // Redshift
real r, // Comovil Position
real t, // Comovil Time
real tb0, // initial time
real rc // Void
)
{
return ((1+z)^2.0)*R(r,t,tb0,rc);
}
real Ev(real z // Redshift
) {
return 0;
}
real mb(real z, // Redshift
real r, // Comovil Position
real t, // Comovil Time
real tb0 , // Local Hubble Rate
real rc
) {
return 5*log10(dL(z,r,t,tb0,rc))+Ev(z)-25.0;
}
real Errorsum(real[] e,int N){
real suma=0;
for (i in 1:N){
suma+=1/(e[i]*e[i]);
}
return suma;
}
real MargFactor(real[] z,real[] m,real[] r, real[] t, real[] e,int N,real tb0,real rc){
real suma=0;
for (i in 1:N){
suma+=(5*log10(dL(z[i],r[i],t[i],tb0,rc))-m[i])/(e[i]*e[i]);
}
return suma;
}
real [] dgz(real z, // Redshift
real r, // Comovil Position
real t, // Comovil Time
real tb0 , // Local Hubble Rate
real rc
) {
real cc = 0.306595;
real z1 = 1.0 + z;
real c1 = (9.0/2.0)^(-1.0/3.0);
real c2 = t - tBB(r,tb0,rc);
real c3 = 1.0 - Nr(r,t,tb0,rc)/2.0;
real f1 = (((c2)^(1.0/3.0))*invmpom(r,tb0,rc))/c3;
return { -1.5*(1.0+Nr(r,t,tb0,rc))*c2/(c3*z1), (9.0/2.0)*cc*c1*f1/z1 };
}
real[] GeodesicEquation(real z, // Redshift
real[] c, // Comovil Coordinates
real[] theta, // Parameters
real[] x_r, // unused data
int[] x_i
) {
real r=c[1];
real t=c[2];
real tb0 = theta[1];
real rc = theta[2];
return dgz(z,r,t,tb0,rc);
}
}
data {
int<lower = 0> N; // number of measurement
real zobs[N]; // Redshifts
real<lower = 0> mobs[N]; // Supernovaes magnitude
real<lower = 0> errorm[N]; // Supernovaes error
}
parameters {
real<lower = 0.1> tb0; // Initial time
real<lower = 0.1> rc; // Void size parameter
}
transformed parameters {
real t0=13.7;
real ComovilCoordinates[N,2]
= integrate_ode_rk45(GeodesicEquation, {t0,0.0}, 0, zobs, {tb0,rc},
rep_array(0.0, 0), rep_array(0, 0),
1e-5, 1e-3, 5e2);
}
model {
tb0 ~ normal(3, 10);
rc ~ normal(2,2);
//sigma ~ lognormal(log(10), 1);
//alpha ~ lognormal(log(10), 1);
//beta ~ lognormal(log(10), 1);
for(i in 1:N){
mobs[i] ~ normal(mb(zobs[i],ComovilCoordinates[i,2],ComovilCoordinates[i,1],tb0,rc), errorm[i]);
}
target+=0.5*(MargFactor(zobs, mobs, ComovilCoordinates[,2],ComovilCoordinates[,1], errorm, N,tb0,rc)*MargFactor(zobs, mobs, ComovilCoordinates[,2],ComovilCoordinates[,1], errorm, N,tb0,rc))/(Errorsum( errorm, N))-log(Errorsum(errorm,N)/(2*pi()));
}
And the error is something like:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: normal_lpdf: Location parameter is nan, but must be finite! (in 'string', line 59, column 1 to column 10)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability ...