# Rejecting initial value:

**URL:** <https://discourse.mc-stan.org/t/rejecting-initial-value/27191>\
**Category:** Modeling\
**Tags:** initialization, fitting-issues\
**Created:** [April 19, 2022, 2:57am UTC](https://discourse.mc-stan.org/t/rejecting-initial-value/27191 "2022-04-19T02:57:20Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![chI-No](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/chi-no/32/12632_2.png) [@chI-No](https://discourse.mc-stan.org/u/chI-No)\
**Post date:** [April 19, 2022, 2:57am UTC](https://discourse.mc-stan.org/t/rejecting-initial-value/27191/1 "2022-04-19T02:57:20Z")

</div>

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:

```nohighlight
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:

```nohighlight
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 ...
```

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [April 19, 2022, 3:16am UTC](https://discourse.mc-stan.org/t/rejecting-initial-value/27191/2 "2022-04-19T03:16:51Z")

</div>

The key part of the error is:

```nohighlight
Chain 1: Exception: normal_lpdf: Location parameter is nan, but must be finite!

```

Which refers to your likelihood here:

```stan
  for(i in 1:N){
    mobs[i] ~ normal(mb(zobs[i],ComovilCoordinates[i,2],ComovilCoordinates[i,1],tb0,rc), errorm[i]); 
  }

```

The error indicates that the value being returned from the `mb()` is `NaN`, which means that either an uninitialised or an invalid input is being used somewhere within the function. I’d recommend using `print()` statements within the stan code to check the arguments being passed, as well as the values being returned within each of the functions that `mb` calls internally.

Additionally, you could use `expose_stan_functions` in `rstan` for testing the individual functions as well

---

<div class="post-metadata">

**Author:** ![chI-No](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/chi-no/32/12632_2.png) [@chI-No](https://discourse.mc-stan.org/u/chI-No)\
**Post date:** [April 21, 2022, 3:43am UTC](https://discourse.mc-stan.org/t/rejecting-initial-value/27191/3 "2022-04-21T03:43:45Z")

</div>

HI! thank you for the response. I put many print statements to see what happens and this are some results that i dont understand:

**1)** The initial value of r in the beggining of the run is giving NAN. This is completely strange to me, cause before the programm call any function the parameter r that enter to the function is already NAN. I could fix this changing the initial condition 0.0 in:

```nohighlight
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);
}

```

to something like 0.00001:

```nohighlight
transformed parameters {
  real t0=13.7;
  real ComovilCoordinates[N,2]
    = integrate_ode_rk45(GeodesicEquation, {t0,0.00001}, 0, zobs, {tb0,rc},
                         rep_array(0.0, 0), rep_array(0, 0),
                         1e-5, 1e-3, 5e2);
}

```

The code now run, but i dont understand why it throws me a NAN value before being 0.0. I have another similar code for other cosmological models with the same 0.0 conditions and it works very well. Also i wish to put 0.0 initial condition cause it seems that some models are very sensitive to little variations around 0.

**2)** Even if the code works, it doesn’t give me reasonable values (speaking in terms of literature and my other samples in Emcee and Zeus for python) and i found that when the code runs, it just throw all the information away cause interpret the first function of the code just to 0 all the time:

```nohighlight
functions {

  real tBB(real r, // Comovil Position
             real tb0, // initial time
             real rc // Void
) {
    print("r:",r);
    print("rc:",rc);
    print("(r/rc)^4.0):",(r/rc)^4.0);
    print("tb0:",tb0); 
  	print("tBB:",(tb0)*exp(-(r/rc)^4.0));
		return (tb0)*exp(-(r/rc)^4.0);
  }

```

The prints output are something like:

```nohighlight
...
tBB:0
r:13.4423
rc:0.122662
(r/rc)^4.0):1.4423e+08
tb0:0.141183
tBB:0
r:13.4423
rc:0.122662
(r/rc)^4.0):1.4423e+08
tb0:0.141183
tBB:0
r:13.4423
...

```

So always tBB act like 0. There is something that i can do to the program doesnt take little values to 0? Cause i am losing all the information.

Thank you for your ncie help

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [April 21, 2022, 4:05am UTC](https://discourse.mc-stan.org/t/rejecting-initial-value/27191/4 "2022-04-21T04:05:35Z")

</div>

The issue here is that the calculation of `tBB` is underflowing, because the result is so small that it can’t be distinguished from 0.

Using the printed values above to illustrate, `tBB` is calculated as:

```nohighlight
tBB = (tb0)*exp(-(r/rc)^4.0)

```

Where:

```nohighlight
tb0 = 0.141183
(r/rc)^4.0 = 1.4423e+08

```

Substituting those values into the calculation of `tBB`:

```nohighlight
tBB = 0.141183 * exp(-1.4423e+08)

```

The underflow will first occur here due to the call to `exp`, since `exp` will [always return 0](https://en.cppreference.com/w/cpp/numeric/math/exp#Notes) for inputs \< -708.

I believe you may need to rewrite your likelihood/calculations to the log scale, so that you’re not exponentiating and multiplying such massively large (either positive or negative) values

---

<div class="post-metadata">

**Author:** ![chI-No](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/chi-no/32/12632_2.png) [@chI-No](https://discourse.mc-stan.org/u/chI-No)\
**Post date:** [April 21, 2022, 5:31am UTC](https://discourse.mc-stan.org/t/rejecting-initial-value/27191/5 "2022-04-21T05:31:21Z")

</div>

Oh, thats unfortunate. Is almost impossible rewrite all in log scale cause the math between this relations is very complicated and tbb is related to a sum between a variable and a really complex integral. I think i’m gonna pass to use STAN in this opportunity.

If you know a solution that doesn’t involve modify all the physics please i’m gonna stay tune! Thank you so much for your time!

I’m gonna come back

---

<div class="post-metadata">

**Author:** ![nhuurre](https://avatars.discourse-cdn.com/v4/letter/n/ad7895/32.png) [@nhuurre](https://discourse.mc-stan.org/u/nhuurre)\
**Post date:** [April 21, 2022, 1:58pm UTC](https://discourse.mc-stan.org/t/rejecting-initial-value/27191/6 "2022-04-21T13:58:50Z")

</div>

`tBB` underflowing is not a problem in itself. Zero is a perfectly acceptable number. The problem comes on line 45

```stan
    real t4 = tBB(r,tb0,rc)/(t - tBB(r,tb0,rc));

```

At the initial point with `t = 0` and this amounts to undefined division 0/0.  
If I change the initial condition for the integral to avoid this problematic point

```stan
  real ComovilCoordinates[N,2]
    = integrate_ode_rk45(GeodesicEquation, {t0,0.0000001}, 0, zobs, {tb0,rc},
                         rep_array(0.0, 0), rep_array(0, 0),
                         1e-5, 1e-3, 5e2);

```

then the model initializes, at least with the small dataset I made up.