fvar distributions!!!

[Daniel Lee]

I think I've got it working using our old code. More to come. We need to do some more testing. Peter, want to help me out? The test distribution tests should test fvars too.

Also, I think it's not efficient and it probably could be made more efficient.

Daniel

[Peter Li]

That's great news. I'm not sure how exactly the distribution testing
framework works, but I'd be glad to help / will look into it.


Peter

> --
> You received this message because you are subscribed to the Google
> Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Ben Goodrich]

This sounds great. The other day I was thinking that we should have a test somewhere that, for each distribution,

-- specified its (f)var parameters
-- drew a large random sample from that distribution given those parameters
-- evaluated the log-likelihood given the true parameters and verified that the average score vector was sufficiently close to a zero vector

That would go a long way toward verifying that everything for that distribution is coherent.

Ben

[Daniel Lee]

I've pushed my changes to feature/fvar_operands_and_partials. I've added a single test for normal at src/test/agrad/fwd/prob/univariate/continuous/normal_test.cpp.

I think the right thing to do is to extend the testing framework to generate every valid combination of tests. For normal_log, that means each of the 3 arguments can be:

1. double

2. var

3. fvar<double>

4. vector<double>

5. vector<var>

6. vector<fvar<double> >

7. Eigen::Matrix<double, 1, -1> (row_vector)
8. Eigen::Matrix<var, 1, -1> (row_vector)
9. Eigen::Matrix<fvar<double>, 1, -1> (row_vector)
10. Eigen::Matrix<double, -1, 1> (vector)
11. Eigen::Matrix<var, -1, 1> (vector)
12. Eigen::Matrix<fvar<double>, -1, 1> (vector)

Except, we won't have all 12^3 combinations. If I understood the induction arguments correctly, we can't have fvar<double> and var in the same argument list.

Should we add fvar<var> in there too? That's going to get complicated.

On Sunday, June 2, 2013 12:31:53 AM UTC-4, Peter Li wrote:

That's great news. I'm not sure how exactly the distribution testing
framework works, but I'd be glad to help / will look into it.


Peter

On 6/2/13 12:18 AM, Daniel Lee wrote:
> I think I've got it working using our old code. More to come. We need
> to do some more testing. Peter, want to help me out? The test
> distribution tests should test fvars too.
>
> Also, I think it's not efficient and it probably could be made more
> efficient.
>
>
> Daniel

[Bob Carpenter]

On 6/4/13 12:14 PM, Daniel Lee wrote:
> I've pushed my changes to feature/fvar_operands_and_partials. I've added a single test for normal at
> src/test/agrad/fwd/prob/univariate/continuous/normal_test.cpp.
>
> I think the right thing to do is to extend the testing framework to generate every valid combination of tests. For
> normal_log, that means each of the 3 arguments can be:
> 1. double
> 2. var
> 3. fvar<double>
> 4. vector<double>
> 5. vector<var>
> 6. vector<fvar<double> >
> 7. Eigen::Matrix<double, 1, -1> (row_vector)
> 8. Eigen::Matrix<var, 1, -1> (row_vector)
> 9. Eigen::Matrix<fvar<double>, 1, -1> (row_vector)
> 10. Eigen::Matrix<double, -1, 1> (vector)
> 11. Eigen::Matrix<var, -1, 1> (vector)
> 12. Eigen::Matrix<fvar<double>, -1, 1> (vector)
>
> Except, we won't have all 12^3 combinations. If I understood the induction arguments correctly, we can't have
> fvar<double> and var in the same argument list.
>
> Should we add fvar<var> in there too? That's going to get complicated.

Yes, and fvar<fvar<var> > and Matrix<fvar<fvar<var> > > and so on.

I think we'll just have to trust that the autodiff
is working because we unit test that independently.

For var, we've customized the autodiff in the functions, so
I think it's fair to test that. I think fvar<double> would
suffice.

The rule is that if there an autodiff type, all the autodiff
types are the same. So I think we can get away with just
testing fvar<double> in addition to var, which given the rule,
should only double the total number of tests if you want to
do it exhaustively.

We could also somehow mock out the vars and fvars using yet
another OperandsAndPartials specialization. That would modularize
the testing for getting the right partials and then dumping out
the right fvar or var.

- Bob

[Daniel Lee]

