Functions block

Laxmi_Pd_Sapkota · June 15, 2022, 9:40am

Dear all, I am just started to use Stan and I am confusing whether we have to enter log of PDF function or the sum of the log of the PDF function in the functions block using the real _lpdf() signature. I have found that some authors have entered _lpdf with sum function in the functions block.(Extreme value analysis and user defined probability functions in Stan). I have shared my stan code below please suggest to me whether my code is correct or not.


functions{
  real OLGE_lpdf(real y, real a, real b, real delta, real theta){
    real z;
    real w;
    z = (1 - exp(-b*y))^a;
    w = (log(a)+log(b)+ theta*log(delta)+log(theta)- b*(y) + (a-1)*(log1m(exp(-b*y)))-
  2*log1m(z) - (theta+1)*log(delta + (z/(1-z)))); //this is my pdf just taking log
  return w;
}
  real OLGE_rng(real a, real b, real delta, real theta){
    real u;
    real v;
    u = uniform_rng(0, 1);
    v = (((1-u)*delta^(-theta))^(-1/theta) - delta);
    return (-1/b)*log1m((v/(1+v))^(1/a));
  }
}


data{
  int N;
  real y[N];
}

parameters{
  real <lower=0> a;
  real <lower=0> b;
  real <lower=0> delta;
  real <lower=0> theta;
}


model{
  for(i in 1 : N){
    y[i]~ OLGE(a, b, delta, theta);
  }
  a~ gamma(0.1, 0.1);
  b~ gamma(0.01, 0.01);
  delta~ gamma(0.2, 0.2);
  theta~ gamma(0.1, 0.1);
}


generated quantities{
  vector [N] yrep;
  for(i in 1 : N){
    yrep[i]= OLGE_rng(a, b, delta, theta);
  }
}

jsocolar · June 15, 2022, 1:12pm

A lpdf function needs to return a single real that differs by a constant from the log of the PDF; that is, it needs to compute the correct log density for its inputs. A common use case that results in the function ending with a sum is when the lpdf is vectorized over inputs corresponding to conditionally independent datapoints. In this case it is appropriate to evaluate the pointwise lpdfs and sum the result.

If the function, like yours, takes inputs of length 1, then in general no sum will be necessary or meaningful (but make sure that the function, if applied multiple times or in a loop, is applied to conditionally independent data!).

If the function takes as input vectors or arrays corresponding to data that are not conditionally independent, then it will need to compute the correct lpdf, which might include a sum, but will not correspond to the sum of some pointwise vectorized computation.

Laxmi_Pd_Sapkota · June 15, 2022, 3:24pm

Thank you so much @jsocolar, you mean my stan codes are technically correct or there is another easy way to write up stan codes, please share if any examples like my problem.

jsocolar · June 15, 2022, 4:16pm

I’m not sure I understand what your question is.

Are you wondering whether your function definition is syntactically valid (as long as it compiles then yes)?
Are you wondering whether you are computing the lpdf correctly (I don’t know, because I don’t have any reference to the function that you are trying to compute other than the stan code that you’ve provided)?
Are you wondering about the syntax for how to use this _lpdf function in the model block?

Laxmi_Pd_Sapkota · June 16, 2022, 7:07am

my question is, the functions, data, parameters, model, and generated quantity blocks are syntactically valid or not? if not or is there another easy way to code up in all blocks please suggest me.
Thanks

jsocolar · June 16, 2022, 1:16pm

The easiest way to check is to compile the model.

If the model doesn’t compile, the errors from the compiler will likely be informative. Post them here if you encounter them and aren’t sure what to do.

Laxmi_Pd_Sapkota · June 16, 2022, 3:17pm

My data set is

y<-c(3, 3, 8, 16, 15, 24, 26, 26, 38, 43, 46, 45, 58, 64, 66, 73, 
     73, 86, 89, 97, 108, 97, 146, 122, 143, 143, 106, 98, 136, 117, 
     121, 113, 100, 158, 81, 64, 37, 58, 65, 54, 73, 67, 85, 83, 102,
     107, 105, 228, 198, 271, 332, 353, 447, 405, 687, 642, 817, 
     972, 1079, 1356, 1625, 1629, 1873, 2381, 2388, 2791, 3271)

Using rstan package with OS windows 10,

fit<- stan("OLGE_Stan.stan", iter = 400, chains = 1,  
          data = list(N=N, y=y))

