Re: 8schools!?!?!?!? + epsilon/mass + epsilon*L = 1

25 views
Skip to first unread message

Bob Carpenter

unread,
Jul 19, 2012, 3:03:36 PM7/19/12
to stan...@googlegroups.com


On 7/19/12 2:44 PM, Andrew Gelman wrote:
> B
>
> I'll try some sims and see what happens. I am happy to have separate mass matrix and epsilon.
> I actually find it helpful to separate these. As Mike explained to me, we can roughly think of the
> mass matrix as corresponding to the variance of the target distribution (that's not really right
> in general but it makes sense in the simple case) and then it makes sense to have epsilon*L
> to be approx 1.

Just the separate mass matrix and epsilon are not
identified. So you can always reascale epsilon
and the mass matrix so that epsilon * L = 1 if you're
doing HMC. For NUTS, L isn't fixed, plus you sample
along the path, not just at the end.

I'd like to know what's special about 1 here. NUTS II tries
to find an epsilon that works best for ALL params in terms
of achieving a 65% accept rate, then it fits a diagonal
mass matrix around that. So I don't see why you'd wind
up with epsilon * L = 1 as a result.

My understanding is that if epsilon * L = 1, the
exact balance between epsilon and L involves
balancing so that L is small and epsilon is large for
efficient movement with epsilon being small enough for
accurate Hamiltonian simulation.

> That's what you're getting: you have epsilon*L as about 6, which isn't
> far from the scale of most of the parameters in the model.
>
> Regarding the point at the end of your email, I'd be even happier if, instead of restarting a chain, we continued it at its current value with better tuning parameters. If that doesn't work, we can restart. But changing tuning parameters without changing the position would be cleaner.

We would be, too. We only do this when we
can't continue the current value because it's NaN or
+/- infinity. Once that's happen, there's no way to
get out of the hole through addition and multiplication.

I suppose we could go back to the initial values
and try again instead of regenerating those, too, but
often the problem is a poor (very low probability,
highly curved space) choice of starting point.

- Bob

Michael Betancourt

unread,
Jul 19, 2012, 3:16:01 PM7/19/12
to stan...@googlegroups.com
> Just the separate mass matrix and epsilon are not
> identified. So you can always reascale epsilon
> and the mass matrix so that epsilon * L = 1 if you're
> doing HMC. For NUTS, L isn't fixed, plus you sample
> along the path, not just at the end.

I would argue that the step size and the mass matrix are NOT
non-identified, even in regular HMC. In the limit \epilson -> 0 that
might be true, but when \epsilon is non-zero you have to worry about
the momenta update,
a la p += -\epsilon * gradV (notice the mass matrix doesn't show up
here). If \epsilon is too large then the update can be unstable, and
there's no mass matrix small enough to suppress that instability (you
know, other than zero).

(Matt and I have discussed this before. He's in Bob's camp, but I'm
pretty sure the finite step size case is not the same as the ideal
case).

> I'd like to know what's special about 1 here. NUTS II tries
> to find an epsilon that works best for ALL params in terms
> of achieving a 65% accept rate, then it fits a diagonal
> mass matrix around that. So I don't see why you'd wind
> up with epsilon * L = 1 as a result.

If you've successfully standardized your distribution then it should
look like a standardized gaussian, where the width of the typical set
is about 1. In that case the optimal distance to travel is 1 =
\epsilon * L.

Bob Carpenter

unread,
Jul 19, 2012, 6:51:28 PM7/19/12
to stan...@googlegroups.com
On 7/19/12 4:15 PM, Andrew Gelman wrote:>

> What is your avg acceptance rate using this hmc? I ran
> it using all your settings including U(-2,2) starting
> points and I'm getting estimates of sigma.theta seem too low,
> not-so-great convergence, and an acceptance rate of 99% or 100%.
> So maybe I have a bug?

I found near 100% acceptance rats for both NUTS and
the fixed-parameter HMC runs.

> Also, in your summaries below, are you using the
> new n.eff definition and split R-hat?

