title: How to define a custom prior in PyMC3
body: <p>I would like to know if it is possible to define a custom prior in PyMC3 (and how to do it).  From <a href="https://stackoverflow.com/questions/23198247/custom-priors-in-pymc">here</a> it seems that in PyMC2 is relatively easy to do (without the need to modified the source code), but in PyMC3 is not that easy (or I am not understanding something).
  I am trying to replicate a prior from the book "Doing Bayesian Data Analysis", that is implemented in BUGS:</p>
  
  <pre><code>model {
  # Likelihood. Each flip is Bernoulli.
  for ( i in 1 : N1 ) { y1[i]  ̃ dbern( theta1 ) }
  for ( i in 1 : N2 ) { y2[i]  ̃ dbern( theta2 ) }
  # Prior. Curved scallo not ps!
  x  ̃ dunif(0,1)
  y  ̃ dunif(0,1)
  N &lt;- 4
  xt &lt;- sin( 2*3.141593*N * x ) / (2*3.141593*N) + x
  yt &lt;- 3 * y + (1/3)
  xtt &lt;- pow( xt , yt )
  theta1 &lt;- xtt
  theta2 &lt;- y
  }
  </code></pre>
  
  <p>The prior does not make too much sense, it just an example of how to define a custom prior and the versatility of BUGS.</p>
  
  <p>My attempt to implement the above custom prior is:</p>
  
  <pre><code>from __future__ import division
  import numpy as np
  import pymc as pm
  from pymc import Continuous
  from theano.tensor import sin, log
  
  # Generate the data
  y1 = np.array([1, 1, 1, 1, 1, 0, 0])  # 5 heads and 2 tails
  y2 = np.array([1, 1, 0, 0, 0, 0, 0])  # 2 heads and 5 tails
  
  class Custom_prior(Continuous): 
  """
  custom prior
  """
      def __init__(self, y, *args, **kwargs):
          super(Custom_prior, self).__init__(*args, **kwargs)
          self.y = y
          self.N = 4
          self.mean = 0.625  # FIXME
      def logp(self, value):
          N = self.N
          y = self.y
          return -log((sin(2*3.141593*N * value)
                       / (2*3.141593*N) + value)**(3 * y + (1/3)))
  
  with pm.Model() as model:
      theta2 = pm.Uniform('theta2', 0, 1)  # prior
      theta1 = Custom_prior('theta1', theta2)  # prior
      # define the likelihood
      y1 = pm.Bernoulli('y1', p=theta1, observed=y1)
      y2 = pm.Bernoulli('y2', p=theta2, observed=y2)
      # Generate a MCMC chain
      start = pm.find_MAP()  # Find starting value by optimization
      trace = pm.sample(5000, pm.NUTS(), progressbar=False)
  </code></pre>
  
  <p>EDIT</p>
  
  <p>Following the answer of <a href="https://stackoverflow.com/a/24739407/815542">chris-fonnesbeck</a></p>
  
  <p>I think I need something like:</p>
  
  <pre><code>with pm.Model() as model:
      theta2 = pm.Uniform('theta2', 0, 1)  # prior
      N = 4
      theta1 = pm.DensityDist('theta1', lambda value: -log((sin(2*3.141593*N * value)
                         / (2*3.141593*N) + value)**(3 * theta2 + (1/3))))
      # define the likelihood
      y1 = pm.Bernoulli('y1', p=theta1, observed=y1)
      y2 = pm.Bernoulli('y2', p=theta2, observed=y2)
  
      # Generate a MCMC chain
      start = pm.find_MAP()  # Find starting value by optimization
      trace = pm.sample(10000, pm.NUTS(), progressbar=False) # Use NUTS sampling
  </code></pre>
  
  <p>The only problem is that I get the same value for all the posterior samples of theta1 and theta2, I guess there is some problem with my custom prior or the combination of priors and likelihoods. A successful definition of a custom prior can be found in this <a href="https://github.com/pymc-devs/pymc/blob/master/pymc/examples/custom_dists.py" rel="nofollow noreferrer">example</a></p>