title: pymc 3.0 Predictive Posterior Distribution
body: <p>I'm converting a very simple example from pymc 2.3 to pymc 3.0, and can't seem to figure out how to sample (or get the MAP) from the predictive posterior distribution. Following the suggestion in the docs (7.3 Model checking and diagonostics: Goodness of Fit), I can sample from this distribution using pymc 2.3 by adding a unobserved stochastic. Here is a <a href="http://nbviewer.ipython.org/url/cs.unm.edu/~bedwards/SimplePYMC2.3.ipynb" rel="nofollow">link to the notebook</a>. All seems to work well.</p>
  
  <p>However when I try to do this in <a href="http://nbviewer.ipython.org/url/cs.unm.edu/~bedwards/SimplePYMC3.0.ipynb" rel="nofollow">pymc 3.0</a>, two strange things happen. </p>
  
  <ol>
  <li>The MAP value is not even close, as if the unobserved stochastic is affecting the minimization, and the MAP value for the unobserved stochastic is wrong?</li>
  <li>The NUTS sampler doesn't alter the value of the unobserved stochastic, and so the trace is simply a single value 10.</li>
  </ol>
  
  <p>Of course this value would make sense if the observed stochastic wasn't present. How has finding MAP values for and sampling from the predictive posterior changed in pymc 3.0?</p>
  
  <p>UPDATE:</p>
  
  <p>Here is a minimal example that illustrates what I was doing wrong:</p>
  
  <pre><code>import pymc as mc
  with mc.Model() as model:
      p = mc.Beta('p',2,2)
      surv_sim = mc.Binomial('surv_sim',n=20,p=p)
      surv = mc.Binomial('surv',n=20,p=p,observed=15)
  
  with model:
      step = mc.step_methods.HamiltonianMC(vars=model.vars) #Again must specificy
                                                            #model.vars or else
                                                            #only continuous values will
                                                            #be sampled
  
  with model:
      trace = mc.sample(100000,step)
  
  mc.traceplot(trace);
  
  with model:
      map_est = mc.find_MAP(vars=model.vars) #Must Specify model.vars, or else
                                             #only continuous stochastics will
                                             #be fit, however this fit will be 
                                             #horrible
  </code></pre>
  
  <p>I think I have partially answered this myself. First the sampling. The NUTS sampler is not equipped to handle discrete variables as far as I can tell. However the Metropolis sampler definitely handles discrete (resolved in issue #235) and HamiltonianMC sampler also seems to handle discrete variables.</p>
  
  <p>MAP estimates of the discrete stochastics are destined to be bad however. Even if you specify that all variables be used (not the default). Because discrete stochastics log probability functions return the log probability of the floor of any discrete stochastic, whichever scipy minimization function is used will be caught in a local optimum. The minimization will only optimize the continuous stochastics while keeping the discrete ones fixed as none of the small steps in the discrete variables results in an improvement of log probability. I am not sure if this is a bug or simply a fundamental limitation of finding MAP estimates when you have both discrete and continuous variables.</p>
  
  <p>UPDATE 2:</p>
  
  <p>Sampling works even better if you use multiple steps. For the above model you can use.</p>
  
  <pre><code>with model:
      step1 = mc.step_methods.NUTS(vars=[p])
      step2 = mc.step_methods.Metropolis(vars=[surv_sim])
  
  with model:
      trace = mc.sample(100000,[step1,step2])
  </code></pre>