PyMC: how can I define a function of two stochastic variables, with no closed-form distribution?

Question

I'm learning PyMC and basically I have a random variable $Z = X + Y$ where (say) $X \sim \mathrm{Normal}(\theta_X)$ and $Y \sim \mathrm{Lognormal}(\theta_Y)$ and $Z$ has no simple closed-form distribution. Now I have observations $z_i,\,i=1...N$ of $Z$ and I want to infer $\theta_X$ and $\theta_Y$. What's the most straight-forward way of doing this with PyMC?

If I had the distribution of $Z$ available, then I think I could do:

Z = DistZ('Z', param_x=theta_x, param_y=theta_y, value=z, observed=True)

and then do inference, but I don't know DistZ. It's also easy to define the sum as:

@pymc.deterministic
def z_sum(x=Y, Y=y):
    return x + y

but then I don't think I can define an observed deterministic function.

I think I could do something like:

@pymc.stochastic(observed=True)
def z_sum(value=z, x=X, y=Y):
    def logp(z, x):
        # return log-likelihood

but I'm not clear on the details. I do know the joint likelihood $\mathcal{L}(z, x)$, but I was hoping it wouldn't be needed.

I was able to do this with a custom Gibbs sampler (using the joint likelihood), but I'm looking for a more "elegant" solution with PyMC.

---

EDIT: found a similar question in the BUGS FAQ that says functions of random variables aren't supported. Not sure if that applies to PyMC, and what the standard approach is.

Answer 1

I would use a latent variable approach, since that's what x an y are. However, its not clear that all four parameters would be identifiable in this case. It would be helpful if you had some prior information for one or two of them. Here's an example:

import pymc as pm

# Priors
mu_x = pm.Normal('mu_x', 0, 0.001, value=0)
sigma_x = pm.Uniform('sigma_x', 0, 100, value=1)
tau_x = sigma_x**-2

mu_y = pm.Normal('mu_y', 0, 0.001, value=1)
sigma_y = pm.Uniform('sigma_y', 0, 100, value=1)
tau_y = sigma_y**-2

# Latent variable
y = pm.Lognormal('y', mu_y, sigma_y, size=len(z_data))

@pm.observed
def Z(value=z_data, mu=mu_x, tau=tau_x, y=y):
    # Likelihood for x (also latent, but fixed given y and z)
    return pm.normal_like(value-y, mu, tau)

Answer 2

I think there are a few approaches here.

First Approach

As far as I know, there is no way to use @deterministic or @stochastic (without the likelihood). An alternative way is to use the potentials class, which is like multiplying your likelihood by a factor. In this case, we should multiply by the pdf of a lognormal given $Z$ and $X$.

import pymc as mc

z = -1.

#instead of 0 and 1, unknowns can be put here. For example:
# mc.Normal( "x", unknown_mu, unknown_std ).

X = mc.Normal( "x", 0, 1, value = -2. ) 


@mc.potential
def Y( x =X, z = z): #similar to my comment above, you can place unknowns here in place of 1, 0.2. 
  return mc.lognormal_like( z-x, 1, 0.2,  )


mcmc = mc.MCMC( [X] )
mcmc.sample(20000,5000)

Notice $Z$ is negative, so this must make $X$ negative too. And we observe this:

By symmetry (since $Y = Z-X$) the posterior of $Y$ is similar:

Z is a vector of observations

If $Z$ is a vector of observations, then the potential function can be modified to look like:

z = [2,3,4]

#...
X = mc.Normal( "x", 0, 1, value = -2., size = 3 ) 

@mc.potential
def Y( x =X, Z = Z):
  return mc.lognormal_like( Z, 1, 0.2,  )

To extend to more than two linear combinations, eg $Z = X_1 + X_2 + ... +X_N$, well to be continued.

---

Second Approach

A more specific approach is to notice that as $X$ is normal, we can think of this task as $Z = Y + \text{noise}$:

import pymc as mc

Z = -1

Y = mc.Lognormal( "y", 1, 0.2 )

obs = mc.Normal( "obs", 0, 1, value = Z, observed = True )


mcmc = mc.MCMC( [Y, obs] )

mcmc.sample( 20000,5000 )

Running this second version did give me some unstable results (was returned a handful incredibly large values )

Answer 3

Conceptually, to do Bayesian inference, one has to serve a conditional likelihood function, with specific pdf. In your case, you have to provide the actual joint distribution of (X,Y). Perhaps potential function can help but the idea is the same, the final entry point to do MCMC should be a specific log pdf.

If X and Y are independent, then the sum of their pdfs is the convolution of the two pdfs. so Z ~ Normal(thetax)*logNormal(thetay). Maybe the convolution integral can be evaluated.