constrain for quadratic function > 0 not working

[Jan]

Hello,

for my gamma(a,b) I am fitting two quadratic functions

a = aa1 * x + aa2 * x^2 + aa3
b = ...

I need them to be > 0 to ensure stability for stan (a, b must be > 0)

However I get the following messages. I attach data and script.

model="
data{
int N;
int P;
vector[P] slices;
vector[P] slices2;
vector[N] y[P];

}

parameters{
  real<lower=0> cc1;
  real<lower=0> cc2;
  real<lower=0> aa2;
  real<lower=0> bb2;
  real<upper=sqrt(4*aa2*cc1)> aa1;
  real<upper=sqrt(4*bb2*cc2)> bb1;

}
transformed parameters{
  vector[P] a;
  vector[P] b;

  a = aa1 * slices + aa2 * slices2 + cc1;
  b = bb1 * slices + bb2 * slices2 + cc2;
}
model{
  cc1 ~ cauchy(0,2.5);
  aa2 ~ cauchy(0,2.5);
  bb2 ~ cauchy(0,2.5);

 for(p in 1:P) target += gamma_lpdf(y[p] | a[p],b[p]);    
}
"

fit <- rstan:::stan(
  model_code=model,
  data = list(N=dim(df)[1], P=length(slices), y=t(df), slices=slices, slices2=slices^2),
  iter=500
)

The following numerical problems occured the indicated number of times after warmup on chain 4 count Exception thrown at line 32: gamma_log: Shape parameter is 0, but must be > 0! 3 Exception thrown at line 32: gamma_log: Shape parameter is -inf, but must be > 0! 2 Exception thrown at line 32: gamma_log: Inverse scale parameter is -0.000696615, but must be > 0! 1 Exception thrown at line 32: gamma_log: Inverse scale parameter is -0.00276573, but must be > 0! 1 Exception thrown at line 32: gamma_log: Inverse scale parameter is -0.0100321, but must be > 0! 1 Exception thrown at line 32: gamma_log: Inverse scale parameter is -0.0104653, but must be > 0! 1 Exception thrown at line 32: gamma_log: Inverse scale parameter is -0.056854, but must be > 0! 1 Exception thrown at line 32: gamma_log: Inverse scale parameter is -0.072508, but must be > 0! 1 Exception thrown at line 32: gamma_log: Inverse scale parameter is -0.473384, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.000634826, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.00295345, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.00355422, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.00904165, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.0117633, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.0135997, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.01663, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.0234932, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.0289879, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.0345296, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.0355037, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.0419938, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.0437273, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.0569922, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.0701146, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.089961, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.099161, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.1152, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -0.533512, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -161.842, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -1.78671e+214, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -1.81511e+12, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -2.42762e-08, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -3.02458e-50, but must be > 0! 1 Exception thrown at line 32: gamma_log: Shape parameter is -6.06167e+52, but must be > 0! 1 When a numerical problem occurs, the Metropolis proposal gets rejected. However, by design Metropolis proposals sometimes get rejected even when there are no numerical problems. Thus, if the number in the 'count' column is small, do not ask about this message on stan-users.

Thanks

[Daniel Lee]

It looks like your computation for a and b aren't staying positive. Did you intend for aa1 and bb1 to have upper bounds instead of lower bounds?

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.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jan]

According to the formula (from e.g., http://www.sosmath.com/algebra/quadraticeq/signquadra/signquadra.html)

aa2 > 0 

aa1^2 + 4ac < 0

now.. why it does not work for my case?

[Jan]

Sorry for all positive values (http://www.biology.arizona.edu/biomath/tutorials/quadratic/roots.html)

aa1^2 - 4ac < 0

aa2 > 0

and is coded correctly (I believe in my code).

But it seems to keep having exceptions.

Thanks for your help!

[Daniel Lee]

Hi Jan,

You should print a and b to verify that all the elements are positive since you're using them for the gamma distribution. I think you'll find that aa1's magnitude is large enough to force a to be negative. 

Daniel