On Sun, Jun 2, 2013 at 1:41 AM, Ben Goodrich <goodri...@gmail.com> wrote:

This sounds great. The other day I was thinking that we should have a test somewhere that, for each distribution,

-- specified its (f)var parameters

?? I'm lost? What would this be used for?

 

-- drew a large random sample from that distribution given those parameters
-- evaluated the log-likelihood given the true parameters and verified that the average score vector was sufficiently close to a zero vector

Did you want to draw samples, then feed that back into the tests? Take a look at the test framework now and see if you can get it going. If that's what you really want to do, if you can get it working or even specified, that'd be great. There are some engineering issues that you're going to need to tackle to make that work. It's a little too much for me to tackle right now, but if you can get it done, I'm all for it.

That would go a long way toward verifying that everything for that distribution is coherent.

The easy alternative is to manually generate those tests by hand. Given that there are only a handful of distributions (< 100), it's feasible to do manually, although it won't scale.

Ben

 

On Sunday, June 2, 2013 12:31:53 AM UTC-4, Peter Li wrote:

That's great news. I'm not sure how exactly the distribution testing
framework works, but I'd be glad to help / will look into it.


Peter

On 6/2/13 12:18 AM, Daniel Lee wrote:
> I think I've got it working using our old code. More to come. We need
> to do some more testing. Peter, want to help me out? The test
> distribution tests should test fvars too.
>
> Also, I think it's not efficient and it probably could be made more
> efficient.
>
>
> Daniel

--

[Peter Li]

I actually just went through and tested all of the fvar functions
(except the matrix ones--those are next) to make sure that fvar<fvar<T>>
and fvar<var> work and give the correct output; this is in the
clean_up_fvar branch.

>
> For var, we've customized the autodiff in the functions, so
> I think it's fair to test that. I think fvar<double> would
> suffice.
>
> The rule is that if there an autodiff type, all the autodiff
> types are the same. So I think we can get away with just
> testing fvar<double> in addition to var, which given the rule,
> should only double the total number of tests if you want to
> do it exhaustively.

I'm not sure this always works; I think all functions called in the
coded derivative must be implemented in our library (with a derivative
too). When I was testing the functions, anything involving the gamma
function did not work because we didn't have the digamma function
implemented (so I added it).

[Daniel Lee]

> On 6/4/13 12:22 PM, Bob Carpenter wrote:
>>
>>
>> On 6/4/13 12:14 PM, Daniel Lee wrote:
>>> I've pushed my changes to feature/fvar_operands_and_partials. I've added a single test for normal at
>>> src/test/agrad/fwd/prob/univariate/continuous/normal_test.cpp.
>>>
>>> I think the right thing to do is to extend the testing framework to generate every valid combination of tests. For
>>> normal_log, that means each of the 3 arguments can be:
>>> 1. double
>>> 2. var
>>> 3. fvar<double>
>>> 4. vector<double>
>>> 5. vector<var>
>>> 6. vector<fvar<double> >
>>> 7. Eigen::Matrix<double, 1, -1> (row_vector)
>>> 8. Eigen::Matrix<var, 1, -1> (row_vector)
>>> 9. Eigen::Matrix<fvar<double>, 1, -1> (row_vector)
>>> 10. Eigen::Matrix<double, -1, 1> (vector)
>>> 11. Eigen::Matrix<var, -1, 1> (vector)
>>> 12. Eigen::Matrix<fvar<double>, -1, 1> (vector)
>>>
>>> Except, we won't have all 12^3 combinations. If I understood the induction arguments correctly, we can't have
>>> fvar<double> and var in the same argument list.
>>>
>>> Should we add fvar<var> in there too? That's going to get complicated.
>>
>> Yes, and fvar<fvar<var> > and Matrix<fvar<fvar<var> > > and so on.
>>
>> I think we'll just have to trust that the autodiff
>> is working because we unit test that independently.
>>

>> For var, we've customized the autodiff in the functions, so
>> I think it's fair to test that. I think fvar<double> would
>> suffice.

I'm not sure it'll work yet. I'm not confident that we have the details worked out for the templated type inside the fvar. I want to verify that the gradients aren't disconnected through our distributions.

[Ben Goodrich]

On Tuesday, June 4, 2013 12:25:12 PM UTC-4, Daniel Lee wrote:

