title: Posterior probability with pymc
body: <p>(This question was originally posted on stats.O. I moved it here because it does relate with <code>pymc</code> and more general matters within it: in fact the main aim is to have a better understanding of how <code>pymc</code> works. If any of the moderators believe it not to be suitable for SO, I'll remove from here.)</p>
  <p>I've been reading the <a href="https://pymc-devs.github.io/pymc/tutorial.html#equation-disaster_model" rel="nofollow noreferrer">pymc tutorial</a> and many other questions both in here and in SO.</p>
  <p>I am trying to understand how to apply <em>Bayes' theorem</em> to compute the posterior probability using certain data. In particular, I have a tuple of independent parameters <img src="https://latex.codecogs.com/gif.latex?$%5Ctheta=(%5Ctheta_1,...,%5Ctheta_n)$" title="$\theta=(\theta_1,...,\theta_n)$" /></p>
  <p><img src="https://latex.codecogs.com/gif.latex?p(%5Ctheta)=%5Cprod_ip(%5Ctheta_1)." title="p(\theta)=\prod_ip(\theta_1)." /></p>
  <p>From the data I'd like to infer the likelihood <img src="https://latex.codecogs.com/gif.latex?$p(E%5Cmid%5Ctheta)$" title="$p(E\mid\theta)$" />, where <img src="https://latex.codecogs.com/gif.latex?E" title="E" /> is a certain event. Then the aim is to compute</p>
  <p><img src="https://latex.codecogs.com/gif.latex?p(%5Ctheta%5Cmid&space;E)=%5Cfrac%7Bp(E%5Cmid%5Ctheta)p(%5Ctheta)%7D%7Bp(E)%7D." title="p(\theta\mid E)=\frac{p(E\mid\theta)p(\theta)}{p(E)}." /></p>
  <p><strong>Some additional comments</strong>:</p>
  <ol>
  <li>This is a sort of unsupervised learning, I know that <img src="https://latex.codecogs.com/gif.latex?E" title="E" /> happened, and I want to find out the parameters <img src="https://latex.codecogs.com/gif.latex?$%5Ctheta$" title="$\theta$" /> that maximise the probability <img src="https://latex.codecogs.com/gif.latex?$p(%5Ctheta&space;|&space;E)$" title="$p(\theta | E)$" />. (*)</li>
  <li>I'd also like to have a parallel procedure where I let <code>pymc</code> compute the likelihood given the data and then <strong>for each tuple of parameters</strong> I want to get the posterior probability.</li>
  </ol>
  <p>In the following I'll assume that <img src="https://latex.codecogs.com/gif.latex?$p(%5Ctheta)=%5Cmathcal%7BU%7D(0,100)%5Ctimes%5Cmathcal%7BN%7D(0,&space;0.0001)$" title="$p(\theta)=\mathcal{U}(0,100)\times\mathcal{N}(0, 0.0001)$" /> and that the likelihood is a multi-dimensional normal distribution with <img src="https://latex.codecogs.com/gif.latex?$%5Cmu=%5Ctheta%5Cmbox%7B&space;and&space;%7D%5Csigma=I_n$" title="$\mu=\theta\mbox{ and }\sigma=I_n$" /> (because of the independence).</p>
  <p>The following is the code I'm using (for simplicity assume there are only two parameters). The code is still under development (I know it can't work!). But I believe it is useful to include it, and then refine it following comments and answers to provide a skeleton for future reference.</p>
  <pre><code>class ObsData(object):
      def __init__(self, params):
          self.theta1 = params[0]
          self.theta2 = params[1]
  
  class Model(object):
      def __init__(self, data):
          # Priors
          self.theta1 = pm.Uniform('theta1', 0, 100)
          self.theta2 = pm.Normal('theta2', 0, 0.0001)
  
          @pm.deterministic
          def model(
              theta1=self.theta1,
              theta2=self.theta2,
          ):
              return (theta1, theta2)
          # Is this the actual likelihood?
          self.likelihood = pm.MvNormal(
              'likelihood',
              mu=model,
              tau=np.identity(2),
              value=data,  # is it correct to put the data here?
              observed=True
          )
  
  def mcmc(observed_data):
      data = ObsData(observed_data)
      pymc_obj = Model(data)
      model = pymc.MCMC(pymc_obj)
      model.sample(10000, verbose=0) # Does this line compute the likelihood and the normalisation factor? 
      # How do I get the posterior distribution?
  </code></pre>
  <p>The questions are the following:</p>
  <ol>
  <li>Does <code>self.likelihood</code> represent the Bayesian likelihood?</li>
  <li>How do I make use of the data? (I suspect <code>value=data</code> is incorrect..)</li>
  <li>Does <code>.sample()</code> actually compute the posterior probability?</li>
  <li>How do I get information from the posterior?</li>
  <li>(*)Should I include anything that relates to the fact that <img src="https://latex.codecogs.com/gif.latex?E" title="E" /> happened at some point?</li>
  <li><strong>As general question</strong>: is it any way to use <code>pymc</code> to only compute the likelihood given the data and the priors?</li>
  </ol>
  <p>Any comments, and also reference to other question or tutorial are welcome!</p>