Re: [stan-users] Problem with truncated data model

[Jiqiang Guo]

Thanks for pointing this out. This might be bug in a stanc as the C++ code for the log-posterior has

      lp__ += stan::prob::normal_log<true>(get_base1(y,n,"y",1), mu, sigma);

      lp__ += log(normal_p(U, mu, sigma));

I think the second += should be -=.  So I am sending this to the stan-dev list. 

--

Jiqiang

On Mon, Sep 3, 2012 at 4:55 PM, Marek Smigielski <marek.sm...@gmail.com> wrote:

Hi,

I have some strange problem with truncated data model from user guide.  It is not producing expected estimation of any parameters neither mu nor sigma. My testing procedure is as follows.

1. Prepare random vector of observations from normal random generator with parameters (mean=0, sigma =1).

2. Next truncate data when it is greater than 0.

3. Compile and fit model.

4. Instead of expecting estimations near real values I get mu=-0.8, sigma=0.5.

My code in R is:

library(rstan) 

model_code <- '

data {

int<lower=0> N;

real U;

real<upper=U> y[N];

}

parameters{

real mu;

real<lower=0> sigma;

}

model {

for (n in 1:N)

  y[n] ~ normal(mu,sigma) T[,U];

}

'

real_data<- rnorm(200)

truncated_data<-real_data[which(real_data<0)]

data <- list(N = length(truncated_data),U=0, y = truncated_data)

fit <- stan(model_code = model_code, data = data)

print(fit)

Could you suggest what am I doing wrong or how can I achieve better estimations?

regards,

Marek Smigielski

[Bob Carpenter]

Good call. Totally my bad. And it affects all of
the truncated distributions, not just the upper-bound
only. None of these will work properly until we release
a patched version of Stan (and RStan), which we will do ASAP.

I just pushed a patch to our code repository that
fixes the issue, which is in stan/gm/generator.hpp,
starting on line 1024 or so.

We should set up a model test as the generator isn't
tested otherwise.

I forgot the negation due to the CDF terms being in
the denominator under a log.

An upper-truncated sampling statement

y ~ foo(...) T[,U];

should translate to the C++:

lp__ += foo_log(y,...);
lp__ -= log(foo_cdf(U,...));

Because I inadvertently made the second line here += instead
of -= in the generator, we'll need to test this
through a model.

The other two cases we have are with just a lower
bound, using log1m(x) == log(1 - x), to set

"" T[L,]
lp__ -= log1m(foo_cdf(L,...));

and with both an upper and lower bound,

"" T[L,U]
lp__ -= log(foo_cdf(U,...) - foo_cdf(L,...));

An earlier thread brings up the issue of whether we
should try to be more careful about underflow of the CDF
to 0 or overflow to 1. We could definitely be much
more careful with the arithmetical precision here.

We could really use three compound log-difference CDF
functions for these three cases.

We could then test if the lower bound was less than
the upper bound, something we're currently not doing.
As is, if they're backwards, you'll try to take the log
of a negative number and get an error that way.

- Bob

> --
>
>

[Nigel/essence]

Never take the log of a normal. The normal can get very small and be set to zero and/or underflow.

Use a specialised algorithm to get log normals.

[forgive me if I am preaching to the converted]

[Marcus Brubaker]

Would it be better to instead of defining

foo_cdf(x)

for each distribution defining a two parameter function

foo_logcdf(L,U)

to avoid issues of over/underflow where needed? By default we can
implement it in terms of foo_cdf(x) but we could also overload it as
needed to be more clever.

Cheers,
Marcus

> --
>
>

[Nigel/essence]

I posted this in another thread, so it may have got lost:

Garsaglia has done some interesting work in this area 'Evaluating the Normal Distribution' in JSS.

But even that was nowhere near enough for my tails, so I used the concept and wrote my own approximation. I was more interested in speed than accuracy.

You can get some ideas from that paper how to get log cum norm approximations at extreme tails.

