error standardizing in transformed parameter block

[atzap]

Hi all

Thanks for your help with my last problem.  I've managed to get past the compiling stage but have predictably run into a new problem.  Unfortunately, I can't seem to figure this one out and reckon other folk might run into this as well...

I'm trying to standardized some parameters to a mean of zero and sd of one for use later in my model.  I've tried doing this as follows (this is a simplified version of my model):

data{
 int<lower=0> N;                  
 int<lower=0> y[N];                
 int<lower=0> y2[N];                
 int<lower=0> y3[N];            
 int<lower=0> n[N];  
 real CA[N];          
}
parameters{
 real lambda[3];
}
transformed parameters{
real<lower = 0, upper = 1> pBound[N];
real<lower = 0, upper = 1> pBound2[N];
real<lower = 0, upper = 1> pBound3[N];
real pBound2_centered[N];
real pBound3_centered[N];
  for(i in 1:N){
   pBound[i] <- fmax(0, fmin(1, inv_logit(lambda[1]) ) );
   pBound2[i] <- fmax(0, fmin(1, inv_logit(lambda[2] * CA[i]) ) );
   pBound3[i] <- fmax(0, fmin(1, inv_logit(lambda[3] * CA[i]) ) );
  }
  for(i in 1:N){
    pBound2_centered[i] <- (pBound2[i] - mean(pBound2))/sd(pBound2);
    pBound3_centered[i] <- (pBound3[i] - mean(pBound3))/sd(pBound3);
  }
}
model{
  for (i in 1:N){
   y[i] ~ binomial(n[i], pBound[i]);
   y2[i] ~ binomial(n[i], pBound2[i]);
   y3[i] ~ binomial(n[i], pBound3[i]);
  }
  lambda ~ normal(0, 1E3);
}

The model compiles but then fails to sample, calling the following error:

Error in sampler$call_sampler(args_list[[i]]) : std::bad_alloc

I can replace:

pBound2_centered[i] <- (pBound2[i] - mean(pBound2))/sd(pBound2)

With:

pBound2_centered[i] <- (pBound2[i] - mean(pBound2))

And the model gets sampled.  But, it doesn't work if I try and replace this with:

pBound2_centered[i] <- (pBound2[i] - sd(pBound2))

So, I have two questions.

1.  What is the error specifically telling me in this case?