The summaries are just what comes out of R2WinBUGs.
So it's whatever Andrew and his successors did there for n.eff
and plain old R-hat. (I know we need to get split R-hat
and whatever Kenny and Andrew decide is a good way to calculate
n.eff into the output analysis -- Jiqiang/Daniel already did this
for the Rcpp-based RStan.)

Here's a new run (so I could calc acceptance rates).

HMC, epsilon = 0.05, leapfrog_steps=128:

8 chains, each with 20000 iterations (first 10000 discarded), n.thin = 10
n.sims = 8000 iterations saved
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
lp__ -16.6 7.3 -27.7 -21.6 -17.9 -12.9 1.6 1 1200
theta[1] 11.4 8.3 -1.8 6.0 10.2 15.6 31.2 1 4500
theta[2] 7.8 6.2 -4.5 3.9 7.7 11.6 20.5 1 6600
theta[3] 6.1 7.8 -11.4 2.0 6.6 10.9 20.5 1 6900
theta[4] 7.7 6.6 -5.2 3.7 7.6 11.6 21.0 1 8000
theta[5] 5.0 6.4 -9.1 1.2 5.5 9.3 16.5 1 5200
theta[6] 6.1 6.7 -8.3 2.2 6.5 10.4 18.6 1 3100
theta[7] 10.7 6.9 -1.3 6.1 10.0 14.6 26.0 1 6900
theta[8] 8.3 7.9 -7.6 3.8 8.1 12.6 24.8 1 8000
mu_theta 7.8 5.1 -2.2 4.6 7.8 11.1 17.8 1 3800
sigma_theta 6.7 5.7 0.3 2.6 5.3 9.2 20.9 1 1200

NUTS II

8 chains, each with 20000 iterations (first 10000 discarded), n.thin = 10
n.sims = 8000 iterations saved
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
lp__ -19.0 4.7 -28.0 -22.2 -19.2 -15.8 -9.1 1 130
treedepth__ 2.3 1.0 1.0 2.0 3.0 3.0 4.0 1 210
stepsize__ 4.1 0.1 3.9 4.1 4.1 4.3 4.3 1 1
theta[1] 12.4 8.6 -2.4 7.4 11.4 16.9 32.6 1 1700
theta[2] 8.2 6.6 -5.4 4.4 8.3 12.2 21.5 1 2000
theta[3] 6.0 8.1 -12.4 1.7 7.1 10.7 20.8 1 1600
theta[4] 7.8 6.8 -6.1 3.9 8.0 11.8 21.3 1 810
theta[5] 4.7 6.7 -9.8 0.6 5.6 9.6 16.3 1 560
theta[6] 6.0 7.1 -9.4 1.6 6.3 10.5 19.0 1 810
theta[7] 11.7 6.8 -0.8 7.6 11.2 15.6 26.5 1 2300
theta[8] 8.8 8.2 -7.6 4.4 8.4 13.6 26.8 1 1400
mu_theta 8.3 5.3 -2.2 5.2 8.4 11.5 18.6 1 850
sigma_theta 7.8 5.3 2.2 3.9 6.4 10.2 21.5 1 160


sigma_theta and lp__ seem way off, and though the posterior
sd is of the magnitude of the parameters, sqrt(n.eff) is high
enough I still wouldn't expect to see this much variation (unless,
of course, n.eff is itself way off!).

And it seems systematic (I did several runs and each HMC run
looked similar and each NUTS run looked similar).
For instance, if I initialize both to
start at 0 on the unconstrained scale in each chain (instead
of each param on each chain being drawn uniformly in [-2,2]).
Here's the output:

HMC(epsilon=0.05,L=128)

