Strange bias in toy model

[James Savage]

Hi all --

I might be trying to do something impossible, but I would like to use 

random data generated in the transformed data block--draws used 

for monte carlo integration within my model. I'm currently using draws 

from R and passing them in as data, so thought I might be able to 

save time by doing it in the model. 

Unsure whether I could use _rng in the transformed parameter block, 

I coded up the following toy example: 

data {

  int T;

}

transformed data {

  vector[T] out;

  

  for(t in 1:T) {

    out[t] = normal_rng(0 , 1);

  }

}

parameters {

  real mu;

  real<lower = 0> sig;

}

model {

  out ~ normal(mu, sig);

}

It compiles and runs fine. But it gives me some really strangely 

biased estimates for mu. When the sample size T is low, the 

estimates of mu are centered on .3. I ran the attached file a couple

of times and got back precisely the same chart. 

It makes me think that the seed for random number generation in 

the transformed data block is fixed. Might this be the case? 

Jim

[Andrew Gelman]

I've never seen normal_rng in the transformed data block.  I thought it could only go in generated quantities.

A

--
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.
<toy_model_rng_in_transformed_data.R>

[James Savage]

To be clear, I thought it wasn't allowed and just coded up this toy model to see if it worked. 

[Bob Carpenter]

Yes, it's allowed, but we have an outstanding issue with it that
it uses different RNGs in each chain, so you can't expect
convergence according to Rhat until we fix that:

https://github.com/stan-dev/stan/issues/2014

- Bob

> On Mar 6, 2017, at 6:03 PM, James Savage <jav...@gmail.com> wrote:
>

[Bob Carpenter]

Hmm. Looks like that bug may have been fixed, leading to
the unexpected behavior you're seeing.

You have to be careful about what you mean by biased. If
I generate 100 draws from normal(0, 1), it might have a mean
of 0.2. The estimator that then estimes a mean of 0.2 isn't
biased.

It looks like the bug is exactly the opposite of what we
listed now --- it doesn't even take the seed into account.
I'll amend the issue.

But it does look like everything's working as intended as far as RNGs
and samplers are concerned. I amended your program to calculate
the sample mean and sd of the pseudo-random numbers and get this:

data {
int<lower=0> N;
}
transformed data {
vector[N] y;
for (n in 1:N)
y[n] = normal_rng(0, 1);
}
parameters {
real mu;
real<lower = 0> sigma;
}
model {
y ~ normal(mu, sigma);
}
generated quantities {
real mean_y = mean(y);
real sd_y = sd(y);
}

When I set N=100, I indeed see that the RNG is generating a sample
on the high side of average:

N=100:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu 0.20 0.00 0.10 -0.01 0.13 0.20 0.27 0.41 4000 1
sigma 1.03 0.00 0.08 0.90 0.98 1.03 1.08 1.19 3676 1
mean_y 0.20 0.00 0.00 0.20 0.20 0.20 0.20 0.20 2 1
sd_y 1.02 0.00 0.00 1.02 1.02 1.02 1.02 1.02 2 1
lp__ -52.42 0.02 1.01 -55.20 -52.82 -52.09 -51.69 -51.44 1796 1


N=1000:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu 0.01 0.00 0.03 -0.05 -0.01 0.01 0.03 0.07 3411 1
sigma 0.99 0.00 0.02 0.94 0.97 0.98 1.00 1.03 3672 1
mean_y 0.01 0.00 0.00 0.01 0.01 0.01 0.01 0.01 2 1
sd_y 0.98 0.00 0.00 0.98 0.98 0.98 0.98 0.98 2 1
lp__ -485.32 0.03 1.02 -488.09 -485.70 -484.98 -484.61 -484.35 1485 1

- Bob


> On Mar 6, 2017, at 6:03 PM, James Savage <jav...@gmail.com> wrote:
>

[Bob Carpenter]

Here's the new issue:

https://github.com/stan-dev/stan/issues/2241

I closed the old one. Thansk for reporting this. We'll
fix it with the refactor in the next release.

- Bob

[James Savage]

Thanks Bob - that makes sense.