?? I'm lost? What would this be used for?

 
Here is what it would look like in R code for the standard normal distribution:

# Pick population parameters
# In Stan, these would be var types rather than doubles
# although the the foo_rng() functions only take doubles I think
mu <- 0
sigma <- 1

# Generate large sample from this population
SIMS <- 100000
X <- rnorm(SIMS, mu, sigma)

# Calculate derivatives
# In Stan, this would entail evaluating the gradient of normal_log(X, mu, sigma)
# In R, I just do the gradient analytically
score <- function(x, mu, sigma) c( (x - mu) / sigma^2, (mu^2 - sigma^2 + x^2 - 2 * mu * x) / sigma^3 )
partials <- t(sapply(X, score, mu = mu, sigma = sigma))
colMeans(partials) # small because at the mode

But even better would be to extract the Hessian when normal_log(X, mu, sigma) is evaluated in Stan, get the inverse of the estimated information matrix, and conduct a score test. In R, this is like

I_inv <- diag(c(sigma^2 / SIMS, 2 * sigma^4 / SIMS)) # known analytically
grad <- colSums(partials)
grad %*% I_inv %*% grad # ~ chi-squared with length(grad) degrees of freedom

So, if that test is rejected, then either

-- SIMS was too small
-- there is a rare Type I error
-- the foo_rng function in Stan is wrong
-- the gradient is wrong
-- the Hessian is wrong

But it would be better if this could be automated as much as possible. Definitely drawing the sample, extracting the gradient and Hessian, and conducting the score test is automatable. Setting the population parameters is a bit tricky because different distributions have different numbers of parameters and different restrictions on them.

Ben

[Bob Carpenter]

On 6/4/13 12:26 PM, Peter Li wrote:
>
> On 6/4/13 12:22 PM, Bob Carpenter wrote:

>> The rule is that if there an autodiff type, all the autodiff
>> types are the same.

This part of my reasoning is OK.

>> So I think we can get away with just
>> testing fvar<double> in addition to var, which given the rule,
>> should only double the total number of tests if you want to
>> do it exhaustively.

> I'm not sure this always works; I think all functions called in the coded derivative must be implemented in our library
> (with a derivative too). When I was testing the functions, anything involving the gamma function did not work because we
> didn't have the digamma function implemented (so I added it).

Good point --- you're absolutely right. And the digamma
function derivative involves a trigramma. It's turtles all the
way down.

- Bob

[Daniel Lee]

On Tue, Jun 4, 2013 at 10:15 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Good point --- you're absolutely right.  And the digamma
function derivative involves a trigramma.  It's turtles all the
way down.

So... how many layers of turtles should we test?
fvar<fvar<double> >

fvar<fvar<var> >

So for something like fvar<var>, the types to test are (without vectorization):

fvar<var>, fvar<var>, fvar<var>

fvar<var>, fvar<var>, var

fvar<var>, var, fvar<var>

fvar<var>, var, var

var, fvar<var>, fvar<var>

var, fvar<var>, var

var, var, fvar<var>

But no doubles in there. If it's fvar<fvar<double> >, no doubles, just fvar<double>?

[Bob Carpenter]

On 6/4/13 10:21 PM, Daniel Lee wrote:

We don't need to mix fvar<var> and var. If there's
an autodiff type, all the autodiff types must be the
same (sort of like vector lengths).

But you do need double. So just replace the bare
var with double and that's what we want to test for
fvar<var>.

I'd say the place to start is with the derivatives
we're going to need for RM-HMC. That's all the
things that are in autodiff.hpp.

I should add some doc for that indicating which
instantiations we need for each function concept.

To get all of these operations with double outputs,
we need these:

autodiff.hpp function autodiff scalar order
----------------------- --------------- -----
gr_tr_mat_times_hessian fvar<fvar<var>> 3
hessian fvar<var> 2
hessian_times_vector fvar<var> 2
gradient_dot_vector var 1
gradient var 1

jacobian var 1
fvar<double> 1

We don't have any use for the Jacobian planned in
the immediate future.

So that means we need tests for

fvar<var>
fvar<fvar<var>>

There is also:

fvar<double>
fvar<fvar<double>>

But we don't really foresee any use for these because
fvar<double> only computes a single derivative (dx1) and
fvar<fvar<double>> only a pair of derivatives and a second
derivative (dx1, dx2 and dx1.dx2).