8 chains, each with 20000 iterations (first 10000 discarded), n.thin = 10
n.sims = 8000 iterations saved
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
lp__ -16.7 7.1 -27.7 -21.6 -18.0 -12.9 0.7 1 4600
theta[1] 11.3 8.3 -2.3 6.0 10.3 15.5 31.4 1 8000
theta[2] 7.9 6.3 -4.9 3.9 7.9 11.8 20.7 1 6000
theta[3] 6.2 7.8 -11.4 2.1 6.7 10.9 20.7 1 4500
theta[4] 7.8 6.7 -5.9 3.7 7.7 11.9 21.2 1 5700
theta[5] 5.2 6.4 -8.9 1.4 5.6 9.4 16.7 1 8000
theta[6] 6.1 6.8 -8.4 2.1 6.5 10.5 18.9 1 6900
theta[7] 10.7 6.7 -1.2 6.1 10.1 14.6 25.5 1 8000
theta[8] 8.4 7.8 -6.9 3.8 8.2 12.6 25.5 1 8000
mu_theta 7.9 5.2 -2.3 4.7 7.9 11.1 18.4 1 8000
sigma_theta 6.6 5.6 0.4 2.5 5.3 9.1 20.8 1 4400


NUTS

8 chains, each with 20000 iterations (first 10000 discarded), n.thin = 10
n.sims = 8000 iterations saved
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
lp__ -19.5 4.4 -28.1 -22.5 -19.5 -16.3 -11.3 1 740
treedepth__ 2.4 0.9 1.0 2.0 3.0 3.0 4.0 1 1200
stepsize__ 4.0 0.1 3.9 4.0 4.0 4.1 4.3 1 1
theta[1] 12.5 9.0 -2.6 6.5 11.6 17.2 33.4 1 1600
theta[2] 8.0 6.7 -5.2 3.7 7.9 12.1 21.8 1 2600
theta[3] 5.7 8.4 -12.9 0.9 6.2 11.0 21.3 1 2100
theta[4] 7.7 6.9 -6.1 3.2 7.7 12.2 21.5 1 2000
theta[5] 4.5 6.7 -9.9 0.3 4.8 9.2 17.0 1 1100
theta[6] 5.8 7.3 -9.7 1.4 6.2 10.7 18.9 1 6400
theta[7] 11.5 7.0 -1.3 6.7 10.9 15.6 27.0 1 1800
theta[8] 8.8 8.5 -7.5 3.6 8.7 13.6 27.6 1 1200
mu_theta 8.1 5.3 -2.0 4.6 8.1 11.5 18.6 1 1300
sigma_theta 8.1 5.4 2.2 4.2 6.7 10.4 21.6 1 740

(And don't be misled by the n.eff here -- I ran each for 20K iterations,
not for the same amount of time -- NUTS is WAY faster as it's doing
many fewer leapfrog steps per iteration than HMC at these settings.)

- Bob

Bob Carpenter

unread,
Jul 19, 2012, 7:17:02 PM7/19/12
to stan...@googlegroups.com


On 7/19/12 3:16 PM, Michael Betancourt wrote:
>> Just the separate mass matrix and epsilon are not
>> identified. So you can always reascale epsilon
>> and the mass matrix so that epsilon * L = 1 if you're
>> doing HMC. For NUTS, L isn't fixed, plus you sample
>> along the path, not just at the end.
>
> I would argue that the step size and the mass matrix are NOT
> non-identified, even in regular HMC. In the limit \epilson -> 0 that
> might be true, but when \epsilon is non-zero you have to worry about
> the momenta update,
> a la p += -\epsilon * gradV (notice the mass matrix doesn't show up
> here). If \epsilon is too large then the update can be unstable, and
> there's no mass matrix small enough to suppress that instability (you
> know, other than zero).
>
> (Matt and I have discussed this before. He's in Bob's camp, but I'm
> pretty sure the finite step size case is not the same as the ideal
> case).

I didn't even know there were camps! I thought
this was just algebra. But now that I re-read
section 5.4.2.4 in Radford's handbook paper,
where he concludes they're equivalent,
I see that he modifies the epsilons on the
p updates so that they are parameter dependent:

p[i] += - epsilon[i] * gradV[i]

He then concludes that the resulting position/momentum
values "can no longer be interpreted as following Hamiltonian
dynamics at consistent time points" and goes on to
say "but that is of no consequence for the use of these
trajectories in HMC".

From the above, Michael's talking about the natural
mass-based physical update, which takes

p[i] += epsilon * gradV[i]