I have got this type of errors, I have checked thelpdf function many times but couldn’t draw samples

 SAMPLING FOR MODEL 'OLGE_Stan' NOW (CHAIN 1).
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: 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: 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: 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: 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: Rejecting initial value
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.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
[1] "error occurred during calling the sampler; sampling not done"

andrjohns · June 16, 2022, 3:52pm

The error is caused by these sections of your likelihood:

z = (1 - exp(-b*y))^a;
...log1m(z)...
...z/(1-z)...

When the exp function is passed values that are < -709, the result is so small that it is treated as zero. This means that calculation for z then becomes: (1 - 0)^a, which is just 1.

The call to log1m(z) then evaluates to log(0) and the call to z/(1-z) then evaluates to 1/0

Laxmi_Pd_Sapkota · June 16, 2022, 4:34pm

Thank you so much for your suggestion, but how to fix these issues, please suggest to me it will be a great help for this type of work in future,

jsocolar · June 16, 2022, 5:07pm

Do everything on the scale of log(z) using numerically stable functions.

log_z = a * log1m_exp(-b*y);
... log1m_exp(log_z) ...
... exp(log_z - log1m_exp(log_z)) ...

I can’t guarantee that that’ll be numerically stable enough (the last line looks ominous to me) but it’s a step in the right direction I think. Hopefully you don’t need to exponentiate the third line above and can continue to work through numerically stable log-scale functions.

jsocolar · June 16, 2022, 5:17pm

Ah, yes, since delta is constrained to be positive, you can replace log(delta + (z/(1-z))) with log(delta) + log1p_exp(log_z - log1m_exp(log_z) - log(delta)). Double check my math! But something like that anyway.

Laxmi_Pd_Sapkota · June 17, 2022, 2:37am

My PDF is like this

f\left( {x;\Phi } \right) = \frac{{ab\delta ^\theta \theta e^{ - bx} \left( {1 - e^{ - bx} } \right)^{a - 1} }} {{\left[ {1 - \left( {1 - e^{ - bx} } \right)^a } \right]^2 }}\left[ {\delta + \frac{{\left( {1 - e^{ - bx} } \right)^a }} {{1 - \left( {1 - e^{ - bx} } \right)^a }}} \right]^{ - \theta - 1} ;x > 0

where \Phi = \left( {a,b,\delta ,\theta } \right)>0 is a vector of parameter space
and I have supposed z = (1 - exp(-b*y))^a after taking log f(x) I think we could not define log of z like log_z = a * log1m_exp(-b*y).in local variable. Please have a look my PDF and suggest me how to define it in function block.

As you suggest

since delta is constrained to be positive, you can replace log(delta + (z/(1-z))) with log(delta) + log1p_exp(log_z - log1m_exp(log_z) - log(delta)) . I have run by correcting this line only but out put is same as earlier that I have post ie could not generate samples.
thank you

jsocolar · June 17, 2022, 7:10am

I don’t understand what you mean when you say

This relationship follows directly from

Given your PDF, I would take the logarithm of it and then try to simplify in such a way that we never use exponentiation directly, and instead use only things like log1m_exp, log1p_exp, log_sum_exp etc that are specifically written to be numerically stable. For example, after taking the logarithm the numerator of the first term becomes log(a) + log(b) + theta * log(delta) + log(theta) - b*x + (a - 1) * log1m_exp(-b*x). The denominator of the first term becomes -2 * log1m_exp(a * log1m_exp(-b * x)). And so forth.

Laxmi_Pd_Sapkota · June 17, 2022, 4:19pm

You are absolutely right i have written PDF as