And it's not just distributions, it's all of our functions
that we need to test.

We don't have second or third order finite diffs for testing.

- Bob

[Peter Li]

Oops. When I added digamma, I didn't include the correct derivative. I'm
not sure how we could add the trigamma function--it's not in boost and
doesn't have a closed form.

Anyone have any ideas?

[Bob Carpenter]

On 6/5/13 6:54 PM, Peter Li wrote:
> Oops. When I added digamma, I didn't include the correct derivative. I'm not sure how we could add the trigamma
> function--it's not in boost and doesn't have a closed form.
>
> Anyone have any ideas?

Use autodiff on digamma? Is it templated well enough
we can do that? You can always test with finite diffs.

- Bob

[Peter Li]

I'm guessing you mean like how we were able to use autodiff for a lot of
the matrix ops, but the reason I thought to implement digamma in the
first place was because the autodiff wasn't working directly when I
tried calling fvar<fvar<T>> or fvar<var> of any of the functions
involving the gamma function. I've been looking for a computational
method for trigamma, but I've had no luck yet.


Peter

[Bob Carpenter]

On 6/6/13 11:36 AM, Peter Li wrote:
> I'm guessing you mean like how we were able to use autodiff for a lot of the matrix ops, but the reason I thought to
> implement digamma in the first place was because the autodiff wasn't working directly when I tried calling fvar<fvar<T>>
> or fvar<var> of any of the functions involving the gamma function. I've been looking for a computational method for
> trigamma, but I've had no luck yet.

There has to be an iterative method for it. But then
you have the issue of differentiating it.

I'd suggest taking a look at what Boost does for
gamma or digamma and seeing if we can't make an
auto-differentiable version. We can make sure it
only gets called when absolutely necessary (e.g., to
compute the trigamma).

- Bob

[Bob Carpenter]

There's also an implementation in R you could check out.

We can't use it directly because of the GPL licensing.

- Bob

On 6/6/13 11:36 AM, Peter Li wrote:

[Michael Betancourt]

You calculate it the same way you calculate digamma -- you use the recurrence relation to relate phi(x) to phi(x') where x' >> 1 (10 would be more than enough, I believe 5 or 6 is recommended in the literature) and then use the Laurent expansion for phi(x), x >> 1.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Peter Li]

Alright. The version that I've added is based on what Michael said and http://people.sc.fsu.edu/~jburkardt/cpp_src/asa121/asa121.cpp, but I'm unsure of how the copyright would work out.


Peter

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Marcus Brubaker]

That code is LGPLd so cannot be included in Stan.

Cheers,

Marcus

[Bob Carpenter]

Copyright is for a specific rendering of the algorithm in
code, not for the algorithm itself.

The author of the code you link cites this paper as the
source of the algorithm:

// BE Schneider,
// Algorithm AS 121:
// Trigamma Function,
// Applied Statistics,
// Volume 27, Number 1, pages 97-99, 1978.

So I'd suggest just writing a fresh C++ version based on
the journal algorithm above.

There's also a remark, which I haven't read, that might be
worth reading:

Remark AS R88: A Remark on Algorithm AS 121: The Trigamma Function
Brian J. Francis
Journal of the Royal Statistical Society. Series C (Applied Statistics)
Vol. 40, No. 3 (1991), pp. 514-515
(article consists of 2 pages)
Published by: Wiley for the Royal Statistical Society

And make sure it's templated and differentiable itself so we can do
higher-order derivatives.

At that point, it's probably not going to look like the other
code very much, so we should be OK on copyright.

- Bob

>> stan-dev+u...@googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.

>> For more options, visit https://groups.google.com/groups/opt_out.
>>
>>
>>
>> --
>> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to

>> stan-dev+u...@googlegroups.com <mailto:stan-dev+u...@googlegroups.com>.

>> For more options, visit https://groups.google.com/groups/opt_out.
>>
>>
>
> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to

> stan-dev+u...@googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.

[Marcus Brubaker]

Of course.  I was just pointing out that if he had copied code from that file, it wouldn't be kosher.  A "clean-room" version based on the paper is perfectly valid.

Cheers,