I'm guessing Matt was talking about Radford's, so you
were considering different definitions of varying step
size.

Now I'll have to dive into our code to see what NUTS is
actually doing, having not recognized this distinction earlier.

> If you've successfully standardized your distribution then it should
> look like a standardized gaussian, where the width of the typical set
> is about 1. In that case the optimal distance to travel is 1 =
> \epsilon * L.

But isn't epsilon * L equal to time, not distance traveled?
With fixed epsilon and L, distance traveled has to depend on
the initial position and initial random momentum.

- Bob

Michael Betancourt

unread,
Jul 19, 2012, 8:07:52 PM7/19/12
to stan...@googlegroups.com, stan...@googlegroups.com
> I didn't even know there were camps! I thought
> this was just algebra. But now that I re-read
> section 5.4.2.4 in Radford's handbook paper,
> where he concludes they're equivalent,
> I see that he modifies the epsilons on the
> p updates so that they are parameter dependent:
>
> p[i] += - epsilon[i] * gradV[i]
>
> He then concludes that the resulting position/momentum
> values "can no longer be interpreted as following Hamiltonian
> dynamics at consistent time points" and goes on to
> say "but that is of no consequence for the use of these
> trajectories in HMC".
>
> From the above, Michael's talking about the natural
> mass-based physical update, which takes
>
> p[i] += epsilon * gradV[i]
>
> I'm guessing Matt was talking about Radford's, so you
> were considering different definitions of varying step
> size.
>
> Now I'll have to dive into our code to see what NUTS is
> actually doing, having not recognized this distinction earlier.

Extreme mass matrices lead to large momenta which must be suppressed by the multiplication by the mass matrix upon the position updates. If we had infinite precision computers then the two approaches would give the same results, but on finite precision machines you can encounter adding/subtracting large numbers. I think the bigger issue, however, is that the large momenta updates can pass the stability threshold and cause the trajectory to diverge.

The mathematical question is whether or not the stability threshold scales with the mass matrix or not. That's probably a better definition of our camps. :-)

> But isn't epsilon * L equal to time, not distance traveled?
> With fixed epsilon and L, distance traveled has to depend on
> the initial position and initial random momentum.

Ah, but remember that the average magnitude of the momentum is the square root of the dimension, which is order unity, and 1*1=1. In RMHMC the idea is that the more complicated momenta distribution helps to compensate so that the total distance traveled is just right for any distribution.

Matt Hoffman

unread,
Jul 19, 2012, 8:08:27 PM7/19/12
to stan...@googlegroups.com
NUTS is using the Radford-style parameter-specific step sizes, which
produces discretized dynamics for the positions (and the energy of the
momenta) that are equivalent to using a diagonal mass matrix with
appropriate variances, which is equivalent to transforming the space
that the variables live in to (hopefully) look more like a standard
normal.