[Marcus Brubaker]

The interface I was proposing would make it plausible to do this at some point.

> --
>
>

[Bob Carpenter]

n 9/4/12 11:18 AM, Marcus Brubaker wrote:
> Would it be better to instead of defining
>
> foo_cdf(x)
>
> for each distribution defining a two parameter function
>
> foo_logcdf(L,U)

Exactly.

I'm thinking we'll either need need three functions, as in

foo_cdf_lb_log(L|...) = log(1 - foo_cdf(L|...))
foo_cdf_ub_log(U|...) = log(foo_cdf(U|...) - 0)
foo_cdf_lub_log(L,U|...) = log(foo_cdf(U|...) - foo_cdf(L|...))

or some way to just use foo_cdf_lub_log(L,U|...)
with -inf and +inf bounds.

I'm reading up on some implementation suggesstions in
the JSS paper Nigel suggested.

- Bob

[Bob Carpenter]

On 9/4/12 10:40 AM, Nigel/essence wrote:
> Never take the log of a normal. The normal can get very small and be set to zero and/or underflow.

Absolutely.

In practice, overflow or underflow in the truncated
distributions shouldn't be much cause for concern. There we have
terms

-log(U - L)
-log(1 - L)
-log(U - 0)

where U and L are the cdf values of the upper and lower
bounds. If these overflow to 1 or underflow to 0, there's
not much harm done :-)

The one area where it could matter is if an upper bound
underflowed or a lower-bound overflowed. That would mean
you have a very very narrow interveral for the constraint.

> Use a specialised algorithm to get log normals.

We already do that for the probability functions and
many other functions.

As I mentioned last time, the frustrating part is that
even if the basic functions are OK, their interactions
can get messy. For instance, we added a log Bernoulli
distribution function that took logit-scaled parameters
because applying inv_logit would underflow or overflow
to 0 or 1 in some cases with initial parameters and lead
to overall 0 probability. (We truncated that to make
sure it didn't happen, but that's a hack, though one
Andrew's written about before under the heading 'Pinocchio
principle') -- a truncated inv_logit that's bounded between
(.01, .99) effectively handles overflow not only of the
arithmetic (for which (1E-15,1-1E-15) suffices.

> [forgive me if I am preaching to the converted]

I like listening to a good sermon on arithmetic
accuracy :-) And thanks for the reference. I'll read
it tonight.

- Bob

[Daniel Lee]

On Tue, Sep 4, 2012 at 6:05 PM, Ben Goodrich <goodri...@gmail.com> wrote:
> On Tuesday, September 4, 2012 3:15:09 PM UTC-4, Bob Carpenter wrote:
>>
>> If you only see a few, as you did here, it's not
>> a problem. We're also working on making that
>> clear by changing "Error" to "Warning" in the report.
>> (What would be really nice is to link them back to
>> a line of the Stan source, but that's MUCH trickier
>> and not on our short-term to-do list.)
>
>
> Can you metaprogram a call to git grep -n $WARNING inside the C++ source?
>

No. git is not on the list of requirements for Stan.

What would that line accomplish? Perhaps we can find another way to do
the same thing.

[Bob Carpenter]

On 9/4/12 5:20 PM, Jiqiang Guo wrote:
> On Tue, Sep 4, 2012 at 3:15 PM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:
> Because the truncation point is data, it's not
> actually doing anything in this model! The data
> are all known and constrained to be < U, so
> T[,U] terms just add some constants of proportionality.
> No, it does not just add some constants. Because the denominator of the pdf
> of truncated function is
>
> normal_cdf(U, mu, sigma),
>
> which involves parameter mu and sigma.

Thanks for the clarification. If it didn't work
as you just pointed out, we'd never be able to
do censored or truncated models!

Daniel and I are in the middle of debugging this
and having some problems getting simple truncated
estimates right. The only problem is we're out of
things to test -- the model's getting translated
properly, all the gradients for the normal density
and CDF seem right, too.

I wrote a model that includes the censored
data:

data {
real U;
int<lower=0> N_censored;
int<lower=0> N_observed;
real<upper=U> y[N_observed];
}
transformed data {
real<lower=0> sigma;
sigma <- 1.0;
}
parameters {
real mu;
}
model {
for (n in 1:N_observed)

y[n] ~ normal(mu,sigma) T[,U];

lp__ <- lp__ + N_censored * log(1 - normal_p(U,mu,sigma));
}

I know that we should have better arithmetic precision
in all this, but I don't think that's the issue here.

I then generate 1000 training points and censor at 0 roughly
the same way as the original poster did.

> N = 1000
> y_raw = rnorm(N,0,1)
> y <- y_raw[y_raw < 0]
> U = 0
> N_observed = length(y)
> N_censored = N - N_observed
> dump(c("N_observed","N_censored","y","U"))

With all this in place, we're getting estimates of
mu around 1.0. And often very ill-behaved adaptation
(massive step sizes).

If I take the line for censored data out of the model,
everything works about how you'd expect (a slight
underestimate of the true mean used to generate the data).
So it seems like it must be a problem with normal_p.

