Getting "Error in function stan::prob::gamma_log(N4stan5agrad3varE)"

[rubi]

Hi (again),

I have a the following hierarchical model with parameter transformations:

At the bottom of the hierarchy I have:

O[g][s][k] ~ Beta(a,b) - that is, Observation k from subgroup s of group g follows a beta distribution.

I reparametrize the a and b parameters in terms of the mean (mu) and variance (sigmaS) of the beta distribution:

a = -mu(sigmaS + mu^2 -mu)/sigmaS

b = (sigmaS + mu^2 -mu)(mu - 1)/sigmaS

The next level up in hierarchy thus has priors for mu and sigmaS:

mu ~ Beta(1,1)

sigmaS ~ Gamma(c,d)

The next level up in hierarchy (the top one) has priors for c and d:

c ~ Exponential(1)

d ~ Exponential(1)

This my code:

data {

  int<lower=1> G; // number of groups

  int<lower=1> S; // number of subgroups

  int<lower=1> K; // number of observations in each subgroup

  matrix<lower=0,upper=1>[S,K] O[G];

}

parameters {

  real mu; 

  real sigmaS;

  real c;

  real d;

}

transformed parameters {

  real a; 

  real b;

  a <- -1*mu*(sigmaS + mu*mu - mu)/sigmaS;

  b <- (sigmaS + mu*mu - mu)*(mu - 1)/sigmaS; 

}

model {

   for (g in 1:G){

c ~ exponential(1);

d ~ exponential(1);

for (s in 1:S){

mu ~ beta(1,1);

sigmaS ~ gamma(c,d);

for (k in 1:K){

O[g,s,k] ~ beta(a,b);

}

}

  }

}

I'm getting this error:

Error : Error in function stan::prob::gamma_log(N4stan5agrad3varE): Shape parameter is -1.2393059165125289:0, but must be > 0!

error occurred during calling the sampler; sampling not done