Marcus

        stan-dev+unsubscribe@googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.

        For more options, visit https://groups.google.com/groups/opt_out.



    --
    You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
    To unsubscribe from this group and stop receiving emails from it, send an email to

    stan-dev+unsubscribe@googlegroups.com <mailto:stan-dev+unsubscribe@googlegroups.com>.

    For more options, visit https://groups.google.com/groups/opt_out.

    --
    You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
    To unsubscribe from this group and stop receiving emails from it, send an email to

    stan-dev+unsubscribe@googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.

    For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Peter Li]

The paper (http://lib.stat.cmu.edu/apstat/121) looks very similar to the code in the link that I sent earlier.


Peter

        stan-dev+u...@googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.

        For more options, visit https://groups.google.com/groups/opt_out.



    --
    You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
    To unsubscribe from this group and stop receiving emails from it, send an email to

    stan-dev+u...@googlegroups.com <mailto:stan-dev+u...@googlegroups.com>.

    For more options, visit https://groups.google.com/groups/opt_out.

    --
    You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
    To unsubscribe from this group and stop receiving emails from it, send an email to

    stan-dev+u...@googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.

    For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Bob Carpenter]

Then I think we're OK if you recoded the algorithm in
C++. Is there a branch with the code in it I can look
at to compare?

- Bob

>> <mailto:li.pe...@gmail.com <mailto:li.pe...@gmail.com>>>

>> wrote:
>>
>> Alright. The version that I've added is based on what
>> Michael said and
>> http://people.sc.fsu.edu/~jburkardt/cpp_src/asa121/asa121.cpp

>> <http://people.sc.fsu.edu/%7Ejburkardt/cpp_src/asa121/asa121.cpp>,

>> but I'm unsure of how the copyright would work out.
>>
>>
>> Peter
>>
>>
>> On 6/6/13 12:44 PM, Michael Betancourt wrote:
>>
>> You calculate it the same way you calculate digamma --
>> you use the recurrence relation to relate phi(x) to phi(x')
>> where x' >> 1 (10 would be more than enough, I believe
>> 5 or 6 is recommended in the literature) and then use the
>> Laurent expansion for phi(x), x >> 1.
>>
>>
>> On Thu, Jun 6, 2013 at 5:41 PM, Bob Carpenter
>> <ca...@alias-i.com <mailto:ca...@alias-i.com>

>> <mailto:stan-dev%2Bunsu...@googlegroups.com
>> <mailto:stan-dev%252Buns...@googlegroups.com>>.

>> For more options, visit
>> https://groups.google.com/groups/opt_out.
>>
>>
>>
>> --
>> You received this message because you are subscribed
>> to the Google Groups "stan development mailing list" group.
>> To unsubscribe from this group and stop receiving
>> emails from it, send an email to
>> stan-dev+u...@googlegroups.com
>> <mailto:stan-dev%2Bunsu...@googlegroups.com>

>> <mailto:stan-dev+u...@googlegroups.com

>> <mailto:stan-dev%2Bunsu...@googlegroups.com>>.
>> For more options, visit
>> https://groups.google.com/groups/opt_out.
>>
>>
>>
>> --
>> You received this message because you are subscribed to
>> the Google Groups "stan development mailing list" group.
>> To unsubscribe from this group and stop receiving emails
>> from it, send an email to
>> stan-dev+u...@googlegroups.com

>> <mailto:stan-dev%2Bunsu...@googlegroups.com>
>> <mailto:stan-dev%2Bunsu...@googlegroups.com
>> <mailto:stan-dev%252Buns...@googlegroups.com>>.

>> For more options, visit
>> https://groups.google.com/groups/opt_out.
>>
>>
>>
>> --
>> You received this message because you are subscribed to the
>> Google Groups "stan development mailing list" group.
>> To unsubscribe from this group and stop receiving emails from
>> it, send an email to stan-dev+u...@googlegroups.com
>> <mailto:stan-dev%2Bunsu...@googlegroups.com>.
>> For more options, visit https://groups.google.com/groups/opt_out.
>>
>>
>>
>> --
>> You received this message because you are subscribed to the Google
>> Groups "stan development mailing list" group.
>> To unsubscribe from this group and stop receiving emails from it,

>> send an email to stan-dev+u...@googlegroups.com

[Peter Li]

It's in feature/clean_up_fvar in stan/math/functions/trigamma.hpp.


Peter