Numeric integration

[Marco Inacio]

Is there any function in Stan that evaluates a definite univariate integral numerically (like integrate in R)? I need to calculate a normalizing factor and reusing existing code would be great.

[Bob Carpenter]

No, but if you want to build such a thing, it'd be
a great addition to Stan! It'd probably have to work
the same way as the integrate_ode() function does. I'd
be happy to help with the language parts of the coding.

There are a couple other things like this that Ben
wants to implement, so maybe he can say more.

- Bob

> On Nov 30, 2015, at 12:24 PM, Marco Inacio <marcoig...@gmail.com> wrote:
>
> Is there any function in Stan that evaluates a definite univariate integral numerically (like integrate in R)? I need to calculate a normalizing factor and reusing existing code would be great.
>

> --
> 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/d/optout.

[Marco Inacio]

I wasn't planning on creating a generic numerical integration function,
but just the one I need a function I need [which comes from orthogonal
series density estimation] which has simple derivative.

But I find that idea interesting, however how would we calculate the
gradient for the integral of a generic functions?

[Ben Goodrich]

On Monday, November 30, 2015 at 2:53:37 PM UTC-5, Bob Carpenter wrote:

There are a couple other things like this that Ben
wants to implement, so maybe he can say more.

I prototyped a little code based on

http://www.codeproject.com/Articles/31550/Fast-Numerical-Integration

which seemed as if it would work for a particular PDF that happens to involve an integral of a function that lends itself to this quadrature algorithm.

Ben

[Bob Carpenter]

I'm not sure, but I think Ben and Michael know how to do
this. We also have the issue of wanting a definite integral
but allowing both a parametric function and bounds of
+/- infinity on the integration.

- Bob

[Michael Betancourt]

https://en.wikipedia.org/wiki/Differentiation_under_the_integral_sign

[Marco Inacio]

My issue was more about implementation than calculus, but I think I can
see how to do it now. Thanks.

[Marco Inacio]

This seems good. I was trying this with the old but good quadpack, but it was hard since I had never worked with importing external Fortran code.

Is your code in a Stan branch?

[Ben Goodrich]

On Monday, November 30, 2015 at 4:39:45 PM UTC-5, Marco Inacio wrote:

This seems good. I was trying this with the old but good quadpack, but it was hard since I had never worked with importing external Fortran code.

Is your code in a Stan branch?

No, I just compiled his C++ code with the PDF I was thinking about implementing.

Ben

[Bob Carpenter]

I don't think we're going to be supporting Fortran linkage
any time soon --- we have enough build headaches as it is.

The simplest thing to do would probably to get a boost
integrator, then just autodiff it to start and make sure
you get reasonable results, then do the analytic derivatives.

Then we can go looking for better integrators, etc. when
the analytic bits are sorted out.

- Bob

[Marco Inacio]

Why not the one Ben Goodrich suggested? It's C++ and seems pretty
simple. Also, I couldn't find boost numerical integrator.

[Bob Carpenter]

All I could find in Boost was the ODE solver, too. Maybe they don't
have a simple definite integrator. The reason I suggested it is licensing
and minimizing dependencies in Stan. If you're just doing this for
a standalone side project, neither of these issues matters.

- Bob

[Marco Inacio]

Licensing is not a problem at all since it's public domain.

What looks troublesome to me now is how to pass the additional arguments to f:

integrate_function(function f, real lower_limit, real upper_limit, ...)

What would "..." be? A single real? A single eigen vector? Variadic template (!!!)?

[Bob Carpenter]

"public domain" is not a license. Whoever wrote the
code owns the copyright and needs to put a proper license on
it.

:-) OK, I followed the link for the "public domain" business
and it links to Creative Commons, who now disavow the whole
public domain business (for good reason!).

It's John Cook, and he seems like a nice guy from his blog,
so maybe he'd be willing to license it under BSD.

- Bob

[Andrew Gelman]

Should I cc Cook on this discussion?