I am pretty much 100% confident that the overall scale of the inverse
mass matrix and the squared step size are completely interchangeable
in the sense that scaling the inverse mass matrix by c and the step
size by 1/sqrt(c) will produce no change in the discretized dynamics
of the position vector or the energy of the momentum vector. (The
momenta will be different, but they'll have the same energy.)

If I may be pedantic, the width of the typical set of a standard
normal is more like 4 or 6 than 1. To continue to be (slightly less)
pedantic, if you've really transformed things so that you're working
with a standard normal, then I'm not sure that you really want to
always travel for epsilon*L=size of typical set, since you may hit a
weird resonance.

Matt

Andrew Gelman

unread,
Jul 20, 2012, 10:54:36 AM7/20/12
to stan...@googlegroups.com
B

Hmmm, if your acc rates were close to 100%, maybe I don't have a bug at all? I'm getting acceptance rates of 99-100% myself:

4 chains, keeping last 5000 iterations out of 10000
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
theta[1] 11.4 8.1 -2.0 6.1 10.3 15.8 30.5 1 1700
theta[2] 7.9 6.1 -4.3 4.1 7.9 11.7 20.6 1 5300
theta[3] 6.2 7.8 -11.6 2.2 6.7 10.9 20.2 1 620
theta[4] 7.7 6.3 -4.9 3.8 7.6 11.6 20.6 1 2300
theta[5] 5.2 6.4 -9.1 1.5 5.7 9.3 16.5 1 1500
theta[6] 6.3 6.6 -8.0 2.5 6.6 10.4 18.7 1 850
theta[7] 10.6 6.6 -1.0 6.2 10.1 14.5 25.5 1 2200
theta[8] 8.5 7.8 -6.5 4.2 8.4 12.7 25.2 1 1400
mu 8.0 5.0 -1.8 4.8 7.9 11.0 17.8 1 1600
log (tau) 1.5 1.0 -0.9 1.0 1.7 2.2 3.0 1 220
Avg acceptance probs: 0.99 0.99 0.99 1

exp(1.7) = 5.5, so the inference for tau is basically the same as you got.

So the bad (or maybe I should say, good) news is that non-adaptive HMC really sucks for this problem. Maybe I'll play around with randomizing epsilon and L at each step. Or perhaps I should just fit a different model than 8-schools for the demo in my book. In the book, I'm happy to be honest about the method's failings but I'd like the demo to work cleanly.

P.S. here's the silly function I wrote to compute split R-hat. I don't have the new n_eff, though; I'm waiting on all of you to create that!

monitor.split <- function (a){
dims <- dim(a)
n <- dims[1]
m <- dims[2]
a1 <- a[1:(n/2),,]
a2 <- a[(n/2+1):n,,]
dims.fat <- c(n/2,2*m,dims[3])
a.fat <- array (NA, dims.fat)
dimnames (a.fat) <- dimnames(a)
a.fat[,1:m,] <- a1
a.fat[,(m+1):(2*m),] <- a2
return (monitor (a.fat))
}

P.P.S. I reread Radford's chapter and it seems to me that when he talks about optimal efficiency, he doesn't at all take into account that computation costs are proportional to L. I'd think this could make a big difference in the optimal acceptance rate even in theory.

Matt Hoffman

unread,
Jul 20, 2012, 1:48:00 PM7/20/12
to stan...@googlegroups.com
RE:P.P.S.: I think Radford's assuming that you're holding εL fixed,
and varying ε and L as necessary to ensure a reasonable acceptance
rate. So cost is proportional to 1/ε, and 1/ε has to grow as d^1/4, so
cost grows as d^5/4.

Bob Carpenter

unread,
Jul 20, 2012, 2:39:28 PM7/20/12
to stan...@googlegroups.com


On 7/20/12 1:48 PM, Matt Hoffman wrote:
> RE:P.P.S.: I think Radford's assuming that you're holding οΏ½L fixed,
> and varying οΏ½ and L as necessary to ensure a reasonable acceptance
> rate. So cost is proportional to 1/οΏ½, and 1/οΏ½ has to grow as d^1/4, so
> cost grows as d^5/4.

I seem to recall two contexts. One is to lower epsilon to get
a reasonable acceptance rate in all parts of the density and the
other is to avoid periodicities.

It's so cool that NUTS just does all this for us!

We're within epsilon of the BUGS estimates on all of
their models with NUTS (other than one parameter of one model, where
we're off by like 1/10th of a percent from the BUGS estimate -- Daniel's
still worried about that one and still looking into it).
We clearly never would've gotten that close with plain old HMC.
What's worrying is that it seems to be biased, not just noisy.

- Bob

Matt Hoffman

unread,
Jul 20, 2012, 2:53:52 PM7/20/12
to stan...@googlegroups.com
"Biased" is such a strong word. If we let it run for long enough with
a small enough randomized step size I'm sure HMC'd figure things out.

On Fri, Jul 20, 2012 at 11:39 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>
>
> On 7/20/12 1:48 PM, Matt Hoffman wrote:
>>
>> RE:P.P.S.: I think Radford's assuming that you're holding εL fixed,
>> and varying ε and L as necessary to ensure a reasonable acceptance
>> rate. So cost is proportional to 1/ε, and 1/ε has to grow as d^1/4, so

Andrew Gelman

unread,
Jul 20, 2012, 3:47:49 PM7/20/12
to stan...@googlegroups.com
I see that working in the private sector hasn't dulled your satirical edge.

Jiqiang Guo

unread,
Jul 20, 2012, 3:56:15 PM7/20/12
to stan...@googlegroups.com
I could not remember which one I am having problem with matching the results to BUGS or JAGS. 
But clearly JAGS has problem in  term of estimation with some models, say vol1/blocker, or vol1/kidney.  
Probably BUGS has the problem as well and it just presented the results from one run.  Or it is just that
the model itself is hard to estimate. 

--
Jiqiang 

Andrew Gelman

unread,
Jul 20, 2012, 4:00:11 PM7/20/12
to stan...@googlegroups.com
For Stan 1.0 we can make life easier by comparing only to Jags and using Jags starting points.  If we have any models where we're confused about which program is giving the right answer (because we seem to be converging to two different distributions), we can run for a zillion iterations using starting points from both estimates.
A

Jiqiang Guo

unread,
Jul 20, 2012, 4:07:48 PM7/20/12
to stan...@googlegroups.com, mdho...@cs.princeton.edu
Just to report that for vol1/dogs example, when 
--seed=1340384924 --chain_id=4,
currently, we ended up at one point.  Not sure what went wrong though. 

Dump of the samples.csv with line NO. 

   1 # Samples Generated by Stan
   2 #
   3 # stan_version_major=1
   4 # stan_version_minor=0
   5 # stan_version_patch=0
   6 # data=dogs.Rdata
   7 # init=random initialization
   8 # append_samples=0
   9 # seed=1340384924
  10 # chain_id=4
  11 # iter=2000
  12 # warmup=1000
  13 # thin=1
  14 # unit_mass_matrix=0
  15 # leapfrog_steps=-1
  16 # max_treedepth=10
  17 # epsilon=-1
  18 # epsilon_pm=0
  19 # delta=0.5
  20 # gamma=0.05
  21 #
  22 lp__,treedepth__,stepsize__,alpha,beta,A,B
  23 # (mcmc::nuts_diag) adaptation finished
  24 # step size=3.9801e+35
  25 # parameter masses:
  26 # 1.99281e-66,1.41421
  27 6.39916e+81,5,3.9801e+35,-0.322508,-1e-05,0.72433,0.99999
  28 6.39916e+81,4,3.9801e+35,-0.322508,-1e-05,0.72433,0.99999
  29 6.39916e+81,4,3.9801e+35,-0.322508,-1e-05,0.72433,0.99999
  30 6.39916e+81,5,3.9801e+35,-0.322508,-1e-05,0.72433,0.99999
  31 6.39916e+81,3,3.9801e+35,-0.322508,-1e-05,0.72433,0.99999
.......
1022 6.39916e+81,4,3.9801e+35,-0.322508,-1e-05,0.72433,0.99999
1023 6.39916e+81,1,3.9801e+35,-0.322508,-1e-05,0.72433,0.99999
1024 6.39916e+81,1,3.9801e+35,-0.322508,-1e-05,0.72433,0.99999
1025 6.39916e+81,2,3.9801e+35,-0.322508,-1e-05,0.72433,0.99999
1026 6.39916e+81,2,3.9801e+35,-0.322508,-1e-05,0.72433,0.99999


--
Jiqiang

Andrew Gelman

unread,
Jul 20, 2012, 4:12:37 PM7/20/12
to stan...@googlegroups.com, mdho...@cs.princeton.edu
1e35, huh?  That's the kind of thing I see in my simulations too.  Isn't Nuts supposed to stop this from happening>>>

Matt Hoffman

unread,
Jul 20, 2012, 4:12:21 PM7/20/12
to stan...@googlegroups.com
This smells like a modeling bug to me. How else could we be getting
log joints of 6.4e81?

In particular, is this line right?

y[dog, trial] ~ bernoulli(exp(alpha * xa[dog, trial] + beta *
xs[dog, trial]));

It seems like there's nothing that prevents that probability from
being greater than 1.

Matt

Andrew Gelman

unread,
Jul 20, 2012, 4:16:11 PM7/20/12
to stan...@googlegroups.com
What happens when we have Bernoulli with prob > 1? Does it just set the density to 0?
Regarding the dogs example: I actually know all about this one, as it's in my book (assuming it's the same dogs example that I know about). The parameters are supposed to be restricted so that both alpha and beta are negative. Also, the model is set up so that xa and xs are always nonnegative. If you actually send in negative values for xa or xs, or positive values for alpha or beta, than I have no idea what will happen.
A

