How do I define mulitple RVs with a joint liklihood?

Questions
1

I’m quite new to pymc3 and have difficulties understanding what happens if I define multiple observed RVs in my model.

I have the following setup:

  • I’m observing a set of complex numbers called visibilities V_i. Every V_i is drawn from a normal distribution \mathcal{N}(\tilde{V}_i, \sigma_i), where \tilde{V}_i is the visibility obtained from a model \tilde{V}_i = FFT(B(\Theta))_i
  • \Theta = \theta_A is the parameter I want to infere
  • I have some uniform prior P(\theta_A)
  • My likelihood is \mathcal{L}(V\vert \theta_A) = \prod_{i=0}^N{\mathcal{N}(\tilde{V}_i(\theta_A), \sigma_i)}. For simplicity let \sigma_i=1.

My question is: How am I supposed to implement this likelihood inside my pymc3 model?

My first solution is this, but I’m not sure if this is correct.

amplitude_model = pymc3.Model()
with amplitude_model:
    amplitude = pymc3.Uniform('amplitude', 0.1, 100) # uniform prior for A
    source_A = pointSourceALM(128,amplitude,32,64) # calculate source for A
    vis_A = theano_RIME_simple(source_A, uvstack, 128) # calculate visibilites using RIME (simple FFT)
    vis_A_real = vis_A[:,:,0].flatten() # flatten both arrays since we don't need positional information & split real/imag
    vis_A_imag = vis_A[:,:,1].flatten()
    vis_flat = vis.flatten() # observed visibilities
    cov = np.array([[1, 0], [0, 1]])
    
    # Generate one normal distirubtion for each visibility
    for i, v in enumerate(zip(vis_flat.real,vis_flat.imag)):
        if v == (0,0): # skip visibilites that have not been measured
            continue
        pymc3.MvNormal('vis'+str(i), mu=[vis_A_real[i],vis_A_imag[i]], cov=cov,observed=v)

I’m not really sure what happens inside of pymc3 when I define multiple observed RVs (vis0…visN).

Thanks in advance for any help.

1 Like
2

Welcome!

I don’t think anything particularly surprising happens when you have multiple, separate observations. During sampling, the log posterior of the sampled parameter values is just the sum of all of the relevant log posteriors.

That being said, you probably want to avoid loops whenever possible because they can’t be optimized under the hood. So it would be ideal to vectorize your loop if at all possible. I didn’t dig into your code too much, but it seems like the parametric components of your likelihood ( vis_A) and your data (vis) may already be in a state that they can be reshaped to allow for straightforward vector operations.

3 Likes
3

Thank very much for your input. You’re right I can actually rewrite it as two multivariate normal distributions for the real and imaginary part I think.

\mathcal{L}(V\vert \theta_A) = \frac{1}{\sqrt{(2\pi)^N}\det{\Sigma}}e^{-\frac{1}{2}(\vec{V}-\vec{\tilde{V}})^T\Sigma^{-1}(\vec{V}-\vec{\tilde{V}})}

with \Sigma being an diagonal matrix (Cov(V_i, V_j)=0)

1 Like