P.S. Can’t we just get public domain by waiting about 150 years and hoping that none of us have descendants who marry anyone in the Disney family?

[Ben Goodrich]

On Monday, November 30, 2015 at 6:15:08 PM UTC-5, Marco Inacio wrote:

Why not the one Ben Goodrich suggested? It's C++ and seems pretty
simple. Also, I couldn't find boost numerical integrator.

Boost doesn't have a general numerical integrator. I think it would be difficult to expose a function like integrate() into the Stan language because no single numerical integration procedure works well on every problem. But Cook's approach works well on functions that decline at a double exponential rate, which covers a lot of things one might want to do in statistics. So, I think something like this could be included in Stan Math and used when appropriate as the back end by functions exposed to the Stan language.

Ben

[Bob Carpenter]

If we can specify what kinds of problems it applies to, I
see no problem introducing something that's limited. We can
add different variants as we go.

I'd think doubly exponential [i.e., exp(-exp(abs(x)))]
would be very limiting, but I have no clue as to what kinds
of applications people would use this for.

Cook points to more of an explanation here:

http://www.johndcook.com/blog/double_exponential_integration/

- Bob

P.S. + :-) on public domain and Disney heirs.

[Ben Goodrich]

On Monday, November 30, 2015 at 10:47:23 PM UTC-5, Bob Carpenter wrote:

I'd think doubly exponential [i.e., exp(-exp(abs(x)))]
would be very limiting, but I have no clue as to what kinds
of applications people would use this for.

Double-exponential after the transformation. Equation 4 in

http://academic2.american.edu/~jpnolan/stable/density.ps

is what I was looking to implement with it.
 

Cook points to more of an explanation here:

  http://www.johndcook.com/blog/double_exponential_integration/

Better explanation with links to papers by Bailey

https://en.wikipedia.org/wiki/Tanh-sinh_quadrature

Ben

 

 

> On Nov 30, 2015, at 10:34 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> On Monday, November 30, 2015 at 6:15:08 PM UTC-5, Marco Inacio wrote:
> Why not the one Ben Goodrich suggested? It's C++ and seems pretty
> simple. Also, I couldn't find boost numerical integrator.
>
> Boost doesn't have a general numerical integrator. I think it would be difficult to expose a function like integrate() into the Stan language because no single numerical integration procedure works well on every problem. But Cook's approach works well on functions that decline at a double exponential rate, which covers a lot of things one might want to do in statistics. So, I think something like this could be included in Stan Math and used when appropriate as the back end by functions exposed to the Stan language.
>
> Ben
>
> --
> 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.

[Bob Carpenter]

Cool --- thanks. I was wondering about transforms.
I'll have to look through the details when I get
a chance.

- Bob

> > 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/d/optout.
>
>

> --
> 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]

OK, I took a quick look. It looks like it works for integrals
on (-1, 1), so I see how it'd work for bounded integrals by scaling
and translating coordinates linearly.

But how would it work with integrals on (-inf, inf) or (0, inf)?
Would those have to be double exponential to start with? Or
could we take a normal(0, 1), transform it to (-1, 1), then
work from there with a Jacobian?

Also, it looks like Bailey himself wrote a BSD version in C++:

http://crd-legacy.lbl.gov/~dhbailey/mpdist/

though it's not clear which of those packages it'd be in.

- Bob

> > 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/d/optout.
>
>

> --
> 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.

[Ben Goodrich]

On Monday, November 30, 2015 at 11:12:11 PM UTC-5, Bob Carpenter wrote:

But how would it work with integrals on (-inf, inf) or (0, inf)?
Would those have to be double exponential to start with?  

If it were double exponential to start with, then it would be okay but I think a different approach would often be needed with already unbounded integrals.
 

Or could we take a normal(0, 1), transform it to (-1, 1), then
work from there with a Jacobian?

Maybe if the tails were well-behaved.
 