Matt Hoffman

unread,
Jul 20, 2012, 4:18:22 PM7/20/12
to stan...@googlegroups.com
Ah, I missed the constraint. And I assume that if we'd been getting
probabilities greater than 1 then bernoulli_log() would have thrown a
domain error. I'll look at this in a little more detail.

Matt

Andrew Gelman

unread,
Jul 20, 2012, 4:33:52 PM7/20/12
to stan...@googlegroups.com
It's one of my favorite examples. I had no idea it made its way into Bugs. It's from an experiment from about 1950 and the data were originally analyzed by Fred Mosteller, who I was a TA for in his last semester teaching at Harvard. It was also an early example of simulation-based model checking (from 1954), Mosteller compared real dogs to "stat-dogs."

Matt Hoffman

unread,
Jul 20, 2012, 5:15:08 PM7/20/12
to stan...@googlegroups.com
Fixed it and pushed. There was a bug in ub_constrain() where it was
subtracting instead of adding, which meant that the model got infinite
points for setting beta=-1e-5. Seems to work fine now--the estimates
agree with the published BUGS estimates.

Perhaps this also explains the 1/10th% problem that Daniel was looking into?

Matt

Bob Carpenter

unread,
Jul 20, 2012, 8:15:43 PM7/20/12
to stan...@googlegroups.com