functions{
  real OLGE_lpdf(real y, real a, real b, real delta, real theta){
    real log_z;
    real w;
    log_z = a*log1m_exp(-b*y);
    w = log(a)+log(b)+ theta*log(delta)+log(theta)- b*(y) + (a-1)*(log1m(exp(-b*y)))-
  2*log1m_exp(a*log1m_exp(-b*y)) - (theta+1)*(log(delta) + log1p_exp(log_z - log1m_exp(log_z)));
  return w;
}

fit<- stan("OLGE_Stan.stan", iter = 4000, chains = 1,  
           data = list(N=N, y=y))

and got the output like this

Inference for Stan model: OLGE_Stan.
1 chains, each with iter=4000; warmup=2000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=2000.

             mean se_mean   sd     2.5%      25%      50%      75%
a            0.61    0.00 0.06     0.49     0.57     0.61     0.66
b            0.00    0.00 0.00     0.00     0.00     0.00     0.00
delta        0.00     NaN 0.00     0.00     0.00     0.00     0.00
theta       10.25    0.68 7.41     1.84     4.80     8.28    13.76
yrep[1]      0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[2]      0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[3]      0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[4]      0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[5]      0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[6]      0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[7]      0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[8]      0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[9]      0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[10]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[11]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[12]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[13]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[14]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[15]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[16]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[17]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[18]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[19]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[20]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[21]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[22]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[23]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[24]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[25]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[26]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[27]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[28]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[29]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[30]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[31]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[32]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[33]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[34]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[35]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[36]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[37]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[38]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[39]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[40]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[41]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[42]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[43]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[44]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[45]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[46]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[47]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[48]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[49]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[50]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[51]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[52]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[53]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[54]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[55]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[56]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[57]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[58]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[59]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[60]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[61]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[62]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[63]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[64]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[65]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[66]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
yrep[67]     0.00     NaN 0.00     0.00     0.00     0.00     0.00
lp__     46672.08    0.14 1.66 46667.81 46671.34 46672.48 46673.22
            97.5% n_eff Rhat
a            0.74   817 1.00
b            0.00   129 1.04
delta        0.00   NaN  NaN
theta       29.71   117 1.03
yrep[1]      0.00   NaN  NaN
yrep[2]      0.00   NaN  NaN
yrep[3]      0.00   NaN  NaN
yrep[4]      0.00   NaN  NaN
yrep[5]      0.00   NaN  NaN
yrep[6]      0.00   NaN  NaN
yrep[7]      0.00   NaN  NaN
yrep[8]      0.00   NaN  NaN
yrep[9]      0.00   NaN  NaN
yrep[10]     0.00   NaN  NaN
yrep[11]     0.00   NaN  NaN
yrep[12]     0.00   NaN  NaN
yrep[13]     0.00   NaN  NaN
yrep[14]     0.00   NaN  NaN
yrep[15]     0.00   NaN  NaN
yrep[16]     0.00   NaN  NaN
yrep[17]     0.00   NaN  NaN
yrep[18]     0.00   NaN  NaN
yrep[19]     0.00   NaN  NaN
yrep[20]     0.00   NaN  NaN
yrep[21]     0.00   NaN  NaN
yrep[22]     0.00   NaN  NaN
yrep[23]     0.00   NaN  NaN
yrep[24]     0.00   NaN  NaN
yrep[25]     0.00   NaN  NaN
yrep[26]     0.00   NaN  NaN
yrep[27]     0.00   NaN  NaN
yrep[28]     0.00   NaN  NaN
yrep[29]     0.00   NaN  NaN
yrep[30]     0.00   NaN  NaN
yrep[31]     0.00   NaN  NaN
yrep[32]     0.00   NaN  NaN
yrep[33]     0.00   NaN  NaN
yrep[34]     0.00   NaN  NaN
yrep[35]     0.00   NaN  NaN
yrep[36]     0.00   NaN  NaN
yrep[37]     0.00   NaN  NaN
yrep[38]     0.00   NaN  NaN
yrep[39]     0.00   NaN  NaN
yrep[40]     0.00   NaN  NaN
yrep[41]     0.00   NaN  NaN
yrep[42]     0.00   NaN  NaN
yrep[43]     0.00   NaN  NaN
yrep[44]     0.00   NaN  NaN
yrep[45]     0.00   NaN  NaN
yrep[46]     0.00   NaN  NaN
yrep[47]     0.00   NaN  NaN
yrep[48]     0.00   NaN  NaN
yrep[49]     0.00   NaN  NaN
yrep[50]     0.00   NaN  NaN
yrep[51]     0.00   NaN  NaN
yrep[52]     0.00   NaN  NaN
yrep[53]     0.00   NaN  NaN
yrep[54]     0.00   NaN  NaN
yrep[55]     0.00   NaN  NaN
yrep[56]     0.00   NaN  NaN
yrep[57]     0.00   NaN  NaN
yrep[58]     0.00   NaN  NaN
yrep[59]     0.00   NaN  NaN
yrep[60]     0.00   NaN  NaN
yrep[61]     0.00   NaN  NaN
yrep[62]     0.00   NaN  NaN
yrep[63]     0.00   NaN  NaN
yrep[64]     0.00   NaN  NaN
yrep[65]     0.00   NaN  NaN
yrep[66]     0.00   NaN  NaN
yrep[67]     0.00   NaN  NaN
lp__     46674.11   137 1.00

Still there is problem, what is happen

jsocolar · June 17, 2022, 4:38pm

It looks like the model is working fine, but there’s still an underflow problem in the RNG and therefore in the generated quantities.

Laxmi_Pd_Sapkota · June 18, 2022, 3:38pm

The RNG that I have written as

real OLGE_rng(real a, real b, real delta, real theta){
    real u;
    real log_v;
    u = uniform_rng(0, 1);
    log_v = log(-delta)+ log1p(exp((-1/theta)*(log1m(u)-
    theta*log(delta)))- log(-delta));
    return (-1/b)*(log1m_exp((1/a)*(log_v - log1p_exp(log_v))));

and my actual Quantile function is
Q\left( p \right) = - \frac{1} {b}\left[ {\log \left\{ {1 - \left( {\frac{v} {{1 +v}}} \right)^{1/a} } \right\}} \right]

where v = \left[ {\left\{ {\left( {1 - p} \right)\delta ^{ - \theta } } \right\}^{ - 1/\theta } - \delta } \right]
is my RNG function is correct according to this?

after compiling with this RNG function I got the same result as earlier

1 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=1000.

             mean se_mean   sd     2.5%      25%      50%      75%
a            0.62    0.01 0.07     0.50     0.57     0.62     0.66
b            0.00    0.00 0.00     0.00     0.00     0.00     0.00
delta        0.00     NaN 0.00     0.00     0.00     0.00     0.00
theta        9.70    1.50 6.66     2.28     4.93     7.83    12.52
yrep[1]       NaN      NA   NA       NA       NA       NA       NA
yrep[2]       NaN      NA   NA       NA       NA       NA       NA
yrep[3]       NaN      NA   NA       NA       NA       NA       NA
yrep[4]       NaN      NA   NA       NA       NA       NA       NA
yrep[5]       NaN      NA   NA       NA       NA       NA       NA
yrep[6]       NaN      NA   NA       NA       NA       NA       NA
yrep[7]       NaN      NA   NA       NA       NA       NA       NA
yrep[8]       NaN      NA   NA       NA       NA       NA       NA
yrep[9]       NaN      NA   NA       NA       NA       NA       NA
yrep[10]      NaN      NA   NA       NA       NA       NA       NA
yrep[11]      NaN      NA   NA       NA       NA       NA       NA
yrep[12]      NaN      NA   NA       NA       NA       NA       NA
yrep[13]      NaN      NA   NA       NA       NA       NA       NA
yrep[14]      NaN      NA   NA       NA       NA       NA       NA

andrjohns · June 18, 2022, 4:35pm

Your RNG is attempting to take the log of a negative number:

log(-delta)

Which returns NaN and causes the final result to also be NaN

Laxmi_Pd_Sapkota · June 19, 2022, 2:33am

I have just correct the Rng function as suggested by @andrjohns
log(-delta) returns NAN.

real OLGE_rng(real a, real b, real delta, real theta){
    real u;
    real log_v;
    u = uniform_rng(0, 1);
    log_v = log_diff_exp((-1/theta)*(log1m(u)-
    theta*log(delta)), log(delta));
    return (-1/b)*(log1m_exp((1/a)*(log_v - log1p_exp(log_v))));

but still get that same output as earlier

andrjohns · June 19, 2022, 12:17pm

You’ll need to do some debugging and testing to see where in your RNG the error is occurring. Try breaking the calculation into separate parts and seeing which part is returning NaN. You can also add print() statements to have Stan print out the values of particular variables/computations during sampling

Topic		Replies	Views
Outputting lpdf functions as vectors instead of reals Developers	6	915	January 3, 2022
Different types of reals in documentation and user defined functions? Modeling	14	1062	November 2, 2017
Custom pdf with 2 variables Modeling	2	577	April 7, 2021
Lambert transformation in user defined function to generate random numbers General rstan	1	348	January 26, 2023
Integrate_1d: fifth argument must be data only Modeling	3	361	July 2, 2021

Functions block

Related Topics