# Initial value error (log(0)) for user defined beta prior distribution in hierarchical generalized extreme value (GEV) distribution

**URL:** <https://discourse.mc-stan.org/t/initial-value-error-log-0-for-user-defined-beta-prior-distribution-in-hierarchical-generalized-extreme-value-gev-distribution/29967>\
**Category:** Modeling\
**Tags:** initialization, rstan\
**Created:** [January 8, 2023, 7:49pm UTC](https://discourse.mc-stan.org/t/initial-value-error-log-0-for-user-defined-beta-prior-distribution-in-hierarchical-generalized-extreme-value-gev-distribution/29967 "2023-01-08T19:49:14Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Alex\_Lee](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/alex_lee/32/13479_2.png) [@Alex\_Lee](https://discourse.mc-stan.org/u/Alex_Lee)\
**Post date:** [January 8, 2023, 7:49pm UTC](https://discourse.mc-stan.org/t/initial-value-error-log-0-for-user-defined-beta-prior-distribution-in-hierarchical-generalized-extreme-value-gev-distribution/29967/1 "2023-01-08T19:49:14Z")

</div>

I have a user defined Beta distribution.

Beta(\theta+0.5|p,q) = (0.5 + \theta)^{p-1} (0.5 - \theta)^{q-1} / B(p,q) between [-0.5,0.5]

where B(p,q) = \Gamma(p) \Gamma(q) / \Gamma(p+q)

```plaintext
real mbeta_lpdf(real theta, real p, real q) {
      real lpp;
      real b;
      lpp = pow(0.5 + theta,p-1) * pow(0.5 - theta,q-1);
      b = (tgamma(p-1)*tgamma(q-1))/tgamma(p+q-1);
      return log(lpp) - log(b);

```

This Beta distribution is implemented into a hierarchical generalized extreme value (GEV) distribution for the shape parameter (`xi`). For the location (`mu`) and scale (`sigma`) parameters I use a normal and gamma prior distributions respectively with uniform hyperpriors (`alpha = uniform(-infinity,+infinity)` and `delta = uniform(0,+infinity)`.

This is the whole model:

```plaintext
functions {
    real gev_lpdf(vector y, real mu, real sigma, real xi) {
        int N = rows(y);
        
        vector[rows(y)] z;
        vector[rows(y)] lp;
        
        
        for(n in 1:N) {
            z[n] = 1 + (y[n] - mu) * xi / sigma;
            lp[n] = (1 + (1 / xi)) * log(z[n]) + pow(z[n], -1/xi);
        }        
        return -sum(lp) - N*log(sigma);
    }
    
    real mbeta_lpdf(real theta, real p, real q) {
      real lpp;
      real b;
      lpp = pow(0.5 + theta,p-1) * pow(0.5 - theta,q-1);
      b = (tgamma(p-1)*tgamma(q-1))/tgamma(p+q-1);
      return log(lpp) - log(b);
    }
}

data {
  int<lower=0> N; //observation and duration vector length
  int<lower=1> J; //number of weather stations
  vector[N] y[J]; //data matrix
}

parameters {
real<lower=-0.5,upper=0.5> xi[J];
real mu[J];
real<lower=0> sigma[J];
real alpha;
real<lower=0> delta;
}

model {
  for(j in 1:J){
    xi[j] ~ mbeta(6.0,9.0);
    mu[j] ~ normal(alpha,25);
    sigma[j] ~ gamma(delta,delta);
    y[j] ~ gev(mu[j],sigma[j],xi[j]);
    }
}

```

When I run this code unfortunately it throws an error:

```plaintext
...
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

```

Where have I gone wrong here?

---

<div class="post-metadata">

**Author:** ![hsusmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hsusmann/32/11775_2.png) [@hsusmann](https://discourse.mc-stan.org/u/hsusmann)\
**Post date:** [January 8, 2023, 8:43pm UTC](https://discourse.mc-stan.org/t/initial-value-error-log-0-for-user-defined-beta-prior-distribution-in-hierarchical-generalized-extreme-value-gev-distribution/29967/2 "2023-01-08T20:43:23Z")

</div>

Are you able to post some example code that runs the model? That would help with debugging. Could you also post the mathematical expression you’re using for the generalized GEV distribution?

At first glance, I’m noticing the `1/xi` terms in `gev_lpdf`, which will blow up when `xi` is close to zero (it looks like `xi` can be from `-0.5` to `0.5`). I’m not sure if that’s the issue, but it should probably be addressed.

---

<div class="post-metadata">

**Author:** ![Alex\_Lee](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/alex_lee/32/13479_2.png) [@Alex\_Lee](https://discourse.mc-stan.org/u/Alex_Lee)\
**Post date:** [January 8, 2023, 9:45pm UTC](https://discourse.mc-stan.org/t/initial-value-error-log-0-for-user-defined-beta-prior-distribution-in-hierarchical-generalized-extreme-value-gev-distribution/29967/3 "2023-01-08T21:45:43Z")

</div>

Thanks for the reply.

Yes xi should not be equal to zero otherwise the GEV becomes the Gumbel distribution.

Maybe I should constrain xi to be `!= 0`, but where and how to do that I don’t know (I’m fairly new to Stan).

Here is my R code:

```plaintext
library(rstan)

vec = readRDS("..../file.rds")

col_len = length(vector[,,1])

y = matrix(NA,nrow = 20, ncol = 41)

for (i in 1:20) {
  y[i,] = IDFlist[,1,i]
}

N = 41
J = 20

dir = ".../stackoverflow/"
model = stan_model(paste0(dir,"GEV_hier.stan"))

fit = sampling(model,list(N=N,D=D,J=J,y=y,d=d),iter=500,chains=1,seed=10)

```

---

<div class="post-metadata">

**Author:** ![hsusmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hsusmann/32/11775_2.png) [@hsusmann](https://discourse.mc-stan.org/u/hsusmann)\
**Post date:** [January 8, 2023, 11:37pm UTC](https://discourse.mc-stan.org/t/initial-value-error-log-0-for-user-defined-beta-prior-distribution-in-hierarchical-generalized-extreme-value-gev-distribution/29967/4 "2023-01-08T23:37:41Z")

</div>

Thanks! I can’t run that code though because it depends on some data files that I don’t have. Maybe you could simulate some sample data as an example? (Fitting to simulated data is a good thing to do anyway, as you can see if the estimated parameters match the simulation parameters.)

As good general debugging strategy when you have an error message but you’re not sure what exactly is causing it is to start removing model components until you get to a simpler model that works. Then, build slowly back up to the original model. For example, you could remove the `mbeta` prior on `xi`. That is, comment out

```plaintext
xi[j] ~ mbeta(6.0,9.0);

```

Then run the code and see if you still have an error. If you do, keep removing model components – maybe assign a simpler response distribution to `y` instead of `gev`, for example. If we can localize the problem it will be much easier to troubleshoot.

---

<div class="post-metadata">

**Author:** ![Alex\_Lee](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/alex_lee/32/13479_2.png) [@Alex\_Lee](https://discourse.mc-stan.org/u/Alex_Lee)\
**Post date:** [January 9, 2023, 8:22am UTC](https://discourse.mc-stan.org/t/initial-value-error-log-0-for-user-defined-beta-prior-distribution-in-hierarchical-generalized-extreme-value-gev-distribution/29967/5 "2023-01-09T08:22:03Z")

</div>

Thanks for the advice.

I’ve rewritten the code and added simulated GEV samples.

Please see below:

```plaintext
library(rstan)

y = matrix(NA,nrow = 20, ncol = 41)

mu = runif(20,15,20)
sigma = runif(20,3,5)
xi = runif(20,-0.5,0.5)

library(extRemes)

for (i in 1:20) {
  y[i,] = revd(41,loc=mu[i],scale=sigma[i],shape=xi[i],type="GEV")
}

N = 41
J = 20

dir = ".../stackoverflow/"
model = stan_model(paste0(dir,"GEV_hier.stan"))

fit = sampling(model,list(N=N,D=D,J=J,y=y),iter=500,chains=1,seed=10)

```

I was able to reproduce the error with this code.

---

<div class="post-metadata">

**Author:** ![hsusmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hsusmann/32/11775_2.png) [@hsusmann](https://discourse.mc-stan.org/u/hsusmann)\
**Post date:** [January 10, 2023, 12:04am UTC](https://discourse.mc-stan.org/t/initial-value-error-log-0-for-user-defined-beta-prior-distribution-in-hierarchical-generalized-extreme-value-gev-distribution/29967/6 "2023-01-10T00:04:22Z")

</div>

Thanks for the example code.

It looks like the issue is that `z[n]` can be negative depending on the value of `mu`, `sigma`, and `xi`, which then causes `log(z[n])` to be undefined in the `gev_lpdf` calculation.

I think this comes down to how the support of the GEV distribution depends on the values of `mu`, `sigma`, and `xi`. (See “Support” in the sidebar here: [Generalized extreme value distribution - Wikipedia](https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution)).

One way this could be fixed by constraining the values of `mu`, `sigma,` and `xi` so that the observations are always within the support of the GEV. You might look to this post: [Bayesian data analysis of extreme values with Stan and rstan. – Georgios Boumis](https://gboumis.com/2022/10/09/750/), which appears to take the approach of constraining `xi`.

---

<div class="post-metadata">

**Author:** ![Alex\_Lee](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/alex_lee/32/13479_2.png) [@Alex\_Lee](https://discourse.mc-stan.org/u/Alex_Lee)\
**Post date:** [January 10, 2023, 1:06pm UTC](https://discourse.mc-stan.org/t/initial-value-error-log-0-for-user-defined-beta-prior-distribution-in-hierarchical-generalized-extreme-value-gev-distribution/29967/7 "2023-01-10T13:06:37Z")

</div>

Excellent! That solved the problem.