On 7/20/12 4:12 PM, Matt Hoffman wrote:
> This smells like a modeling bug to me. How else could we be getting
> log joints of 6.4e81?
>
> In particular, is this line right?
>
> y[dog, trial] ~ bernoulli(exp(alpha * xa[dog, trial] + beta *
> xs[dog, trial]));
>
> It seems like there's nothing that prevents that probability from
> being greater than 1.

That's funny -- I said exactly the same thing when I first
looked at the model. The answer is that each term in the sum is constrained
to be negative because alpha and beta are constrained < 0 and the
xa and xs are all positive.

The Bernoulli would also gripe and return NaN if the chance of
success param is outside of (0,1).

There's also the issue that the mass for alpha adapted to 2e-66
and the step size adapted to 4e+35! Something's definitely
off in the model. I'll take another side-by-side squint at it.

- Bob

Bob Carpenter

unread,
Jul 20, 2012, 8:24:36 PM7/20/12
to stan...@googlegroups.com
Damn those signs. This was almost certainly my
fault in the original ub_constrain() code.

- Bob
>>>>>>>>> RE:P.P.S.: I think Radford's assuming that you're holding οΏ½L fixed,
>>>>>>>>> and varying οΏ½ and L as necessary to ensure a reasonable acceptance
>>>>>>>>> rate. So cost is proportional to 1/οΏ½, and 1/οΏ½ has to grow as d^1/4, so

Bob Carpenter

unread,
Jul 20, 2012, 8:47:42 PM7/20/12
to stan...@googlegroups.com
I just pushed a fix for the unit test, which
matched the buggy definition of ub (the danger
of the same person using the same definition for
the function and the unit test).

- Bob


On 7/20/12 5:15 PM, Matt Hoffman wrote:

Bob Carpenter

unread,
Jul 24, 2012, 4:07:59 AM7/24/12
to stan...@googlegroups.com


On 7/19/12 8:08 PM, Matt Hoffman wrote:
> NUTS is using the Radford-style parameter-specific step sizes,

In light of this, I changed the command-line specification
"unit_mass_matrix" to "equal_step_sizes" and I changed the
code variables and doc to match.

Thanks for the clarification, Matt.

- Bob
Reply all
Reply to author
Forward
0 new messages