Hierarchical model involving complex numbers

[Dennis Conway]

Hello all,

I'm currently working on a model with a custom probability density calculation which involves complex numbers - which, as far as I can tell, is not supported by Stan. I was wondering however if I were to write my own density distribution in C++, would the automatic differentiator be able to handle the complex values, or would it break down?

Thanks for your help,

Dennis

EDIT: Sorry I don't think I made it clear, none of my model parameters need be complex - the complex numbers are just used in calculating the logp__ from my parameters.

[Daniel Lee]

Hi Dennis,

You should be able to write your own probability density calculation with complex numbers by representing the real and imaginary parts with two parameters, unless I'm missing something fundamental.

Re: C++, do you have a class that represents a complex number? I don't think there's a straightforward way to go from the Stan language, extend it to a new type, have that type represented in C++ somehow, pass that through functions with different signatures, then get that back out to the interfaces. But... if you wanted to write C++, that's easy. If you also represented complex numbers as 'vector[2]' in the Stan language, that should work as a proof of concept.

If you have a description of this density function, that would help. Either code or latex will do.

Daniel

On Fri, Apr 17, 2015 at 10:19 AM, Dennis Conway <denni...@gmail.com> wrote:

Hello all,

I'm currently working on a model with a custom probability density calculation which involves complex numbers - which, as far as I can tell, is not supported by Stan. I was wondering however if I were to write my own density distribution in C++, would the automatic differentiator be able to handle the complex values, or would it break down?

Thanks for your help,

Dennis

--
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.
For more options, visit https://groups.google.com/d/optout.

[Dennis Conway]

Hey Daniel, thanks for your help on this.

I think it would be quite tricky to split up the imaginary and real parts into two parameters, as there are a few multiplications in my formula which I would need to write in terms of stretchings and rotations. The calculation involves a recursive formula, which I've put at the bottom of this message, and the probability is given by a misfit between the output of the recursive formula and some data that I have. If this is straightforward enough to write up in C++ then I will have a go at it, as I'm curious to see how my model (high dimensionality with highly correlated parameters) performs with the NUTS sampler.

Thanks,

Dennis

t is real

c&q are complex valued

c_{n} = \frac{c_{n-1}\cdot q_{n}+\tanh(tq_{n})}{1+c_{n-1}\cdot q_{n}\cdot \tanh(tq_{n})}\frac{1}{q_{n}}

EDIT: Oh right, I just realised that multiplication and division are not as hard as I was imagining... I can just break down the real and imaginary parts for each operation and follow where they go - I guess this is what you were saying but my maths is a bit rusty!

[Bob Carpenter]

To do that in Stan you'd need to define functions for all the
complex operations and treat the vector as a 2 x N vector with
a real and imaginary component. Then you could write that recursive
equation (you'd need a base case, though!).

I don't know that it'll be easier to do in C++. The point is
that Stan's going to need gradients with respect to all the parameters
and expects a templated functor it can automatically differentiate,
which means something implemented with our library of functions where
the partials are known.

If you define the functions in Stan, this'll all happen automatically,
by the way.

- Bob

> On Apr 18, 2015, at 2:22 AM, Dennis Conway <denni...@gmail.com> wrote:
>
> Hey Daniel, thanks for your help on this.
>
> I think it would be quite tricky to split up the imaginary and real parts into two parameters, as there are a few multiplications in my formula which I would need to write in terms of stretchings and rotations. The calculation involves a recursive formula, which I've put at the bottom of this message, and the probability is given by a misfit between the output of the recursive formula and some data that I have. If this is straightforward enough to write up in C++ then I will have a go at it, as I'm curious to see how my model (high dimensionality with highly correlated parameters) performs with the NUTS sampler.
>
> Thanks,
> Dennis
>
> t is real
> c&q are complex valued
> c_{n} = \frac{c_{n-1}\cdot q_{n}+\tanh(tq_{n})}{1+c_{n-1}\cdot q_{n}\cdot \tanh(tq_{n})}\frac{1}{q_{n}}
>

[Dennis Conway]

Thanks for the help Bob, it's becoming much clearer now, and the manual looks well documented for the rest of the problem.

Just to clarify, is the vector 2 x N because there are N steps and we need to keep track of the parameters at each step for the differentiation?

[Bob Carpenter]

If there are N complex variables, you need (2 * N) real variables,
N for the real parts and N for the complex parts.

- Bob

[Dennis Conway]

Oh of course, I'm with you now.

Cheers,

Dennis

[Dennis Conway]

Hello,

I realise that because I have non-linear transformations of my parameters I need to increment logp by the log Jacobian determinant of the complete inverse transform, but I am really confused about where to start on finding this.

I have a Nx1 vector of resistivities, R, which I transform into a complex value Z (which I am representing as a 1x2 vector). What I think I need to do is to rewrite the transformation and take the derivative of each R[i] w.r.t. Z - using complex differentiation - but that doesn't leave me with a square Jacobian.

If anyone could put me on the right track it would be much appreciated!

[Michael Betancourt]

There is not going to be any need for complex differentiation — in fact it

will probably be easier to write out your entire model in terms of the

real-valued components (especially if your likelihood depends on only

the magnitude of the complex numbers).

In any case, you cannot map an N-dimensional distribution into

an (M != N)-dimensional one as probability will not be conserved.

What you have to do is engineer an N->N map and then marginalize.

For example,

p(x, y)

-> p( f(x, y), y) | d (f, y) / d (x, y) | (f, y) 

-> \int dy p( f, y) | d (f, y) / d (x, y) | (f, y)

--
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.
For more options, visit https://groups.google.com/d/optout.

<code.stan>

[Dennis Conway]

Thanks for your help Michael!
 

There is not going to be any need for complex differentiation — in fact it

will probably be easier to write out your entire model in terms of the

real-valued components (especially if your likelihood depends on only

the magnitude of the complex numbers).

My likelihood depends on the distance between two complex numbers, but I could rewrite it in terms of a magnitude and phase difference.
 

In any case, you cannot map an N-dimensional distribution into

an (M != N)-dimensional one as probability will not be conserved.

What you have to do is engineer an N->N map and then marginalize.

Would you be able to explain how I can find an N->N transformation? Can I do something arbitrary like have each entry equal to 1/N of the end value and sum over them afterwards?

Thank you for your help,
Dennis
 

[Michael Betancourt]

There an infinite number of possible transforms — it’s a bit of

an art finding one that is computationally convenient.  Setting

everything to the same output value and then summing may

not be particularly useful as the Jacobian will be nontrivial and

likely prevent an analytic sum.  The strategy that I’ve most often

seen is one I originally described — fill in the (M - N) missing

entries with the last (M - N) original parameters.

[Dennis Conway]

Thanks, that makes a lot of sense! Hopefully I should be able to work through the rest of it now.