Prioridades personalizadas en PyMC

Digamos que quiero poner una costumbre previa en dos variables a y b en PyMC, por ejemplo:

p(a,b)∝(a+b)^(−5/2)

(Para la motivación detrás de esta elección de prior, vea esta respuesta )

¿Se puede hacer esto en PyMC? ¿Si es así, cómo?

A modo de ejemplo, me gustaría definir tales previas en a y b en el modelo a continuación.

 import pymc as pm # ... # Code that defines the prior: p(a,b)∝(a+b)^(−5/2) # ... theta = pm.Beta("prior", alpha=a, beta=b) # Binomials that share a common prior bins = dict() for i in xrange(N_cities): bins[i] = pm.Binomial('bin_{}'.format(i), p=theta,n=N_trials[i], value=N_yes[i], observed=True) mcmc = pm.MCMC([bins, ps]) 

Actualizar

Siguiendo el consejo de John Salvatier, bash lo siguiente (tenga en cuenta que estoy en PyMC2, aunque me encantaría cambiar a PyMC3), pero mis preguntas son:

  1. ¿Qué debo importar para poder heredar correctamente de Continuous ?
  2. En PyMC2, ¿todavía necesito atenerme a la notación Theano?
  3. Finalmente, ¿cómo puedo decir más tarde a mi distribución Beta que alpha y beta tienen un previo de esta distribución multivariable?

    import pymc.Multivariate.Continuous

    clase CustomPrior (Continuo): “” “p (a, b) ∝ (a + b) ^ (- 5/2)

     :Parameters: None :Support: 2 positive floats (parameters to a Beta distribution) """ def __init__(self, mu, tau, *args, **kwargs): super(CustomPrior, self).__init__(*args, **kwargs) def logp(self, a,b): return np.log(math.power(a+b),-5./2) 

En PyMC2, el truco es juntar los parámetros a y b :

 # Code that defines the prior: p(a,b)∝(a+b)^(−5/2) @pm.stochastic def ab(power=-2.5, value=[1,1]): if np.any(value <= 0): return -np.Inf return power * np.log(value[0]+value[1]) a = ab[0] b = ab[1] 

Este cuaderno tiene un ejemplo completo.

¡Sip! Es bastante posible, y de hecho bastante sencillo.

Si está en PyMC 2, consulte la documentación sobre la creación de variables estocásticas .

 @pymc.stochastic(dtype=int) def switchpoint(value=1900, t_l=1851, t_h=1962): """The switchpoint for the rate of disaster occurrence.""" if value > t_h or value < t_l: # Invalid values return -np.inf else: # Uniform log-likelihood return -np.log(t_h - t_l + 1) 

Si estás en PyMC 3, echa un vistazo a multivariate.py . Tenga en cuenta que los valores pasados ​​a init y logp son todas variables teano, no matrices numpy. ¿Es eso suficiente para empezar?

Por ejemplo, esta es la distribución normal multivariable.

 class MvNormal(Continuous): """ Multivariate normal :Parameters: mu : vector of means tau : precision matrix .. math:: f(x \mid \pi, T) = \frac{|T|^{1/2}}{(2\pi)^{1/2}} \exp\left\{ -\frac{1}{2} (x-\mu)^{\prime}T(x-\mu) \right\} :Support: 2 array of floats """ def __init__(self, mu, tau, *args, **kwargs): super(MvNormal, self).__init__(*args, **kwargs) self.mean = self.median = self.mode = self.mu = mu self.tau = tau def logp(self, value): mu = self.mu tau = self.tau delta = value - mu k = tau.shape[0] return 1/2. * (-k * log(2*pi) + log(det(tau)) - dot(delta.T, dot(tau, delta)))