PyMC: ¿Cómo puedo definir una función de dos variables estocásticas, no distribución de forma cerrada?

Question

Estoy aprendiendo PyMC y básicamente tengo una variable aleatoria $Z = X + Y$ donde (dicen) $X \sim \mathrm{Normal}(\theta_X)$ $Y \sim \mathrm{Lognormal}(\theta_Y)$ $Z$ no tiene simple y cerrada de forma de la distribución. Ahora tengo observaciones $z_i,\,i=1...N$ $Z$ y quiero inferir $\theta_X$$\theta_Y$. ¿Cuál es la más sencilla forma de hacer esto con PyMC?

Si yo tuviera la distribución de $Z$ disponible, entonces creo que podría hacer:

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

y, a continuación, hacer inferencia, pero no sé DistZ. También es fácil definir la suma como:

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

pero luego pienso que no se puede definir una función determinista.

Yo creo que yo podría hacer algo como:

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

pero no me queda claro en los detalles. Sé el conjunto probabilidad de $\mathcal{L}(z, x)$, pero tenía la esperanza de que no sería necesario.

Yo era capaz de hacer esto con una costumbre muestreador de Gibbs (mediante la articulación de probabilidad), pero estoy buscando una forma más "elegante" con la solución de PyMC.

---

EDIT: se encontró una pregunta similar en los ERRORES de preguntas frecuentes que dice: funciones de variables aleatorias no son compatibles. No sé si esto se aplica a PyMC, y lo que el enfoque estándar.

Answer 1

Usar un enfoque de variable latente, ya que eso es lo que x una y. Sin embargo, no claro que todos los cuatro parámetros sería identificables en este caso. Sería útil si tuvieras alguna información previa de uno o dos de ellos. Aquí está un ejemplo:

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

Creo que hay un par de enfoques.

Primera Aproximación

Que yo sepa, no hay ninguna forma de uso @deterministic o @stochastic (sin la probabilidad). Una forma alternativa es el uso de la potentials de la clase, que es como multiplicando la probabilidad por un factor. En este caso, se debe multiplicar por el pdf de un lognormal dado $Z$ $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)

Aviso de $Z$ es negativo, por lo que este debe hacer $X$ negativo. Y nos encontramos con esto:

Por simetría (desde $Y = Z-X$) de la parte posterior de $Y$ es similar:

Z es un vector de observaciones

Si $Z$ es un vector de observaciones, a continuación, la función potencial que puede ser modificado para que parezca:

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,  )

Para extender a más de dos combinaciones lineales, por ejemplo,$Z = X_1 + X_2 + ... +X_N$, bien para ser continuado.

---

Segundo Método

Un enfoque más específico, es notar que como $X$ es normal, se puede pensar en esta tarea como $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 )

En ejecución de esta segunda versión me dio algunos resultados inestables (se devuelve un puñado increíblemente grandes valores )