I also tried constraining a and b to <lower=0> and mu to <lower=0,upper=1> (since it is the mean of O[g][s][k] for each [g][s] vector and O[g][s][k]'s are constrained <lower=0,upper=1>) but the problem persists.

Any help would be greatly appreciated.

Thanks a lot,

rubi

[Marcus Brubaker]

a, b, c and d must be constrained to be positive.  I.E., you also need:

real<lower=0> c;

real<lower=0> d;

Cheers,

Marcus

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[rubi]

Thanks for the help.

I set:

parameters {

  real<lower=0,upper=1> mu; 

  real<lower=0> sigmaS;

  real<lower=0> c;

  real<lower=0> d;

}

transformed parameters {

  real<lower=0> a; 

  real<lower=0> b;

  a <- -1*mu*(sigmaS + mu*mu - mu)/sigmaS;

  b <- (sigmaS + mu*mu - mu)*(mu - 1)/sigmaS; 

}

 However now I'm getting this error:

Error : Invalid value of a: Error in function validate transformed params N4stan5agrad3varE: a is -0.4640491230132906:0, but must be greater than 0

I'm puzzled as to how negative parameter values are sampled where all my limits are set to positive and all the distributions only have positive supports.

Any idea?

Thanks again,

rubi

[Bob Carpenter]

On 2/21/13 11:01 PM, rubi wrote:
> Thanks for the help.
> I set:
> parameters {
> real<lower=0,upper=1> mu;
> real<lower=0> sigmaS;
> real<lower=0> c;
> real<lower=0> d;
> }
> transformed parameters {
> real<lower=0> a;
> real<lower=0> b;
> a <- -1*mu*(sigmaS + mu*mu - mu)/sigmaS;
> b <- (sigmaS + mu*mu - mu)*(mu - 1)/sigmaS;
> }
> However now I'm getting this error:
> Error : Invalid value of a: Error in function validate transformed params N4stan5agrad3varE: a is -0.4640491230132906:0,
> but must be greater than 0
>
> I'm puzzled as to how negative parameter values are sampled where all my limits are set to positive and all the
> distributions only have positive supports.
>
> Any idea?
>
> Thanks again,
> rubi

Stan requires all transformed parameters to satisfy
their declared constraints at the end of the transformed
parameter statements.

Your transformed parameter a is declared to be positive.
But its definition will produce a negative result whenever

sigmaS > mu - mu*mu.

If you need to keep parameter a positive, then declare sigmaS to
satisfy the necessary constraint with

real<lower=0, upper=(mu - mu * mu)> sigmaS;

Stan's constraint language is limited by the fact that the
constraints have to be declared in order. But that'll work
here.

Better yet, find an alternative simpler way to parameterize
the model.

- Bob

[Ben Goodrich]

On Thursday, February 21, 2013 6:29:52 PM UTC-5, rubi wrote:

Hi (again),

I have a the following hierarchical model with parameter transformations:

At the bottom of the hierarchy I have:

O[g][s][k] ~ Beta(a,b) - that is, Observation k from subgroup s of group g follows a beta distribution.

I reparametrize the a and b parameters in terms of the mean (mu) and variance (sigmaS) of the beta distribution:

a = -mu(sigmaS + mu^2 -mu)/sigmaS

b = (sigmaS + mu^2 -mu)(mu - 1)/sigmaS

The next level up in hierarchy thus has priors for mu and sigmaS:

mu ~ Beta(1,1)

sigmaS ~ Gamma(c,d)

What Bob said will fix your problems with the parameter support for sigmaS. Another thing is that if you are going to use a uniform prior on mu, that is fine but you should just comment out the

mu ~ beta(1,1);

line. It won't make a huge difference in terms of speed, but there is some overhead in evaluating that constant density, which doesn't affect the shape of the posterior distribution. In the absence of an explicitly specified prior, parameters are uniformally distributed over their support.

Ben

[rubi]

Thanks a lot for the help.

I followed your advice and set the constraint on sigmaS and deleted the sampling from beta(1,1). I also changed sigmaS to be sampled from beta instead of from gamma, because sigmaS is the variance of beta distributed observations. (See code below)

As Bob's response alluded to perhaps the reparameterization is not ideal, because now I'm getting this error:

Initialization failed after 100 attempts.  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

The thing is that my observations are rates (i.e., 0 <= O[g][s][k] <= 1 for every k,s,g) so sampling them from beta makes sense, and my aim is to model their variance across subgroups and groups as I believe it is shared across these levels of the hierarchy. That's why my solution was to reparametrize the the variance of the beta distribution as a function of a and b and put a prior on that variance. If you easily see a straight forward alternative solution I'd be really grateful to hear about it.

Thanks again,

rubi

Code

data {

  int<lower=1> G; // number of groups

  int<lower=1> S; // number of subgroups

  int<lower=1> K; // number of observations in each subgroup

  matrix<lower=0,upper=1>[S,K] O[G];

}

parameters {

  real<lower=0,upper=1> mu; 

  real<lower=0,upper=(mu - mu * mu)> sigmaS;

  real<lower=0> c;

  real<lower=0> d;

}

transformed parameters {

  real<lower=0> a; 

  real<lower=0> b;

  a <- -1*mu*(sigmaS + mu*mu - mu)/sigmaS;

  b <- (sigmaS + mu*mu - mu)*(mu - 1)/sigmaS; 

}

model {

   for (g in 1:G){

c ~ exponential(1);

d ~ exponential(1);

for (s in 1:S){

sigmaS ~ beta(c,d);

for (k in 1:K){

O[g,s,k] ~ beta(a,b);

}

}

  }

}

[Ben Goodrich]

On Friday, February 22, 2013 4:43:33 PM UTC-5, rubi wrote:

I followed your advice and set the constraint on sigmaS and deleted the sampling from beta(1,1). I also changed sigmaS to be sampled from beta instead of from gamma, because sigmaS is the variance of beta distributed observations. (See code below)

Since sigmaS is now bounded from above by mu * (1 - mu) < 1, a beta distribution for sigmaS is problematic because it puts positive mass in the region from mu * (1 - mu) to 1. Since the amount of mass is not a constant but a function of the unknown mu, this mismatch affects the log-posterior and the partial derivatives thereof. This may be related to the error you are getting.

Thus, at a minimum, you should truncate your beta prior to the [0,mu * (1 - mu)] interval, which may entail a noticeable slowdown because the beta CDF is somewhat expensive to evaluate. Alternatively, you could just omit the explicit prior on sigmaS if an implicit uniform prior between 0 and mu * (1 - mu) is acceptable to you.

Ben

[rubi]

Thanks a lot Ben,

I assumed what you mean in the code is something like:

sigmaS ~ beta(c,d) T[,mu - mu * mu];

I couldn't understand from the examples in the manual if this syntax is legitimate or not, and neither from the compilation error message:

no matches for function name="beta_cdf"

    arg 0 type=real

    arg 1 type=real

    arg 2 type=real

upper truncation not defined for specified arguments to beta

So is my syntax wrong or is a truncated beta not possible?

Thanks

[Ben Goodrich]

On Friday, February 22, 2013 6:34:17 PM UTC-5, rubi wrote:

Thanks a lot Ben,

I assumed what you mean in the code is something like:

sigmaS ~ beta(c,d) T[,mu - mu * mu];

I couldn't understand from the examples in the manual if this syntax is legitimate or not, and neither from the compilation error message:

no matches for function name="beta_cdf"

    arg 0 type=real

    arg 1 type=real

    arg 2 type=real

upper truncation not defined for specified arguments to beta

So is my syntax wrong or is a truncated beta not possible?

I guess we hadn't implemented the incomplete beta function to do the CDF as of v1.1.1. It is implemented on GitHub now. In the meantime, I put an example of implementing the Kumaraswamy distribution at

https://groups.google.com/d/msg/stan-users/4gv3fNCqSNk/4dwmJj1YwbsJ

It has an easy CDF too. You would just need to do

lp__ <- lp__ - log1m(CDF);

where CDF is the Kumaraswamy CDF evaluated at mu * (1 - mu).

Ben

[Daniel Lee]

All the univariate cdfs will be available in the next version, v1.2.0. The cdfs aren't available on the master branch of github yet, but will be soon.

Daniel

[rubi]

Thanks again for help.

I hope I understood Ben's suggestion correctly. My model is now:

model {

   for (g in 1:G){

c ~ exponential(1);

d ~ exponential(1);

for (s in 1:S){

lp__ <- lp__ - log1m(1-pow((1-pow(mu - mu * mu,c)),d)); //where 1-pow((1-pow(mu - mu * mu,c)),d) is the Kumaraswamy(c,d) CDF evaluated at mu  - mu*mu, which is supposed to correct lp for no mass in the interval [mu - mu * mu,1]

for (k in 1:K){

O[g,s,k] ~ beta(a,b);

}

}

  }

}

I'm still getting this error though:

Initialization failed after 100 attempts.  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

error occurred during calling the sampler; sampling not done

But I think that this is due to the fact the in O[g,s,k] ~ beta(a,b), the beta(a,b) is truncated (as beta(c,d) is for sigmaS ~ beta(c,d)).

Would it make sense to perform the same correction for O[g,s,k] ~ beta(a,b) (i.e., lp__ <- lp__ - log1m(CDF); where in this case CDF is the Kumaraswamy CDF evaluated at at some point which I should be close to the truncation point (say mu+3*sqrt(sigmaS))?

Thanks again,

rubi

[Bob Carpenter]

Can you send the full model, and not just the model
block?

You almost certainly don't want to evaluate
the sampling statement for c and d a total
of G times:

> for (g in 1:G){
> c ~ exponential(1);
> d ~ exponential(1);

That's equivalent to saying

p(c) propto pow(Exponential(c|1), G)

What is c?

That error arises when the sampler can't move at all
from 100 different tries at random initialization.

It can sometimes come from values being too extreme
which can cause arithmetic issues.

- Bob

> stan-users+...@googlegroups.com <javascript:>.
> For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

[rubi]

Sure.

Attached is the model and a toy data set. 

I create a list (called O) of size G (# of groups) of SxK matrices (S is # subgroups and K is the # of observations in each subgroup) and call rstan with:

stan(file=model.stan,data=dat,iter=1000,chains=4)
where dat = list(G=G,S=S,K=K,O=O)

The model is:

At the bottom of the hierarchy I have:

O[g][s][k] ~ Beta(a,b) - that is, Observation k from subgroup s of group g follows a beta distribution.

I reparametrize the a and b parameters in terms of the mean (mu) and variance (sigmaS) of the beta distribution:

a = -mu(sigmaS + mu^2 -mu)/sigmaS

b = (sigmaS + mu^2 -mu)(mu - 1)/sigmaS

The next level up in hierarchy models subgroup variance and thus has a prior sigmaS and a flat prior for mu:

mu ~ Beta(1,1)

sigmaS ~ Gamma(c,d)

The next level up in the hierarchy (the top one) models the group level with c and d:

c ~ Exponential(1)

d ~ Exponential(1)

Since I have G groups I thought I need sample c and d inside the loop over g. Am I wrong?

Thanks,

rubi

[Ben Goodrich]

On Saturday, February 23, 2013 6:10:46 PM UTC-5, rubi wrote:

Sure.

Attached is the model and a toy data set. 

Did you upload the wrong thing for data.txt? It doesn't seem to be in Stan format:

goodrich@CYBERPOWERPC:/tmp/SU/rubi$ cat data.txt
G       S       Ogs
1.2     293T    0.7,0.72,0.69,0.68,0.71
1.2     MCF7    0.5,0.48,0.55,0.51,0.52
1.2     U937    0.3,0.4,0.35,0.38,0.39
3.4     293T    0.2,0.19,0.26,0.24,0.35
3.4     MCF7    0,0.01,0.05,0.04,0.03
3.4     U937    0.8,0.85,0.78,0.81,0.82
5.6     293T    0.55,0.58,0.54,0.52,0.51
5.6     MCF7    0.42,0.45,0.44,0.46,0.43
5.6     U937    0.74,0.75,0.71,0.7,0.72
7.8     293T    0.2,0.18,0.19,0.21,0.22
7.8     MCF7    0.4,0.48,0.49,0.5,0.51
7.8     U937    0.6,0.68,0.69,0.7,0.71

Ben

[rubi]

Sorry about that.

Attached is an R wrapper script that parses data.txt and calls model.stan with it's appropriate input list.

usage: wrapper.R <data.txt> <model.stan>

Thanks

[Bob Carpenter]

You need to scope each of your variables properly
in loops.

for (g in 1:G){
c ~ exponential(1);
d ~ exponential(1);
for (s in 1:S){

//mu ~ beta(1,1);
//sigmaS ~ beta(c,d);

lp__ <- lp__ - log1m(1-pow((1-pow(mu - mu * mu,c)),d));

for (k in 1:K){
O[g,s,k] ~ beta(a,b);
}
}
}

All of the sampling statements other than the one
for O[g,s,k] should be outside the loop.

I don't understand why you're making the priors for
a and b so complicated. A more natural way to do this
with an exponential scale would be:

parameters {
real<lower=0> lambda;
real<lower=0,upper=1> theta;

real<lower=0> a;
real<lower=0> b;
}

model {
lambda ~ exponential(1); // use any prior with support on (0,inf)
theta ~ beta(1,1); // can remove because it's uniform and theta bounded

for (g in 1:G)

for (s in 1:S)

for (k in 1:K)

O[g,s,k] ~ beta(lambda*theta,lambda*(1 - theta));
}

The standard deviation isn't particularly natural for
beta because it's bounded. The variable theta will
get the hierarchical mean and lamda the hierarchical
prior count. Andrew suggested a Pareto prior on lambda
in Bayesian Data Analysis (chapter 5, which starts with
a very similar model), but I don't think he likes that choice
any more.

You can make it faster by vectorizing the last loop, but
I wouldn't worry about that until you have a version
working.

- Bob

On 2/23/13 6:10 PM, rubi wrote:
> Sure.
> Attached is the model and a toy data set.
> I create a list (called O) of size G (# of groups) of SxK matrices (S is # subgroups and K is the # of observations in
> each subgroup) and call rstan with:
> stan(file=model.stan,data=dat,iter=1000,chains=4)
> where dat = list(G=G,S=S,K=K,O=O)
> The model is:
> At the bottom of the hierarchy I have:

> O[g][s][k] ~ Beta(/a/,/b/) - that is, Observation k from subgroup s of group g follows a beta distribution.
> I reparametrize the /a/ and /b/ parameters in terms of the mean (/mu/) and variance (/sigmaS/) of the beta distribution:
> /a/ = -/mu/(/sigmaS/ + /mu/^2 -/mu/)//sigmaS/
> /b/ = (/sigmaS/ + /mu/^2 -/mu/)(/mu/ - 1)//sigmaS/
>
> The next level up in hierarchy models subgroup variance and thus has a prior /sigmaS /and a flat prior for /mu/:
> /mu/ ~ Beta(1,1)
> /sigmaS/ ~ Gamma(/c/,/d/)
>
> The next level up in the hierarchy (the top one) models the group level with /c/ and /d/:
> /c/ ~ Exponential(1)
> /d/ ~ Exponential(1)

> <https://groups.google.com/groups/opt_out> <https://groups.google.com/groups/opt_out

[Ben Goodrich]

On Sunday, February 24, 2013 2:57:47 AM UTC-5, Bob Carpenter wrote:

I don't understand why you're making the priors for 
a and b so complicated. 

There is that. But I can't get any beta model to work on these data.

goodrich@CYBERPOWERPC:/tmp/SU/rubi$ ./test --test_grad --data=dat.txt --init=0

TEST GRADIENT MODE

 Log probability=-nan

 param idx           value           model     finite diff           error

         0               0            -inf            -inf            -nan

         1               0         15.8735            -nan            -nan

goodrich@CYBERPOWERPC:/tmp/SU/rubi$ ./test --test_grad --data=dat.txt

Initialization failed after 100 attempts.  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

 

See attachments.

Ben

[rubi]

Thanks a lot for the advice.

I was thinking of trying to simplify my problem by assuming that my observations are normally distributed, and seeing what would be the proper model for my objective given these assumptions.

So I have O[g][s][k]~N(µ_gs,σ^2gs) (as before, g denotes group, s denotes subgroup, and k is the observation index)

What I believe to be the data generating process is that the different groups share variance (but not necessarily the means), and obviously the different subgroups share variance as well (but again not necessarily the means). So my objective is to estimate the variance shared among groups (or parameters of it).

Following Andrew Gelman's "Prior distributions for variance parameters in hierarchical models" (2005), I would have:

σ^2gs ~ Cauchy(0,As) (half Cauchy with a scale parameter A for subgroup s)

But since I have the group layer in the hierarchy I think I would like As to be a hyperparameter and a prior for it. Perhaps:

As ~ exponential(1) (or some empirical value instead of 1)

Does this make sense?

I also need to put a prior for µ_gs but I don't plan it to have a hyperparameter since I don't believe different subgroups and different groups share means, so I was thinking of:

µ_gs ~ N(0,1) (or some empirical values different from 0 and 1) or a uniform with empirical parameters.

Any advice would be greatly appreciated.

Best,

rubi

On Thursday, February 21, 2013 6:29:52 PM UTC-5, rubi wrote:

[Bob Carpenter]

After some discussion on the dev list, Matt correctly
guessed there was a zero in the data being modeled as
beta distributed, which I repeat here, with the zero in the
first line of data highlighted.

O <-
structure(c(0.7, 0.2, 0.55, 0.2, 0.5, 0, 0.42, 0.4, 0.3, 0.8, 0.74, 0.6, 0.72,
^
0.19, 0.58, 0.18, 0.48, 0.01, 0.45, 0.48, 0.4, 0.85, 0.75, 0.68, 0.69, 0.26,
0.54, 0.19, 0.55, 0.05, 0.44, 0.49, 0.35, 0.78, 0.71, 0.69, 0.68, 0.24, 0.52,
0.21, 0.51, 0.04, 0.46, 0.5, 0.38, 0.81, 0.7, 0.7, 0.71, 0.35, 0.51, 0.22, 0.52,
0.03, 0.43, 0.51, 0.39, 0.82, 0.72, 0.71),
.Dim = c(4, 3, 5))

I don't know there's much we can do to fix this issue. An actual 0 is
outside the proper support for the beta distribution.

- Bob

On 2/25/13 5:38 PM, Matt Hoffman wrote:
> Ah. Might the problem be that there's a 0 in the data? That 0 means you compute
> –(a – 1) * ∞
> which is big trouble for any value of a that's not 1. On the
> unconstrained scale, there's actually no valid setting of a. But even
> on the constrained scale, the prior assigns zero probability to a=1,
> so you're guaranteed that the log-density of the data will be either
> infinity or negative infinity. What happens if you replace the 0 with
> 0.001?






>
> Matt
>
> On Mon, Feb 25, 2013 at 2:27 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>> On Monday, February 25, 2013 4:57:42 PM UTC-5, Matt Hoffman wrote:
>>>
>>> But that's not ideal. I'm not sure what we can do about this within
>>> our current architecture, which really wants to transform beta r.v.s
>>> to [0, 1] and gamma r.v.s to [0, inf].
>>
>>
>> All true, but in this particular case, the RV was data and not being
>> transformed. So, it seems that for alpha very small and beta very large and
>> x close to zero, the gradient of the log density blows up. Possibly the beta
>> distribution is inappropriate for these data but it would be nice if we
>> could come up with a solution to our numerical problems.

>>
>> Ben
>>
>> --
>> You received this message because you are subscribed to the Google Groups

>> "stan development mailing list" group.

>> To unsubscribe from this group and stop receiving emails from it, send an

>> email to stan-dev+u...@googlegroups.com.

[rubi]

Thanks a lot again.

These data are only a sample of my real data (attached in just case along with the wrapper R script that parses it and calls the stan model), however I'm not sure that adding pseudocounts to avoid 0 observations will help.

Running stan in test_grad mode with init=0 (with the model at the bottom) gives this output:

TEST GRADIENT MODE

 Log probability=nan

 param idx           value           model     finite diff           error

         0               0        -3735.74             nan             nan

         1               0            -inf            -inf             nan

TESTING GRADIENT FOR MODEL 'test' NOW (CHAIN 2).

TEST GRADIENT MODE

 Log probability=nan

 param idx           value           model     finite diff           error

         0               0        -3735.74             nan             nan

         1               0            -inf            -inf             nan

TESTING GRADIENT FOR MODEL 'test' NOW (CHAIN 3).

TEST GRADIENT MODE

 Log probability=nan

 param idx           value           model     finite diff           error

         0               0        -3735.74             nan             nan

         1               0            -inf            -inf             nan

TESTING GRADIENT FOR MODEL 'test' NOW (CHAIN 4).

TEST GRADIENT MODE

 Log probability=nan

 param idx           value           model     finite diff           error

         0               0        -3735.74             nan             nan

         1               0            -inf            -inf             nan

I guess this means that a beta distribution for the observations doesn't really cut it.

I was wondering what you think about the model for normally distributed observations I wrote in the previous post?

Thanks a lot,

rubi

stan file used:

data {

  int<lower=1> G; // number of groups

  int<lower=1> S; // number of subgroups

  int<lower=1> K; // number of observations in each subgroup

  matrix<lower=0,upper=1>[S,K] O[G];

  matrix<lower=0,upper=1>[S,K] M[G];

}

parameters {

  real<lower=0> a;

  real<lower=0> b;

}

model {

   for (g in 1:G){

for (s in 1:S){

for (k in 1:K){ 

if(M[g,s,k] == 1){

O[g,s,k] ~ beta(a,b);

}

}

}

   }

}

[Ben Goodrich]

On Tuesday, February 26, 2013 12:27:20 AM UTC-5, rubi wrote:

Thanks a lot again.

These data are only a sample of my real data (attached in just case along with the wrapper R script that parses it and calls the stan model), however I'm not sure that adding pseudocounts to avoid 0 observations will help.

If the data have 0s (or 1s), then the plain beta distribution will always fail. But a normal distribution truncated to [0,1] will probably not work well either. There is some work on zero inflated beta models; for example see

http://arxiv.org/pdf/1103.2372.pdf

Ben

 

[M.F. Jonker]

Hi Ruby,

Another option could be to slightly re-scale your data prior to the estimation. See e.g. the section "2.2.1 Handling boundary observations" in the attached paper of Verkuilen and Smithson (2011).

Good luck!
Marcel

---

Date: Tue, 26 Feb 2013 03:36:24 -0800
From: goodri...@gmail.com
To: stan-...@googlegroups.com
Subject: Re: [stan-users] beta distribution

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[rubi]

Thanks a lot Marcel and Ben.

My idea for the normal distribution model is to use it on variance-preserving transformation of the data, making it normally distributed. My main dilemma is which prior to use for the top level hierarchy.

As for the beta distribution, I adopted the epsilon/2 + O[g][s][k]*((1-epsilon)/2) correction to rescale my observations which thus avoids the 0 and 1 boundaries.

I ran stan with the test_grad mode and was wondering whether I could get your help interpreting the output:

TESTING GRADIENT FOR MODEL 'test' NOW (CHAIN 1).

TEST GRADIENT MODE

 Log probability=0

 param idx           value           model     finite diff           error

         0               0        -6746.12        -6746.12    -3.78069e-07

         1               0          1970.6          1970.6    -1.44837e-07

TESTING GRADIENT FOR MODEL 'test' NOW (CHAIN 2).

TEST GRADIENT MODE

 Log probability=0

 param idx           value           model     finite diff           error

         0               0        -6746.12        -6746.12    -3.78069e-07

         1               0          1970.6          1970.6    -1.44837e-07

TESTING GRADIENT FOR MODEL 'test' NOW (CHAIN 3).

TEST GRADIENT MODE

 Log probability=0

 param idx           value           model     finite diff           error

         0               0        -6746.12        -6746.12    -3.78069e-07

         1               0          1970.6          1970.6    -1.44837e-07

TESTING GRADIENT FOR MODEL 'test' NOW (CHAIN 4).

TEST GRADIENT MODE

 Log probability=0

 param idx           value           model     finite diff           error

         0               0        -6746.12        -6746.12    -3.78069e-07

         1               0          1970.6          1970.6    -1.44837e-07

Thanks again,

rubi

[Bob Carpenter]

The output of test_grad contains the unconstrained params in the
order they are defined and their gradient. So with init=0, they
all have initial value 0. It's evaluating whether the auto-dif
gets the same thing as finite diffs -- those are the two columns
and their diff is the error.

- Bob

> ------------------------------------------------------------------------------------------------------------------------

> Date: Tue, 26 Feb 2013 03:36:24 -0800

> From: goodri...@gmail.com <javascript:>
> To: stan-...@googlegroups.com <javascript:>

> Subject: Re: [stan-users] beta distribution
>
> On Tuesday, February 26, 2013 12:27:20 AM UTC-5, rubi wrote:
>
> Thanks a lot again.
>
> These data are only a sample of my real data (attached in just case along with the wrapper R script that parses
> it and calls the stan model), however I'm not sure that adding pseudocounts to avoid 0 observations will help.
>
>
> If the data have 0s (or 1s), then the plain beta distribution will always fail. But a normal distribution truncated
> to [0,1] will probably not work well either. There is some work on zero inflated beta models; for example see
>

> http://arxiv.org/pdf/1103.2372.pdf <http://arxiv.org/pdf/1103.2372.pdf>

>
> Ben
>
> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com

> <javascript:>.
> For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

[Olive Buhule]

I am new to using stan but have been reading through other people's posts about an error message that am also  getting when I run the code below. Unfortanatelly, I have failed to fix the problem myself. In the code, I am trying to fit a joint model to a longitudinal outcome (y1), a binary outcome (y2)  and survival  outcome (survt), which is not observed directly. I am therefore treating it as a trasformed parameter because it involved parameters "delta" and random effects "b[, 4]. I am trying to fitt a piecewise exponetial  for the  survival component. Any help  to fix the error message "  stan::math::poisson_log: Rate parameter is -73.6109, but must be >= 0!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."   will be highly appreciated.

 

Thanks in advance,

 

Olive

 

 data{
int<lower=1> Nobs; // number of observations
real time[Nobs];
real y1[Nobs];
int<lower=-3> K;
int<lower=0> N; // number of subjects
int<lower=1,upper=N> id[Nobs]; // subject ID
real time2[N];
int<lower=0,upper=1>  y2[N];
int<lower=0,upper=1>  fail[N];
vector[4] zero; // a vector of zeros
int<lower=0> T;
real g;
real q;
}

parameters{
real delta; // fixeffs
vector[3] beta;// fixeffs
real alpha;// fixeffs
vector[3] gama;// fixeffs
vector[3] r;// fixeffs
vector[4] b[N]; // subj raneffs
real b5[N]; // subj raneffs
real<lower=0> sigma_e; // residual variance
vector<lower=0>[4] sigma_b; // subj variance
corr_matrix[4] Omega_b; // corr matrix for
// subj raneffs
vector<lower=0>[4] h0; // hazard for each interval

}

transformed parameters{
cov_matrix[4] Sigma_b; // varcov matrix for subj ranefs
real survt[N]; // survival time

for (j in 1:N){
survt[j] <- time2[j]-(delta + b[id[j],4]);

}

 
for(l in 1:4){
for(c in 1:4){
Sigma_b[l,c] <- sigma_b[l] * sigma_b[c] * Omega_b[l,c];

}
}
}

model{
real mu[Nobs];
real A[Nobs];
int idt[N,T]; 
int dN[N,T];
vector[T] D[N];

real  t[T+1];
real  maximum ;

for(i in 1:Nobs){
A[i] <- time[i]-(delta + b[id[i],4]);

mu[i] <-K;

if (A[i]==0){
mu[i] <- mu[i]; 

} else if (A[i]>0){
mu[i] <- mu[i] + beta[1]*(time[i]-(delta + b[id[i],4])) + beta[2]*((time[i]-(delta + b[id[i],4]))*(time[i]-(delta + b[id[i],4]))) + b[id[i],1]*(time[i]-(delta + b[id[i],4]))+ b[id[i],2]*((time[i]-(delta + b[id[i],4]))*(time[i]-(delta + b[id[i],4])));
}

else{
mu[i] <- mu[i] + beta[3]*(time[i]-(delta + b[id[i],4])) + b[id[i],3]*(time[i]-(delta + b[id[i],4]));
}  

}

y1 ~ normal(mu,sigma_e);
b ~ multi_normal(zero,Sigma_b);

 maximum<-max(survt);

 t[1]<-0;
 t[2]<-maximum/4;
 t[3]<-(maximum/4)*2;
 t[4]<-(maximum/4)*3;
 t[5]<- maximum/4;

for (j in 1:N){
 for (i in 1:T){
    dN[j,i]<- int_step(survt[j]-t[i])*fail[j]*int_step(t[i+1]-survt[j]);
   D[j,i]<-(fmin(survt[j], t[i+1])-t[i])*int_step(survt[j]-t[i]);
   }
 
}

 

for (j in 1:N){
 for (i in 1:T){
h0[i]~ gamma(q*g, g);
   //if (D[j,i]!=0)
 increment_log_prob(poisson_log(dN[j,i], h0[i]*D[j,i]*exp(gama[1]*b[id[j],1]+ gama[2]*b[id[j],2]+gama[3]*b[id[j],3])));
  }

}
 

for(j in 1:N){
y2 ~ bernoulli_logit(alpha + r[1]*b[id[j],1]+ r[2]*b[id[j],2]+r[3]*b[id[j],3]);
}

 

sigma_b ~ gamma(1.5,1.0E-4);
sigma_e ~ gamma(1.5,1.0E-4);
Omega_b ~ lkj_corr(2.0);
delta ~ normal(0, 100);
gama ~ normal(0, 100);
}

 

 

[Daniel Lee]

Hi Olive,

The error message has to do with line 86: increment_log_prob(poisson_log(dN[j,i], h0[i]*D[j,i]*exp(gama[1]*b[id[j],1]+ gama[2]*b[id[j],2]+gama[3]*b[id[j],3])));

The second argument to poisson_log, h0[i]*D[j,i]*exp(gama[1]*b[id[j],1]+ gama[2]*b[id[j],2]+gama[3]*b[id[j],3], is ending up negative; in your particular example, it was -73.6. There might be a bug in that expression or something that affects the calculation of the rate parameter for the poisson distribution.

Daniel

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Olive Buhule]

Thank you so much, Daniel!

I have decided to add a constrain and it seem to be running now. Here
is what I did

hl[j,i]<- h0[i]*D[j,i]*exp(gama[1]*b[id[j],1]+
gama[2]*b[id[j],2]+gama[3]*b[id[j],3]);
if (hl[j,i]>0)
increment_log_prob(poisson_log(dN[j,i], hl[j,i]));

Hope it makes sense.

Thanks again!

Olive

[Bob Carpenter]

That is almost certainly not what you want to do. It means
you'll just be ignoring a bunch of data in a way that is conditional
on a parameter value.

What Daniel was suggesting is that you do something like constrain
h0 > 0 and D > 0 in their declarations.

- Bob

[Olive Buhule]

Thank you so much , Bob! I have tried doing that but this is what I get:

failed to parse Stan model 'rstan-model-8-12-15v4' with error message:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:

require unconstrained. found range constraint.

ERROR at line 58

57: vector[T] hl[N];
58: vector<lower=0>[T] D[N];
^
59: real t[T+1];


Thanks,

Olive

[Jonah]

This is because constraints can't be declared in the model block. You'd have to do it in parameters or transformed parameters.

[Bob Carpenter]

You can't put constraints on local variables.

- Bob