But all the tests for normal_p show all systems go.

I had to add in a signature declaration for normal_p in
stan/gm to get it to run.

- Bob

[Ben Goodrich]

On Tuesday, September 4, 2012 6:10:01 PM UTC-4, Daniel Lee wrote:

On Tue, Sep 4, 2012 at 6:05 PM, Ben Goodrich <goodri...@gmail.com> wrote:
> On Tuesday, September 4, 2012 3:15:09 PM UTC-4, Bob Carpenter wrote:
>>
>> If you only see a few, as you did here, it's not
>> a problem.  We're also working on making that
>> clear by changing "Error" to "Warning" in the report.
>> (What would be really nice is to link them back to
>> a line of the Stan source, but that's MUCH trickier
>> and not on our short-term to-do list.)
>
>
> Can you metaprogram a call to git grep -n $WARNING inside the C++ source?
>

No. git is not on the list of requirements for Stan.

I know. But if we were to metaprogram or makefile it into the source code when make dist happens, then the binaries would be able to output line numbers.

What would that line accomplish? Perhaps we can find another way to do
the same thing.

Rather than metaprogramming, it may be easier to try to do that with shell programming when make dist is called. Like run a script that augments the warning with the line number of the warning.

Ben

[Bob Carpenter]

On 9/4/12 7:07 PM, Ben Goodrich wrote:
> On Tuesday, September 4, 2012 6:10:01 PM UTC-4, Daniel Lee wrote:
>

> On Tue, Sep 4, 2012 at 6:05 PM, Ben Goodrich <goodri...@gmail.com <javascript:>> wrote:
> > On Tuesday, September 4, 2012 3:15:09 PM UTC-4, Bob Carpenter wrote:
> >>
> >> If you only see a few, as you did here, it's not
> >> a problem. We're also working on making that
> >> clear by changing "Error" to "Warning" in the report.
> >> (What would be really nice is to link them back to
> >> a line of the Stan source, but that's MUCH trickier
> >> and not on our short-term to-do list.)
> >
> >
> > Can you metaprogram a call to git grep -n $WARNING inside the C++ source?
> >
>
> No. git is not on the list of requirements for Stan.
>
>
> I know. But if we were to metaprogram or makefile it into the source code when make dist happens, then the binaries
> would be able to output line numbers.

I see.

The problem I see is that I want to point
back to the Stan model code, not the generated
C++ code or underlying C++ code! The text "WARNING"
shows up in the Stan C++ code, not even the code
generated for a model.

From C++, it'd be nice to have something like
a stack trace for where errors happen. But on the C++
side, a simple grep -R will work, too. Most of our
users won't have the git repository installed anyway,
I suspect.

I've never ever downloaded a repository
from a project I wasn't a developer for.

- Bob