Also, it looks like Bailey himself wrote a BSD version in C++:

  http://crd-legacy.lbl.gov/~dhbailey/mpdist/

though it's not clear which of those packages it'd be in.

Cook's version is fine, although Bailey may have some other good stuff in ARPREC. The problem with public domain is that nothing prevents someone else (like us) from slapping their own license on it. We'd have to rewrite a good bit of it anyway to make it work with Stan.

Ben

> On Nov 30, 2015, at 11:00 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> On Monday, November 30, 2015 at 10:47:23 PM UTC-5, Bob Carpenter wrote:
> I'd think doubly exponential [i.e., exp(-exp(abs(x)))]
> would be very limiting, but I have no clue as to what kinds
> of applications people would use this for.
>
> Double-exponential after the transformation. Equation 4 in
>
> http://academic2.american.edu/~jpnolan/stable/density.ps
>
> is what I was looking to implement with it.
>  
> Cook points to more of an explanation here:
>
>   http://www.johndcook.com/blog/double_exponential_integration/
>
> Better explanation with links to papers by Bailey
>
> https://en.wikipedia.org/wiki/Tanh-sinh_quadrature
>
> Ben
>  
>  
> > On Nov 30, 2015, at 10:34 PM, Ben Goodrich <goodri...@gmail.com> wrote:
> >
> > On Monday, November 30, 2015 at 6:15:08 PM UTC-5, Marco Inacio wrote:
> > Why not the one Ben Goodrich suggested? It's C++ and seems pretty
> > simple. Also, I couldn't find boost numerical integrator.
> >
> > Boost doesn't have a general numerical integrator. I think it would be difficult to expose a function like integrate() into the Stan language because no single numerical integration procedure works well on every problem. But Cook's approach works well on functions that decline at a double exponential rate, which covers a lot of things one might want to do in statistics. So, I think something like this could be included in Stan Math and used when appropriate as the back end by functions exposed to the Stan language.
> >
> > Ben
> >
> > --
> > 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/d/optout.
>
>
> --
> 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.

[Bob Carpenter]

Yeah, I figured it'd work with already double exponential but
not otherwise. But it looks awesome enough to normalize a beta
distribution.

For licensing, public domain is problematic. See:

https://en.wikipedia.org/wiki/Public_domain#Dedicating_works_to_the_public_domain

and

https://creativecommons.org/about/cc0

The latter, I think, is meant to replace Cook's original intent. The
thing he pointed to is now pointing to a "retired legal tool":

http://creativecommons.org/licenses/publicdomain/

- Bob

> > > 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/d/optout.
> >
> >
> > --
> > 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/d/optout.
>
>
> --
> 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.

[Marco Inacio]

I was able to create the numerical integration function when the limits of integral are double, and the additional parameter is a double or a var.