2.  What is the proper way to standardize a parameter in STAN?  I note that I tried the "stacks" example from BUGS (see http://code.google.com/p/stan/source/browse/src/models/bugs_examples/vol1/stacks/stacks_a_normal.stan), where standardization is done in the transformed data block, and it runs fine.  So clearly, I'm missing something with respect to trying to standardize parameters. 

Many thanks for your help and patience,

Andrew

[atzap]

Sorry, I should have added that I not only replace "pBound2_centered" but also "pBound3_centered" in each example.

[Bob Carpenter]

Let me reiterate that it's much easier for us to
help you debug if you send us two things as attachments:

1. The model code.

2. Data file that works with the model.

We want to help, but it's very hard without having
an example we can run.

- Bob

[Bob Carpenter]

On 10/2/12 2:47 PM, atzap wrote:
> Hi all
>
> Thanks for your help with my last problem. I've managed to get past the compiling stage but have predictably run into a
> new problem. Unfortunately, I can't seem to figure this one out and reckon other folk might run into this as well...
>
> I'm trying to standardized some parameters to a mean of zero and sd of one for use later in my model. I've tried doing

> this as follows (/this is a simplified version of my model/):

Let me try to make clearer what I mean by simplified.
I mean reduce the model to the point where it's as simple
as possible and still causes the error. For instance,
if you're having an issue with normalizing, write a model
that just does the normalization step.

That is, take something like this:

To simplify, you could reduce to only one of y, y2 and y3 and
pBound, pBound2 and pBound3.

It's not doing the same thing, but it should have the same
error behavior and it's MUCH easier to debug smaller
models because you can easily see what each piece is
doing.

The second way to simplify is to reduce the data set size.

This isn't BUGS -- we've taken more care in our inv_logit
implementation so you don't need to do things like:

> pBound[i] <- fmax(0, fmin(1, inv_logit(lambda[1]) ) );

Just

pBound[i] <- inv_logit(lambda[1]);

should suffice.

It's also not BUGS in that things like this are inefficient
in Stan:

> for(i in 1:N){
> pBound2_centered[i] <- (pBound2[i] - mean(pBound2))/sd(pBound2);
> pBound3_centered[i] <- (pBound3[i] - mean(pBound3))/sd(pBound3);
> }

It's computing the mean and sd every time through the
loop! Change this to this:

real mean_pbound_2 = mean(pBound2);
real sd_bound_2 = sd(pBound2);
for (n in 1:N)
pBound2_centered[n] <- (pBound2[i] - mean_pbound_2) / sd_bound_2;

Then, move all of this to the derived quantities block,
because nothing depends on pBound2_centered in the model.
It'll be MUCH faster with these two changes.

Having said all this, there's no reason to believe that
pBound2_centered[n] will remain between 0 and 1. If
sd_bound_2 is real low, then this value can be arbitrarily
big.

And do the same for pBound3. It'll be MUCH faster, because
the transformed params get executed each time.

I'll wait until you send us data to try your model. See
the last message I posted on asking for help.

- Bob

[atzap]

Dear Bob

I'm terribly sorry for the useless example.  I thought it would be easier to follow than my much more convoluted model, but I now understand that trying to find short-cuts just ends up wasting your time.  Again, I'm sorry.

I've tried to simplify my model by focusing on the normalization step.  I wasn't sure if I could remove the rest of what I've included, because I thought it was important to show the context of how the normalization is being used.  I apologize in advance if there are any errors in the code.  By the way, the efficiencies you've pointed out definitely help, thanks!

So, I am still finding the same problem when I try and calculate the standard deviation.  I attach the  model that runs, where I don't calculate a sd in the transformed parameters block (see model1_no_sd.R).  Where I've added the sd step, I'm now getting the following error (see model2_with_sd.R): 

Error in sampler$call_sampler(args_list[[i]]) : 

  Error in function stan::prob::binomial_log(N4stan5agrad3varE): Probability parameter is 1.#QNAN:0, but must be finite!

Presumably this is because something goes awry with the sd?  The only thing I can think of is that the sd is somehow hitting zero and so the normalization is throwing up infinite values...  I tried specifying initial values for the parameter for which I am calculating the sd, but got the same error (these are included in the model2 script).

Hopefully I've done a better job of following the instructions for help.

Andrew  

[atzap]

I've now tried various ways of bounding the sd estimate, such as through the "real" definition and using the fmax function, but this doesn't help.  

I've also tried manually calculating the sd rather than relying on the internal function.  I've added this to my model2_with_sd.R:

transformed parameters{
...
real df1[N];
for(i in 1:N){
    df1[i] <- pow((foodavail[i] - mean_fooda),2);
  }
  sd_fooda <- sqrt(sum(df1) * (1/(N-1)));
...
}

This still throws up:

SAMPLING FOR MODEL 'm3' NOW (CHAIN 1).

Error in sampler$call_sampler(args_list[[i]]) : 

  Error in function stan::prob::binomial_log(N4stan5agrad3varE): Probability parameter is 1.#QNAN:0, but must be finite!

I even get this when I further simplify the manual SD calculation to just sd_fooda <- sum(df1).  Am I missing something really obvious here?  I don't understand how all values of df1 are going to zero!

Thanks in advance for your help and patience.

[Daniel Lee]

On Wed, Oct 3, 2012 at 11:16 AM, atzap <ata...@gmail.com> wrote:
> I've now tried various ways of bounding the sd estimate, such as through the
> "real" definition and using the fmax function, but this doesn't help.
>
> I've also tried manually calculating the sd rather than relying on the
> internal function. I've added this to my model2_with_sd.R:
>
> transformed parameters{
> ...
> real df1[N];
> for(i in 1:N){
> df1[i] <- pow((foodavail[i] - mean_fooda),2);
> }
> sd_fooda <- sqrt(sum(df1) * (1/(N-1)));
> ...
> }
>
>
>
> This still throws up:
>
> SAMPLING FOR MODEL 'm3' NOW (CHAIN 1).
> Error in sampler$call_sampler(args_list[[i]]) :
> Error in function stan::prob::binomial_log(N4stan5agrad3varE): Probability
> parameter is 1.#QNAN:0, but must be finite!
>
>
> I even get this when I further simplify the manual SD calculation to just
> sd_fooda <- sum(df1). Am I missing something really obvious here? I don't
> understand how all values of df1 are going to zero!
>
> Thanks in advance for your help and patience.

Those values aren't 0's. That's how NaN's are being printed. That
message is indicating that the probability parameter parameter getting
into binomial_log() is 1) an auto-dif variable (either defined in the
parameter or transformed parameter block) and 2) that the value is
NaN, so something went off the rails before that call.

[atzap]

Aware that something went off before calling the binomial_log.  

If you have a chance to look at my model, you'll see that I think the cause is I'm dividing by zero to derive foodavail_centered, which then goes into the probability parameter.  So, I think the zero denominator is the cause.  I just don't understand why that denominator (equal to sd_fooda) is going to zero...  I suspect it is something to do with my butchering of the STAN language.

[Daniel Lee]

Thanks for your patience. I don't think it's a zero denominator. I'll
try to get to why we're seeing NaN's for pBound when we start
sampling.

Daniel

[Daniel Lee]

This issue was fixed in the last couple weeks within Stan. The fix
will be available with our next release.



Daniel