Local Level Model Issue
188 views
Skip to first unread message
dynotec
Aug 22, 2012, 4:57:04 AM8/22/12
to stan-...@googlegroups.com
I'm trying to fit a simple local level model in Stan that I've attatched, and it compiles fine. But then I get output of:
OMPILING MODEL 'anon_model' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'anon_model' NOW.
In file included from /home//x86_64-redhat-linux-gnu-library/2.15/rstan/include/rstan/stan_fit.hpp:37,
from /home//x86_64-redhat-linux-gnu-library/2.15/rstan/include/rstan/rstaninc.hpp:4,
from file6ce37c460aea.cpp:6:
/home//x86_64-redhat-linux-gnu-library/2.15/rstan/include//stansrc/stan/io/dump.hpp: In member function ‘bool stan::io::dump_reader::scan_struct_value()’:
/home//x86_64-redhat-linux-gnu-library/2.15/rstan/include//stansrc/stan/io/dump.hpp:592: warning: comparison between signed and unsigned integer expressions
/home//x86_64-redhat-linux-gnu-library/2.15/rstan/include//stansrc/stan/io/dump.hpp:595: warning: comparison between signed and unsigned integer expressions
/home//x86_64-redhat-linux-gnu-library/2.15/rstan/include//stansrc/stan/io/dump.hpp:624: warning: comparison between signed and unsigned integer expressions
/home//x86_64-redhat-linux-gnu-library/2.15/rstan/include//stansrc/stan/io/dump.hpp:627: warning: comparison between signed and unsigned integer expressions
file6ce37c460aea.cpp: In member function ‘void anon_model_namespace::anon_model::write_array(std::vector<double, std::allocator<double> >&, std::vector<int, std::allocator<int> >&, std::vector<double, std::allocator<double> >&)’:
file6ce37c460aea.cpp:165: warning: unused variable ‘function__’
file6ce37c460aea.cpp:179: warning: unused variable ‘lp__’
file6ce37c460aea.cpp: In member function ‘virtual void anon_model_namespace::anon_model::write_csv(std::vector<double, std::allocator<double> >&, std::vector<int, std::allocator<int> >&, std::ostream&)’:
file6ce37c460aea.cpp:211: warning: unused variable ‘function__’
file6ce37c460aea.cpp:223: warning: unused variable ‘lp__’
/home//x86_64-redhat-linux-gnu-library/2.15/rstan/include//stansrc/stan/agrad/agrad.hpp: At global scope:
/home//x86_64-redhat-linux-gnu-library/2.15/rstan/include//stansrc/stan/agrad/agrad.hpp:2215: warning: ‘void stan::agrad::free_memory()’ defined but not used
/home//x86_64-redhat-linux-gnu-library/2.15/rstan/include//stansrc/stan/math/special_functions.hpp:585: warning: ‘double stan::math::log_sum_exp(const std::vector<double, std::allocator<double> >&)’ defined but not used
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Iteration: 2 / 5 [ 40%] (Adapting)
********************************************************************************
Error in model.grad_log_prob:
Dynamic exception type: std::domain_error
std::exception::what: Error in function stan::prob::uniform_log<N4stan5agrad3varE>(N4stan5agrad3varE): y is nan:0, but must not be nan!
Iteration: 3 / 5 [ 60%] (Sampling)
"
The set-up seems to run the schools example fine without any issues. Does anyone have any advice?
The set-up seems to run the schools example fine without any issues. Does anyone have any advice?
Marcus Brubaker
Aug 22, 2012, 10:42:05 AM8/22/12
to stan-...@googlegroups.com
The compiler warnings can be ignored safely. The domain_errors mean
the sampler tried a parameter value which resulted in a bad value
somewhere along the line. This is normal, particularly in the
adaptation stage.
Cheers,
Marcus
the sampler tried a parameter value which resulted in a bad value
somewhere along the line. This is normal, particularly in the
adaptation stage.
Cheers,
Marcus
dynotec
Aug 23, 2012, 3:38:57 AM8/23/12
to stan-...@googlegroups.com
The issue is that I get domain errors at every iteration, as well as non-sense parameter estimates even after very long runs on multiple chains on simulated chains that disagree with estimates obtained from BUGS(Giving estimated values above 1 or below zero when the range of observations is between .48 and .51). Is there an error with the model code that is causing this?
Bob Carpenter
Aug 23, 2012, 4:48:51 AM8/23/12
to stan-...@googlegroups.com
The issue is that the model
------------
data {
int[0,] J; // number of days
real y[J]; // Polls
}
parameters {
real alpha[J];
real[0,.1] tau;
}
model {
alpha[1]~uniform(0,1);
for(i in 2:J)
alpha[i]~normal(alpha[i-1],tau);
for(i in 1:J)
y[i]~normal(alpha[i],0.01);
}
---------------
tries to sample
alpha[1] ~ uniform(0,1);
but doesn't constrain alpha[1] to fall between 0 and 1.
So what happens is that in initialization or when
randomizing momentum in each iteration, alpha[1] is going
outside its support in the 0 to 1 range, and thus
uniform() returns not-a-number (NaN).
Fix 1.
Put a different prior on alpha[1].
For instance,
alpha[1] ~ normal(0,100);
If you just drop the sampling statement for
alpha[1], you implicitly have an improper uniform
prior on alpha[1]. Because alpha[1] is constrained
through the sampling statement for alpha[2] in the
loop, this might work. In general, the posterior
must be proper.
Fix 2.
Define alpha[1] separately with a constraint.
real<lower=0,upper=1> alpha_1;
and then make alpha one element shorter:
real alpha[J-1];
and then you can sample
alpha_1 ~ uniform(0,1);
We've changed the constraint syntax on you since
your first post. This is probably the last syntax change
before the 1.0 release, and we'll try very hard to
maintain backward compatibility after that.
Then you'll need to change the rest of the model:
for(i in 2:J)
alpha[i]~normal(alpha[i-1],tau);
for(i in 1:J)
y[i]~normal(alpha[i],0.01);
to something where the new indexing is one less
alpha[1] ~ normal(alpha_1,tau);
for (i in 2:(J-1))
alpha[i] ~ normal(alpha[i-1],tau);
y[1] ~ normal(alpha_1,0.01);
for (i in 2:J)
y[i] ~ normal(alpha[i-1],0.01);
This kind of loop unfolding and consequent index
fiddling's a pain, but we don't have a way to
define a constraint on a single member of
an array.
An only slight less efficient alternative would be to
define a derived
real<lower=0,upper=1> alpha_1
real alpha_raw[J-1];
and then a transformed parameter
transformed parameters {
real alpha[J];
alpha[1] <- alpha_1;
for (j in 2:J)
alpha[j] <- alpha_raw[j-1];
}
Then you can use alpha as before.
- Bob
On 8/23/12 3:38 AM, dynotec wrote:
> The issue is that I get domain errors at every iteration, as well as non-sense parameter estimates even after very long
> runs on multiple chains on simulated chains that disagree with estimates obtained from BUGS(Giving estimated values
> above 1 or below zero when the range of observations is between .48 and .51). Is there an error with the model code that
> is causing this?
>
> On Wednesday, August 22, 2012 9:42:05 AM UTC-5, Marcus wrote:
>
> The compiler warnings can be ignored safely. The domain_errors mean
> the sampler tried a parameter value which resulted in a bad value
> somewhere along the line. This is normal, particularly in the
> adaptation stage.
>
> Cheers,
> Marcus
>
------------
data {
int[0,] J; // number of days
real y[J]; // Polls
}
parameters {
real alpha[J];
real[0,.1] tau;
}
model {
alpha[1]~uniform(0,1);
for(i in 2:J)
alpha[i]~normal(alpha[i-1],tau);
for(i in 1:J)
y[i]~normal(alpha[i],0.01);
}
---------------
tries to sample
alpha[1] ~ uniform(0,1);
but doesn't constrain alpha[1] to fall between 0 and 1.
So what happens is that in initialization or when
randomizing momentum in each iteration, alpha[1] is going
outside its support in the 0 to 1 range, and thus
uniform() returns not-a-number (NaN).
Fix 1.
Put a different prior on alpha[1].
For instance,
alpha[1] ~ normal(0,100);
If you just drop the sampling statement for
alpha[1], you implicitly have an improper uniform
prior on alpha[1]. Because alpha[1] is constrained
through the sampling statement for alpha[2] in the
loop, this might work. In general, the posterior
must be proper.
Fix 2.
Define alpha[1] separately with a constraint.
real<lower=0,upper=1> alpha_1;
and then make alpha one element shorter:
real alpha[J-1];
and then you can sample
alpha_1 ~ uniform(0,1);
We've changed the constraint syntax on you since
your first post. This is probably the last syntax change
before the 1.0 release, and we'll try very hard to
maintain backward compatibility after that.
Then you'll need to change the rest of the model:
for(i in 2:J)
alpha[i]~normal(alpha[i-1],tau);
for(i in 1:J)
y[i]~normal(alpha[i],0.01);
to something where the new indexing is one less
alpha[1] ~ normal(alpha_1,tau);
for (i in 2:(J-1))
alpha[i] ~ normal(alpha[i-1],tau);
y[1] ~ normal(alpha_1,0.01);
for (i in 2:J)
y[i] ~ normal(alpha[i-1],0.01);
This kind of loop unfolding and consequent index
fiddling's a pain, but we don't have a way to
define a constraint on a single member of
an array.
An only slight less efficient alternative would be to
define a derived
real<lower=0,upper=1> alpha_1
real alpha_raw[J-1];
and then a transformed parameter
transformed parameters {
real alpha[J];
alpha[1] <- alpha_1;
for (j in 2:J)
alpha[j] <- alpha_raw[j-1];
}
Then you can use alpha as before.
- Bob
On 8/23/12 3:38 AM, dynotec wrote:
> The issue is that I get domain errors at every iteration, as well as non-sense parameter estimates even after very long
> runs on multiple chains on simulated chains that disagree with estimates obtained from BUGS(Giving estimated values
> above 1 or below zero when the range of observations is between .48 and .51). Is there an error with the model code that
> is causing this?
>
> On Wednesday, August 22, 2012 9:42:05 AM UTC-5, Marcus wrote:
>
> The compiler warnings can be ignored safely. The domain_errors mean
> the sampler tried a parameter value which resulted in a bad value
> somewhere along the line. This is normal, particularly in the
> adaptation stage.
>
> Cheers,
> Marcus
>
dynotec
Aug 23, 2012, 1:43:25 PM8/23/12
to stan-...@googlegroups.com
Thank you for the quick response! I ended up just restricting the entire vector to [0,1] since I'm working with proportions anyway.
I'm having issues using the mean function. The root issue seems to be when I define fixed effects in the model block:
for(i in 1:NORG)
house[i]~normal(0,HouseSigma);
and then try to normalize them in the transformed parameter block before the model uses them:
for (j in 1:NORG)
house2[j]<-house[j]-housemean;
and then the model uses them as
##Measurement Model
for(i in 1:NPOLLS)
y[i]~normal(alpha[i]+house2[org[i]]+alpha[date[i]] ,
sqrt(alpha[date[i]]*(1-alpha[date[i]])/SampleSize[i])*PIE1)
;
.
The model runs fine if I don't normalize the fixed effects, but when I try it with them, the model returns pretty lengthy error output and doesn't generate any iterations. This sort of code works on BUGS, does the mean function cause problems when used this way in STAN?
Bob Carpenter
Aug 23, 2012, 2:38:53 PM8/23/12
to stan-...@googlegroups.com, stan...@googlegroups.com
On 8/23/12 1:43 PM, dynotec wrote:> Thank you for the quick response! I ended up just restricting the entire vector to
keep all the values restricted to [0,1].
You might want to switch to our newer syntax, which
will go out in version 1.0. It uses named arguments
for the bounds, but is otherwise the same.
> I'm having issues using the mean function. ....
Thanks for the report. Your model should work as
written. I'll fix it and push a patch.
The key line in the error output is the following, which
says we're missing an operator+=() specialization or
have called something wrong:
...
include//stansrc/stan/math/matrix.hpp:731: error: no match for ‘operator+=’
in ‘sum += ((const std::vector<stan::agrad::var, std::allocator<stan::agrad::var> >*)v)
->std::vector<_Tp, _Alloc>::operator[]
[with _Tp = stan::agrad::var, _Alloc = std::allocator<stan::agrad::var>](i)’
- Bob
[0,1] since I'm working with
> proportions anyway.
Well that's the right thing to do, then. It'll
> proportions anyway.
keep all the values restricted to [0,1].
You might want to switch to our newer syntax, which
will go out in version 1.0. It uses named arguments
for the bounds, but is otherwise the same.
> I'm having issues using the mean function. ....
Thanks for the report. Your model should work as
written. I'll fix it and push a patch.
The key line in the error output is the following, which
says we're missing an operator+=() specialization or
have called something wrong:
...
include//stansrc/stan/math/matrix.hpp:731: error: no match for ‘operator+=’
in ‘sum += ((const std::vector<stan::agrad::var, std::allocator<stan::agrad::var> >*)v)
->std::vector<_Tp, _Alloc>::operator[]
[with _Tp = stan::agrad::var, _Alloc = std::allocator<stan::agrad::var>](i)’
- Bob
Reply all
Reply to author
Forward
0 new messages