truncated data works with mean or sd, but not both
41 views
Skip to first unread message
Bob Carpenter
Sep 5, 2012, 2:54:55 AM9/5/12
to stan...@googlegroups.com
I've been exploring some models for truncated data
and am more convinced than ever we need some diagnostics
for what NUTS/HMC is doing.
The simplest model estimates the mean only:
------------
data {
real U;
int<lower=1> N;
real<upper = U> y[N];
}
parameters {
real mu;
}
model {
for (n in 1:N)
y[n] ~ normal(mu,1) T[,0];
}
-------------
That T[,0] restricts the normal to (-inf,0]. There's
no prior on mu.
Around 500 data points were generated in R in the obvious way:
U <- 0
y_raw <- rnorm(1000,0,1)
y <- y_raw[y_raw < 0]
N <- length(y)
It works just fine and all the gradients check out and
the error's about what you'd expect.
If I fix mu at its true value 0 and add a sigma parameter
for sd, everything still works OK, even without a prior
on sigma:
--------------
parameters {
real<lower=0> sigma;
}
model {
for (n in 1:N)
y[n] ~ normal(0,sigma) T[,0];
}
---------------
But if I also treat sigma as a parameter, even with
a proper prior (though one that's uniform on
a small fixed range) on both sigma and mu, things
go wonky.
------------------
parameters {
real<lower=0,upper=2> sigma;
real mu;
}
model {
mu ~ normal(0,1);
for (n in 1:N)
y[n] ~ normal(mu,sigma) T[,0];
}
------------------
I get some warnings about sigma=0 and mu=inf to begin,
then things settle down, then both the
estimates for sigma and mu wind up too large (around 1
for mu instead of 0 and around 1.25
for sigma instead of 1).
With sigma defined on [0,2], its unconstrained zero-init
will be 1, and with mu unconstrained, its zero-init
will be 0.
What's somewhat disturbing is the value of the gradient:
TEST GRADIENT MODE
Log probability=79.0865
param idx value model finite diff error
0 0 -3.00002 -3.00002 -4.88227e-07
1 0 16.9871 16.9871 -4.42688e-08
What I'm wondering is why there's a positive gradient for sigma (the
second parameter)? Things seem to be going the right way in
general:
TEST GRADIENT MODE
Log probability=-462.455
param idx value model finite diff error
0 -1.82834 826.009 826.009 4.62399e-06
1 -0.409461 414.698 414.698 -3.6559e-06
This init has mu way too low, so there's a big positive gradient. But
there's also a big positive gradient on sigma. I assume that's to
account for the bad prediction of mu at that point?
I'm having similar issues with censored models, for what
it's worth.
I just added unit tests for normal_p, and everything checks out
vs. R., so I don't think the problem is there. And the finite
diffs all check out, so it seems unlikely to be the gradients, either.
- Bob
and am more convinced than ever we need some diagnostics
for what NUTS/HMC is doing.
The simplest model estimates the mean only:
------------
data {
real U;
int<lower=1> N;
real<upper = U> y[N];
}
parameters {
real mu;
}
model {
for (n in 1:N)
y[n] ~ normal(mu,1) T[,0];
}
-------------
That T[,0] restricts the normal to (-inf,0]. There's
no prior on mu.
Around 500 data points were generated in R in the obvious way:
U <- 0
y_raw <- rnorm(1000,0,1)
y <- y_raw[y_raw < 0]
N <- length(y)
It works just fine and all the gradients check out and
the error's about what you'd expect.
If I fix mu at its true value 0 and add a sigma parameter
for sd, everything still works OK, even without a prior
on sigma:
--------------
parameters {
real<lower=0> sigma;
}
model {
for (n in 1:N)
y[n] ~ normal(0,sigma) T[,0];
}
---------------
But if I also treat sigma as a parameter, even with
a proper prior (though one that's uniform on
a small fixed range) on both sigma and mu, things
go wonky.
------------------
parameters {
real<lower=0,upper=2> sigma;
real mu;
}
model {
mu ~ normal(0,1);
for (n in 1:N)
y[n] ~ normal(mu,sigma) T[,0];
}
------------------
I get some warnings about sigma=0 and mu=inf to begin,
then things settle down, then both the
estimates for sigma and mu wind up too large (around 1
for mu instead of 0 and around 1.25
for sigma instead of 1).
With sigma defined on [0,2], its unconstrained zero-init
will be 1, and with mu unconstrained, its zero-init
will be 0.
What's somewhat disturbing is the value of the gradient:
TEST GRADIENT MODE
Log probability=79.0865
param idx value model finite diff error
0 0 -3.00002 -3.00002 -4.88227e-07
1 0 16.9871 16.9871 -4.42688e-08
What I'm wondering is why there's a positive gradient for sigma (the
second parameter)? Things seem to be going the right way in
general:
TEST GRADIENT MODE
Log probability=-462.455
param idx value model finite diff error
0 -1.82834 826.009 826.009 4.62399e-06
1 -0.409461 414.698 414.698 -3.6559e-06
This init has mu way too low, so there's a big positive gradient. But
there's also a big positive gradient on sigma. I assume that's to
account for the bad prediction of mu at that point?
I'm having similar issues with censored models, for what
it's worth.
I just added unit tests for normal_p, and everything checks out
vs. R., so I don't think the problem is there. And the finite
diffs all check out, so it seems unlikely to be the gradients, either.
- Bob
Matt Hoffman
Sep 5, 2012, 1:24:14 PM9/5/12
to stan...@googlegroups.com
Odd, I got perfectly reasonable results when I ran this model—samples
for mu in the range of about (-0.5, 1.5) and samples for sigma in the
range of about (0.8, 1.5).
Matt
> --
>
>
for mu in the range of about (-0.5, 1.5) and samples for sigma in the
range of about (0.8, 1.5).
Matt
>
>
Bob Carpenter
Sep 5, 2012, 1:27:09 PM9/5/12
to stan...@googlegroups.com
Thanks!
I've attached the data running mu in (-3,3) and
sigma in (0.1,2) in 0.1 increments. I included
a column for the (unnormalized) log probability. I could plot
the density, but I don't even know what you mean by a vector
field of the gradients, much less how to plot it.
Here's the model on which the numbers are based:
data {
int<lower=1> N;
real<upper=0> y[N];
}
parameters {
real mu;
real<lower=0,upper=2> sigma;
}
N is around 500 and generated from normal(0,1)
and truncated above at 0 (I just generated 1000
and rejected the ones above 0).
- Bob
P.S. I just reinstated your posting priveleges. I think
Daniel turned a bunch off by accident in a move to
cut down on spam.
On 9/5/12 9:04 AM, Michael Betancourt wrote:
> Bob,
>
> Since it's two dimensional, can you make a surface plot of the density
> and a vector field of the gradients? If you scan a grid and save x,
> y, dV/dx, and dV/dy I can throw the plots together quickly, if that's
> more convenient.
>
> -Mike
>
> (I tried posting directly to stan-dev but for some reason it bounced me).
> --
>
I've attached the data running mu in (-3,3) and
sigma in (0.1,2) in 0.1 increments. I included
a column for the (unnormalized) log probability. I could plot
the density, but I don't even know what you mean by a vector
field of the gradients, much less how to plot it.
Here's the model on which the numbers are based:
data {
int<lower=1> N;
real<upper=0> y[N];
}
parameters {
real mu;
real<lower=0,upper=2> sigma;
}
model {
for (n in 1:N)
y[n] ~ normal(mu,sigma) T[,0];
model {
for (n in 1:N)
}
N is around 500 and generated from normal(0,1)
and truncated above at 0 (I just generated 1000
and rejected the ones above 0).
- Bob
P.S. I just reinstated your posting priveleges. I think
Daniel turned a bunch off by accident in a move to
cut down on spam.
On 9/5/12 9:04 AM, Michael Betancourt wrote:
> Bob,
>
> Since it's two dimensional, can you make a surface plot of the density
> and a vector field of the gradients? If you scan a grid and save x,
> y, dV/dx, and dV/dy I can throw the plots together quickly, if that's
> more convenient.
>
> -Mike
>
> (I tried posting directly to stan-dev but for some reason it bounced me).
>
Bob Carpenter
Sep 5, 2012, 2:04:40 PM9/5/12
to stan...@googlegroups.com
I also realized why I was seeing so many rejections.
The simple models I've been looking at have had
median tree depth of one, so rejection rate is roughly
equal to average rejection rate on the second half of
the trees.
Now onto the real issue.
On 9/5/12 1:24 PM, Matt Hoffman wrote:
> Odd, I got perfectly reasonable results when I ran this model�samples
Maybe I messed up the data simulation, but it looks
OK to me plotted in R. I checked the model and data
in at:
src/models/basic_estimators/truncated_normal.stan
/truncated_normal.Rdata
Would you (Matt) mind running the checked in version to see
if it gives you OK results?
Here's what I get after 50K samples in a single
chain (though running longer or more chains
doesn't matter -- it settles on this answer
pretty quickly). Here's what the output looks
like, with simulated mu = 0, sigma = 1:
> quantile(tns$mu)
0% 25% 50% 75% 100%
-0.1016810 0.5404598 0.8620630 1.2935150 3.7038100
> quantile(tns$sigma)
0% 25% 50% 75% 100%
0.961492 1.225290 1.332955 1.457465 1.974080
> colMeans(tns)
lp__ treedepth__ stepsize__ mu sigma
80.4058466 1.3740000 0.1573470 0.9808178 1.3553201
I have switched to the latest compilers on Mountain
Lion. I never used to question the compiler or system
before I moved to C++!
- Bob
The simple models I've been looking at have had
median tree depth of one, so rejection rate is roughly
equal to average rejection rate on the second half of
the trees.
Now onto the real issue.
On 9/5/12 1:24 PM, Matt Hoffman wrote:
> Odd, I got perfectly reasonable results when I ran this model�samples
> for mu in the range of about (-0.5, 1.5) and samples for sigma in the
> range of about (0.8, 1.5).
I like that answer better than mine :-)
> range of about (0.8, 1.5).
Maybe I messed up the data simulation, but it looks
OK to me plotted in R. I checked the model and data
in at:
src/models/basic_estimators/truncated_normal.stan
/truncated_normal.Rdata
Would you (Matt) mind running the checked in version to see
if it gives you OK results?
Here's what I get after 50K samples in a single
chain (though running longer or more chains
doesn't matter -- it settles on this answer
pretty quickly). Here's what the output looks
like, with simulated mu = 0, sigma = 1:
> quantile(tns$mu)
0% 25% 50% 75% 100%
-0.1016810 0.5404598 0.8620630 1.2935150 3.7038100
> quantile(tns$sigma)
0% 25% 50% 75% 100%
0.961492 1.225290 1.332955 1.457465 1.974080
> colMeans(tns)
lp__ treedepth__ stepsize__ mu sigma
80.4058466 1.3740000 0.1573470 0.9808178 1.3553201
I have switched to the latest compilers on Mountain
Lion. I never used to question the compiler or system
before I moved to C++!
- Bob
Matt Hoffman
Sep 5, 2012, 2:08:55 PM9/5/12
to stan...@googlegroups.com
I'm not seeing the checked in model—did you push?
On Wed, Sep 5, 2012 at 11:04 AM, Bob Carpenter <ca...@alias-i.com> wrote:
> I also realized why I was seeing so many rejections.
> The simple models I've been looking at have had
> median tree depth of one, so rejection rate is roughly
> equal to average rejection rate on the second half of
> the trees.
>
> Now onto the real issue.
>
>
> On 9/5/12 1:24 PM, Matt Hoffman wrote:
>>
>> Odd, I got perfectly reasonable results when I ran this model—samples
> --
>
>
On Wed, Sep 5, 2012 at 11:04 AM, Bob Carpenter <ca...@alias-i.com> wrote:
> I also realized why I was seeing so many rejections.
> The simple models I've been looking at have had
> median tree depth of one, so rejection rate is roughly
> equal to average rejection rate on the second half of
> the trees.
>
> Now onto the real issue.
>
>
> On 9/5/12 1:24 PM, Matt Hoffman wrote:
>>
>
>
Bob Carpenter
Sep 5, 2012, 2:30:31 PM9/5/12
to stan...@googlegroups.com
Sorry -- forgot to push. Just pushed it.
Daniel just ran it under Windows and got
the same results as I did.
- Bob
On 9/5/12 2:08 PM, Matt Hoffman wrote:
> I'm not seeing the checked in model�did you push?
Daniel just ran it under Windows and got
the same results as I did.
- Bob
On 9/5/12 2:08 PM, Matt Hoffman wrote:
> I'm not seeing the checked in model�did you push?
>
> On Wed, Sep 5, 2012 at 11:04 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>> I also realized why I was seeing so many rejections.
>> The simple models I've been looking at have had
>> median tree depth of one, so rejection rate is roughly
>> equal to average rejection rate on the second half of
>> the trees.
>>
>> Now onto the real issue.
>>
>>
>> On 9/5/12 1:24 PM, Matt Hoffman wrote:
>>>
>>> Odd, I got perfectly reasonable results when I ran this model�samples
> On Wed, Sep 5, 2012 at 11:04 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>> I also realized why I was seeing so many rejections.
>> The simple models I've been looking at have had
>> median tree depth of one, so rejection rate is roughly
>> equal to average rejection rate on the second half of
>> the trees.
>>
>> Now onto the real issue.
>>
>>
>> On 9/5/12 1:24 PM, Matt Hoffman wrote:
>>>
Matt Hoffman
Sep 5, 2012, 6:15:46 PM9/5/12
to stan...@googlegroups.com
Ok, finally got around to running this. When I run it on your data, I
get the same results you did. When I run it on my data, I get more
reasonable results. The posterior density estimates (courtesy of
ggplot2) are attached.
Could it just be a bad draw? What happens if you sample a new set of data?
Matt
On Wed, Sep 5, 2012 at 11:30 AM, Bob Carpenter <ca...@alias-i.com> wrote:
> Sorry -- forgot to push. Just pushed it.
>
> Daniel just ran it under Windows and got
> the same results as I did.
>
> - Bob
>
>
> On 9/5/12 2:08 PM, Matt Hoffman wrote:
>>
>> I'm not seeing the checked in model—did you push?
> --
>
>
get the same results you did. When I run it on my data, I get more
reasonable results. The posterior density estimates (courtesy of
ggplot2) are attached.
Could it just be a bad draw? What happens if you sample a new set of data?
Matt
On Wed, Sep 5, 2012 at 11:30 AM, Bob Carpenter <ca...@alias-i.com> wrote:
> Sorry -- forgot to push. Just pushed it.
>
> Daniel just ran it under Windows and got
> the same results as I did.
>
> - Bob
>
>
> On 9/5/12 2:08 PM, Matt Hoffman wrote:
>>
>>
>> On Wed, Sep 5, 2012 at 11:04 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>>>
>>> I also realized why I was seeing so many rejections.
>>> The simple models I've been looking at have had
>>> median tree depth of one, so rejection rate is roughly
>>> equal to average rejection rate on the second half of
>>> the trees.
>>>
>>> Now onto the real issue.
>>>
>>>
>>> On 9/5/12 1:24 PM, Matt Hoffman wrote:
>>>>
>>>>
>>>> Odd, I got perfectly reasonable results when I ran this model—samples
>> On Wed, Sep 5, 2012 at 11:04 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>>>
>>> I also realized why I was seeing so many rejections.
>>> The simple models I've been looking at have had
>>> median tree depth of one, so rejection rate is roughly
>>> equal to average rejection rate on the second half of
>>> the trees.
>>>
>>> Now onto the real issue.
>>>
>>>
>>> On 9/5/12 1:24 PM, Matt Hoffman wrote:
>>>>
>>>>
>
>
Daniel Lee
Sep 5, 2012, 6:40:28 PM9/5/12
to stan...@googlegroups.com
Can you check your data into the repository?
> --
>
>
>
>
Matt Hoffman
Sep 5, 2012, 6:46:34 PM9/5/12
to stan...@googlegroups.com
Sure. Pushed.
Matt
> --
>
>
Matt
>
>
Bob Carpenter
Sep 5, 2012, 7:43:08 PM9/5/12
to stan...@googlegroups.com
But doesn't that mean (sorry again for
all the puns) that the posterior mean of
mu is above 0.5 and the posterior mean for
sigma is around 1.3? That's close to what
I got:
quantile(tns$mu)
0% 25% 50% 75% 100%
-0.1016810 0.5404598 0.8620630 1.2935150 3.7038100
quantile(tns$sigma)
0% 25% 50% 75% 100%
0.961492 1.225290 1.332955 1.457465 1.974080
The "right" answer is mu=0 and sigma=1,
which don't get much mass in the posterior!
This is another good example where the MAP estimate
(which is around mu=0 and sigma=1 according to our
model) is pretty far from the posterior mean on the
scale we care about.
- Bob
all the puns) that the posterior mean of
mu is above 0.5 and the posterior mean for
sigma is around 1.3? That's close to what
I got:
quantile(tns$mu)
0% 25% 50% 75% 100%
-0.1016810 0.5404598 0.8620630 1.2935150 3.7038100
quantile(tns$sigma)
0% 25% 50% 75% 100%
0.961492 1.225290 1.332955 1.457465 1.974080
which don't get much mass in the posterior!
This is another good example where the MAP estimate
(which is around mu=0 and sigma=1 according to our
model) is pretty far from the posterior mean on the
scale we care about.
- Bob
Matt Hoffman
Sep 5, 2012, 9:23:01 PM9/5/12
to stan...@googlegroups.com
The right answer doesn't get any mass under the posterior!
Pedantry aside, the means I got aren't too unreasonable: they're .08
and 1.07 for mu and sigma. The ranges are a little misleading because
the distribution is skewed. Here are quantiles:
> quantile(temp5$mu)
> quantile(temp5$sigma)
Matt
> --
>
>
Pedantry aside, the means I got aren't too unreasonable: they're .08
and 1.07 for mu and sigma. The ranges are a little misleading because
the distribution is skewed. Here are quantiles:
> quantile(temp5$mu)
0% 25% 50% 75% 100%
-0.5246780 -0.1060650 0.0391880 0.2226787 2.6395300
> quantile(temp5$sigma)
0% 25% 50% 75% 100%
0.817321 0.996989 1.060570 1.133740 1.785260
Matt
>
>
Bob Carpenter
Sep 6, 2012, 1:30:29 PM9/6/12
to stan...@googlegroups.com
On 9/5/12 9:23 PM, Matt Hoffman wrote:
> The right answer doesn't get any mass under the posterior!
> Pedantry aside,
Never! I want to get this right.
> the means I got aren't too unreasonable: they're .08
> and 1.07 for mu and sigma.
density plot you sent.
Anyway, I just wanted to work this out to the point
where we were confident it wasn't a bug in the
way truncated densities were calculated.
- Bob
Matt Hoffman
Sep 6, 2012, 1:55:03 PM9/6/12
to stan...@googlegroups.com
> That's very surprising given the posterior
> density plot you sent.
I meant the means for the dataset that I generated
> density plot you sent.
(normal_truncated2.Rdata)—I get less pleasing results for the original
data (normal_truncated.Rdata).
Matt
Nigel/essence
Sep 9, 2012, 7:29:27 AM9/9/12
to stan...@googlegroups.com, mdho...@cs.princeton.edu
Are you saying there is no junk in the trunk?
Reply all
Reply to author
Forward
0 new messages