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

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.

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)

• Not sure what your first sentence means. Can you clarify? PyMC takes all of the relevant objects and builds a DAG out of them (the graph of the model) using the dependency structure implied by the arguments of each object. For example, if b has a as an argument, then it means a is a parent of b. PyMC then assigns a StepMethod (some MCMC algorithm, often Metropolis-Hastings) to each stochastic in the model, then iterates through each of them at every step of the sampling run. This generates traces of conditional (on the current value of the other parameters) sample for each parameter. – fonnesbeck Apr 28 '13 at 17:37
• (Sorry, made a typo:) Seems like @pm.observed has the same effect as @pm.potential? I understand what potential does, but I'm still a bit confused about how PyMC does its sampling. – roger_ May 2 '13 at 16:18
• @observed has a similar effect as @potential in that they both contribute log-probability terms to the posterior. The latter does not have a value associated with it, and is generally used to constrain the model. The former is for modeling data. Check out the user's guide for more info. – fonnesbeck May 8 '13 at 3:00

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 )

• Thanks, but isn't that assuming that $Z$ is fixed, and not random (of which observations are available)? In reality there's a vector $[z_i]_{i=1}^N$ (and actually it's $Z = X_1 + X_2 + X_3 + ...$), so I'm not sure this would generalize. – roger_ Apr 19 '13 at 4:07
• isn't $Z$ fixed? i.e. you mention you have samples observations $z_i$. – Cam.Davidson.Pilon Apr 19 '13 at 4:12
• Sorry what I meant was I have multiple realizations of it, so it's not just a constant. My goal is infer the parameters of $X$ and $Y$ from $z_i$. – roger_ Apr 19 '13 at 15:29
• but they are still observations. Whereas I have $Z = 1$, you may have $Z = [1,2,3]$, correct? I made some edits to the first approach to reflect this. – Cam.Davidson.Pilon Apr 19 '13 at 15:42
• I made some edits to the 1st approach fyi – Cam.Davidson.Pilon Apr 20 '13 at 1:21

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.