To further generalize it to vector of double or var I would like to know if there is a function in Stan to convert std::vector<var> to std::vector<double> just like value_of does (if there's isn't, no problem, I'll write one, just asking to avoid reinventing the wheel).

Also, I wrote the code using OperandsandPartials to make it easier to vectorize, and my function do get the value of gradient correct, but somehow after I pass it to OperandsandPartials it gets messed up. Any ideas on why?


    template <typename F, typename T_param>
    inline
    T_param integrate_function(const F& f,
                               const double a,
                               const double b,
                               const T_param& param) {
...
        for (size_t n = 0; n < N; n++) {

            double tempp =
            get_integrate_val(boost::bind<double>(get_integrate_gradient<F, T_param>, f, _1, param, param_vec[n]), a, b, 1e-6);
            //THIS DOES PRINTS THE CORRECT VALUE!
            std::cout << "CALCULATED DERIVATIVE" << tempp << std::endl;

             //Then it gets messed up here
          operands_and_partials.d_x1[n] +=    tempp;    }

        return operands_and_partials.to_var(val_, param);




Output of finite_diff test:

CALCULATED DERIVATIVE0.5
test/unit/math/rev/arr/functor/integrate_function_test.cpp:67: Failure
Value of: g[0]
  Actual: 1.5
Expected: (integrate_function(if1, .2, .7, .6+1e-6) - integrate_function(if1, .2, .7, .6-1e-6))/2e-6
Which is: 0.5

CALCULATED DERIVATIVE1.63731
test/unit/math/rev/arr/functor/integrate_function_test.cpp:81: Failure
Value of: g[0]
  Actual: 4.6111894
Expected: (integrate_function(if2, 0., 1.1, .68+1e-6) - integrate_function(if2, 0., 1.1, .68-1e-6))/2e-6
Which is: 1.6373117


(full code attached)

[Bob Carpenter]

> On Dec 1, 2015, at 7:12 PM, Marco Inacio <marcoig...@gmail.com> wrote:
>
> I was able to create the numerical integration function when the limits of integral are double, and the additional parameter is a double or a var.
>
> To further generalize it to vector of double or var I would like to know if there is a function in Stan to convert std::vector<var> to std::vector<double> just like value_of does (if there's isn't, no problem, I'll write one, just asking to avoid reinventing the wheel).

Not that I know of. The right thing to do would be to overload
value_of, I'd think.

I'm not sure what you mean by passing it to OperandsAndPartials.
You create an OperandsAndPartials then populate it. If the
values are right that you are trying to put into OperandsAndPartials,
it's then probably an issue with sizing. It's a complicated
object that increments by default.

- Bob

> --
> 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/d/optout.

> <integrate_function.hpp><integrate_function_test.cpp>

[Marco Inacio]

> The problem with public domain is that nothing prevents someone else
> (like us) from slapping their own license on it.
>

But would this forbid anyone else to use the original public domain code?

[Bob Carpenter]

Read the links I used in the previous reply --- it'll fill
you in on why "public domain" is a bad idea when someone
actually could hold the copyright (that is, it's not someone
who's been dead a long time).

- Bob

> On Dec 1, 2015, at 10:52 PM, Marco Inacio <marcoig...@gmail.com> wrote:
>
>
>> The problem with public domain is that nothing prevents someone else (like us) from slapping their own license on it.
>>
> But would this forbid anyone else to use the original public domain code?
>

> --
> 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.

[Marco Inacio]

I understand that it creates a lot of legal complications, so how do we
proceed here?

[Bob Carpenter]

What do you want to do? Add it to Stan or just run
your own example? If it's the latter, it's no big deal.
If the former, and you really want to use this integrator,
then we should write to John Cook and see if he'll put
out a version under BSD.

We need to keep the licensing clean on Stan itself.
Sorry this keeps coming up with things you want to do!

- Bob

[Marco Inacio]

The problem is actually somewhere else: every time I run this function,
it returns the correct value plus the value previously returned. So if
the correct value is 1.40186, we get:

1.40186
2.80371
4.20557
5.60742
7.00928

and so on.

template <typename F, typename T_param>
inline

double get_integrate_gradient(const F& f,
const double x,
const T_param& param, //T_param is a
std::vector<var>
const var& param_n) { //One of the
elements of param

std::vector<double> gradient;
try {
stan::math::start_nested();
var dparams_vars = f(x, param);
std::vector<var> avec(1, param_n);
dparams_vars.grad(avec, gradient);
set_zero_all_adjoints_nested();
} catch (const std::exception& e) {
stan::math::recover_memory_nested();
throw;
}
stan::math::recover_memory_nested();
return gradient[0];
}

...
get_integrate_gradient<F, T_param>(f, .4, param, param_vec[0]);
get_integrate_gradient<F, T_param>(f, .4, param, param_vec[0]);
get_integrate_gradient<F, T_param>(f, .4, param, param_vec[0]);
get_integrate_gradient<F, T_param>(f, .4, param, param_vec[0]);
get_integrate_gradient<F, T_param>(f, .4, param, param_vec[0]);

[Marco Inacio]

I want to add it to Stan.

In fact, it's almost complete for fixed upper/lower integration limits
and rev-differentiation for (vector of) parameters.

Does integrate_ode has support for fwd-diff? I looked for it to learn
how to do nested auto-diff for fwd, but I couldn't find it.

[Krzysztof Sakrejda]

On Wednesday, December 2, 2015 at 1:27:30 AM UTC-5, Bob Carpenter wrote:
> What do you want to do? Add it to Stan or just run
> your own example? If it's the latter, it's no big deal.
> If the former, and you really want to use this integrator,
> then we should write to John Cook and see if he'll put
> out a version under BSD.

I've corresponded a little with John Cook before so I took the liberty of sending him a message with a pointer to this discussion. In the past he's been surprisingly responsive so we might just hear from him. Otherwise I'm happy to follow up or let you (Marco) do it. Let me know.

>
> We need to keep the licensing clean on Stan itself.
> Sorry this keeps coming up with things you want to do!

>
> - Bob

[Bob Carpenter]

That'd be great.

No, integrate_ode doesn't have support for forward-mode autodiff yet.

Andrew and I can write to John Cook and see if he'll release
the code via BSD.

- Bob

[Marco Inacio]

Weird. I'm out of clues. I don't understand why the adjoint of param
does not get cleared if param is std::vector<var>, but it does get
cleared if param is var.

try {
start_nested();
f(x, param).grad();
result = param[0].adj();

set_zero_all_adjoints_nested();
} catch (const std::exception& e) {

recover_memory_nested();
throw;
}
recover_memory_nested();
param[0].adj(); //not zero. Why?

//But param.adj() would be zero if param was var


(if I change set_zero_all_adjoints_nested() to set_zero_all_adjoints(),
then obviously it works.)

[Bob Carpenter]

It looks like you're defining param outside
of the start_nested(), so it's not going to get
recovered when you set the nested ones to zero.
Only the parameters defined after the start_nested()
will be set to zero.

- Bob

[Marco Inacio]

If I just copy somevar to another var inside nested this won't work,
right? What should I do instead?


var somevar;

start_nested()
var anothervar = somevar;
//do stuff
recover_memory_nested();

[Bob Carpenter]

I'm not sure what your intent is. The point of the
nested evaluation is to start a whole new autodiff
on top of the current autodiff stack, evaluate it,
then pop it off the stack. It doesn't make sense to
involve variables that are declared before start_nested()
is called --- think of it as forming a scope just like
in a program.

- Bob

[Marco Inacio]

I was actually interested only in the value of param. Anyway, I managed
to do it.

Could you tell me where should I put Cook's "DEIntegrator.h" and
"DEIntegratorConstants.h" files? Or if you preffer, organizing them
yourself, I created a clean branch
feature/issue-210-numerical-integration in stan-dev/math.

[Ben Goodrich]

On Thursday, December 3, 2015 at 10:14:50 AM UTC-5, Marco Inacio wrote:

Could you tell me where should I put Cook's "DEIntegrator.h" and
"DEIntegratorConstants.h" files? Or if you preffer, organizing them
yourself, I created a clean branch
feature/issue-210-numerical-integration in stan-dev/math.

They should go in the Math library under

stan/math/prim/scal/fun

for the general templated versions and change prim to rev and / or fwd if there are specializations.

The DEIntegrator.h should be renamed DEIntegrator.hpp in Stan so that it will be tested by make test-headers. The .h extension is used for things like DEIntegratorConstants.h that cannot be used in isolation.

Ben

[Bob Carpenter]

I think:

stan/math/prim/arr/functor

might make more sense. It'll take a functor argument and probably
an array (std::vector) for parameters. And
that's where the rest of the ODE stuff's at.

Anyway, don't worry about it too much --- it's easy enough
to move things around later. The bigger deal is putting
a double-only version there, then a version in fwd and rev in
the same path.

- Bob