cov_matrix_constrain

[Bob Carpenter]

I looked at wishart and recomputed the log
absolute Jacobian determinant for the wishart2x2
case, and my algebra matches the one in the model.

I'm still not getting the same answers from Stan.
The parameters are slightly different in that
wishart uses a covariance parameter and wishart2x2
a correlation parameter.

wishart: cov12,var1,var2 and wishart2x2
wishart2x2: corr12,var1,var2.

I'm so lost with Ben's transform that I don't know
what to do next for debugging.

- Bob

[Ben Goodrich]

This is with or without the Jacobian adjustment that you commented out before?

[Bob Carpenter]

If the log prob isn't being incremented in
cov_matrix_constrain() with the log Jacobian determinant
of the full transform, the steps need to be added back in.
The only transform called out in the read_cov_L() doc is the
correlation -> covariance one and none are called out
in read_corr_L().

I also don't understand why there is a double loop there
of size K when the doc says it's computing a diagonal
Jacobian.

I tried with all four options of commenting out/including
the Jacobian adjustments for the transforms of the first
two parts:

+/- : (-inf,inf) -> (-1,1) [CPCs]
+/- : (-inf,inf) -> (0,inf) [sds]

We get a random init of the parameters in the
unconstrained space. But the parameterizations in
wishart and wishart2x2 are not the same!

In wishart, the unconstrained params are whatever the
LKJ thingy uses. I think that's a CPC (presumably just the
correlation) and two standard deviations?

In the wishart2x2.hpp case, the unconstrained parameters
are exp() of the two variances and atanh() of the correlation.

I'd expect the correlation to have the same
gradient. But maybe not because of the dependency?
I haven't worked through the algebra.

- Bob

[Ben Goodrich]

On Monday, March 12, 2012 9:48:43 PM UTC-4, Bob Carpenter wrote:

If the log prob isn't being incremented in
cov_matrix_constrain() with the log Jacobian determinant
of the full transform, the steps need to be added back in.

I am pretty sure those 2 are supposed to be in.
 

The only transform called out in the read_cov_L() doc is the
correlation -> covariance one and none are called out
in read_corr_L().

I also don't understand why there is a double loop there
of size K when the doc says it's computing a diagonal
Jacobian.

We should cut out the loops since they are equivalent to

log_prob += log_sds.sum() * K * 2.0;

but the Jacobian (sub)matrix is diagonal of order K choose 2. The diagonal has the product of every pair of standard deviations because the the derivative of the i,j covariance with respect to the i,j correlation is the product of the ith and jth standard deviation.

 

I tried with all four options of commenting out/including
the Jacobian adjustments for the transforms of the first
two parts:

   +/- : (-inf,inf) -> (-1,1)    [CPCs]
   +/- : (-inf,inf) -> (0,inf)   [sds]

We get a random init of the parameters in the
unconstrained space.  But the parameterizations in
wishart and wishart2x2 are not the same!

I started to change wishart2x2 so it would have the same parameterization, but I still can't make them line up.

Under wishart.stan with seed = 1331410622

W =
 0.0549377:0   0.155889:0
  0.155889:0    1.57899:0

whereas in my new wishart2x2.stan

W =
  0.0549377:0   0.0843354:0
  0.0843354:0     1.57899:0

so the (squared) standard deviations are right but the correlation is different somehow.

 Ben

[Ben Goodrich]

On Monday, March 12, 2012 10:46:18 PM UTC-4, Ben Goodrich wrote:

We should cut out the loops since they are equivalent to

log_prob += log_sds.sum() * K * 2.0;

That was a mistake. It should be (and now is)

log_prob += (log_sds.sum() + log(2.0)) * K;

The whole thing may well be fixed as of commit 68c0587 . I also reparameterized wishart2x2.stan so that it matches the parameterization of wishart.stan and now they both yield the same output in --test-grad mode. However, the posterior means differ slightly from each other --- albeit more than I would expect from finite precision issues --- when you run them both for 2000 iterations starting with the same seed. So there may be some other bug lurking still.

Ben

[Bob Carpenter]

That's fantastic.

The effective sample size may still be pretty low for
this. I'll take a look at a bunch of longer runs and
see how multiple runs of wishart vary compared to wishart2x2.

- Bob

[Bob Carpenter]

How are you even seeing W? The test_grad function
only prints the unconstrained value. I should go
back and make sure that the first sample is the init
value if there is no warmup. I don't think it is
right now -- I only think I'm printing the first sample
computed via HMC/NUTS.

Should I make the changes or do you want to?

1. replace the two loops in read_cov_L() with:

log_prob += log_sds.sum() * K * 2.0;

But I don't see why that's right. If it's just the
diagonal, then isn't that just transforming var -> sd?
In which case I see the 2.0 * log_sds.sum(), but why
the K?

2. Add the Jacobians back in for the cov_matrix_constrain()

Now I have to go back and follow your unbranch instructions.

- Bob

[Ben Goodrich]

On Monday, March 12, 2012 11:17:38 PM UTC-4, Bob Carpenter wrote:

How are you even seeing W?  The test_grad function
only prints the unconstrained value.  

I added a std::cout << "W = " << std::endl << W << std::endl;

to the C++ file temporarily.
 

Should I make the changes or do you want to?

I did before seeing this post.
 

1.  replace the two loops in read_cov_L() with:

     log_prob += log_sds.sum() * K * 2.0;

But I don't see why that's right.  If it's just the
diagonal, then isn't that just transforming var -> sd?
In which case I see the 2.0 * log_sds.sum(), but why
the K?

It wasn't quite right. The two Jacobians in question here are both diagonal. The first Jacobian is for the transformation from correlations to covariances (of order K choose 2). The derivative of a the i,j covariance with respect to the i,j correlation is sd_i * sd_j, and the product of each K choose 2 pairs of standard deviations goes on the diagonal. So, that log determinant is (K - 1) * log_sds.sum().

The second Jacobian is for the transformation from standard deviations to variances (of order K). The derivative of the ith variance is 2 * sd_i. So, the log determinant is (2.0 * sd).log().sum(), which is K * log(2.0) + sd.log().sum(). So, the correct line (which is in the code now) is

log_prob += (sds.log().sum() + log(2.0)) * K;
 

2.  Add the Jacobians back in for the cov_matrix_constrain()

Yes, that was necessary.
 

Now I have to go back and follow your unbranch instructions.

 
Hopefully, you don't actually have to unmerge a branch that you already merged into master. If you made a branch that you haven't merged into master, then you can just throw the branch away.

Ben

[Bob Carpenter]

Ben:

Could you make the changes?

Indeed, I doin't need to unmerge, just discard,
which is easy. It's just not at my fingertips.

I'm getting sucked into revising the Stan NSF proposal,
so haven't actually coded anything recently or written
any doc.

- Bob

[Ben Goodrich]

On Tuesday, March 13, 2012 12:42:26 PM UTC-4, Bob Carpenter wrote:

Ben:

Could you make the changes?

They have been in since last night. Did you pull this morning?

git pull
git show 882409^1...68c0587

Ben