How to use ODE solver to convolve two random variables

[Christopher Fisher]

Hi all-

The likelihood function for my model is a convolution of a normal and lognormal distribution. I am having difficulty setting up the function to pass to integrate_ode_rk45() and figuring out which arguments are necessary for my purposes. I think the setup is unclear to me because its generally used for differential equations. 

I need to pass a vector of parameters (theta), data point (rt) and something to indicate the limits of y (0 to + Inf) to Convolve() via ntegrate_ode_rk45() . Here is what I have so far:

real[] Convolve(real[] y,real[] theta,real rt) {
          real dydt[1];
          real mu; //mean of normal
          real sigma; //sd of normal
          real act; //negative mean of lognormal
          real noise; //sd of lognormal
          mu = theta[1];
          sigma = theta[2];
          act = -theta[3];
          noise = (.15*3.1415926535897)/sqrt(3);
          //Convolve normal and lognormal distributions
          dydt[1] = normal_lpdf(rt - y[1] | mu, sigma) + lognormal_lpdf(y[1] |act,noise);
          return dydt;
        }

Any help would be greatly appreciated.

-Chris 

[Sebastian Weber]

Hi!

I always thought this is a cool application for the ODE integrator. Attached is a working example of what I think you intend to do.

Sebastian

[Christopher Fisher]

Thanks for the worked example. Much appreciated!

I haven't tested it yet, but it looks about right. And for future reference, I will post updated code if I need to modify anything. 

I'm currently working with computational models that boil down to mixtures of convolutions. So the ODE solver should be useful and make for an interesting application of Stan. I hope to leverage the forthcoming discrete Fourier transform functionality for scaling up to more complex models. 

Thanks again,

Chris 

[Christopher Fisher]

Hi Sebastian-

Out of curiosity, do you know whether its possible to integrate across multiple variables? I tried extending my model to the convolution of 3 variables and received the following error: 

first argument to integrate_ode_rk45 must be the name of a function with signature (real, real[], real[], real[], int[]) : real[] 

I tried several things. Here is one thing I tried:

real[] kernel(real[] x, real[] y, real[] theta, real[] x_r, int[] x_i) {
real dydx[1];
real mu = theta[1]; //mean of normal
real sigma = theta[2]; //sd of normal
real act1 = -theta[3]; //negative mean of lognormal
real act2 = -theta[4]; //negative mean of lognormal
real noise = (.15*3.1415926535897)/sqrt(3); //sd of lognormal

//Convolve normal and lognormal distributions

dydx[1] = exp(normal_lpdf(x_r[1] - x[1] -x[2]| mu, sigma) + lognormal_lpdf(x[1] |act1,noise) + lognormal_lpdf(x[2] |act2,noise));
return dydx;
} real convolve_kernel(real a, real b, real[] theta, real x) {
int x_i[0];
return(integrate_ode_rk45(kernel, rep_array(0., 2), a, rep_array(b, 1), theta, rep_array(x, 1), x_i)[1,1]);
}

I'm guessing its not possible or I am not using the right inputs. 

Thanks,

Chris 

[Bob Carpenter]

See below.

> On Feb 7, 2017, at 5:00 PM, Christopher Fisher <fish...@miamioh.edu> wrote:
>
> Hi Sebastian-
> Out of curiosity, do you know whether its possible to integrate across multiple variables? I tried extending my model to the convolution of 3 variables and received the following error:
>
> first argument to integrate_ode_rk45 must be the name of a function with signature (real, real[], real[], real[], int[]) : real[]

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Here's the example from page 262 of the manual:

real[] sho(real t,

real[] y,
real[] theta,
real[] x_r,
int[] x_i) {

real dydt[2];
dydt[1] = y[2];
dydt[2] = -y[1] - theta[1] * y[2];
return dydt;
}

sho is a function with the appropriate signature. I don't know what your x is
supposed to be, and here's the call in the code:


y_hat = integrate_ode_rk45(sho, y0, t0, ts, theta, x_r, x_i);


- Bob

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

[Sebastian Weber]

Hi Chris!

You can integrate over multiple variables at the same time, yes. You just need to define a larger dydx vector. However, I can only recommend doing this if and only if those three variables are somehow dependent on one another.

I think what you want to do is just call the routine three times with different parameters. The reason is that you always want to keep the sensitivity system as small as possible; have a read on ODEs in the manual and if you need more details, please read the arxiv paper on autodiff in Stan; it has a specific chapter on ODEs. Once you read this and made yourself familiar with the crazy numerical cost of ODEs, you will automatically use ODEs in Stan very carefully.

Sebastian

[Christopher Fisher]

Thanks for your replies Sebastian and Bob. It seems like the ODE lacks the functionality I need and, as Sebastian noted, it would not be very scalable. What I have done so far with Cubature seems to provide further evidence of very poor scalability. 

Ultimately, discrete Fourier transform (and its inverse) seems to be the ideal solution based on what I have gathered.  I'm trying to understand if Stan will have the tools for convolution via DFT.  A few ideas were proposed in the development forum for handling complex numbers, but the plan to implement one, if any, seemed ambiguous. So I was wondering if either of you know more about the plans? If Stan won't have this functionality in the near future, is there a way to call C++ code or a library to perform the DFT? Ultimately all I need is to return a probability density. My model would not need complex numbers directly. 

Thanks for your help,

-Chris 

[Bob Carpenter]

Yes, we're working on it right now and some form will be out and
we'll add the other forms as requested. So it sounds like we want
something to return the pair of complex and real components.
It may take a couple releases before all the pieces are there.

- Bob

[Sebastian Weber]

rstan and cmdstan now allow you to call your own C++ functions; see the vignettes and the respective manual. Of course, this is for the very experienced ones... and going down this route comes at the price of potentially having to change you code with a new release - no guarantees are given in terms of release-to-release compatibility here. However, should you embark on this, then I guess that you are more than welcome to comment on the current efforts to get a DFT into Stan (comments, doc, test cases, code). Specific examples are always a welcome benchmark.

Sebastian

[Bob Carpenter]

I'm working on includes for Stan programs, but we also need to
work on making includes of C++ part of the core language.

The tricky part is going to be a reliable and safe way to what's
being hacked now --- declaring a function header in the functions
block and then hoping it gets compiled and linked with the C++.
The declaration adds it to the symbol table so you can use it in
the Stan program; the C++ adds in the implementation. I don't have
a good idea on how to improve what we have now, so if anyone
has any suggestions, I'd love to hear them.


- Bob