title: Adding measurement errors to pymc model
body: <p>I have the following model in pymc2:</p>
  
  <pre><code>import pymc
  from scipy.stats import gamma
  
  alpha = pymc.Uniform('alpha', 0.01, 2.0)
  scale = pymc.Uniform('scale', 1.0, 4.0)
  
  @pymc.deterministic(plot=False)
  def beta(scale=scale):
      return 1.0 / scale
  
  @pymc.potential
  def p_factor(alpha=alpha, scale=scale, lmin=lmin, n=len(sample)):
      dist = gamma(alpha, loc=0., scale=scale)
      fp = 1.0 - dist.cdf(lmin)
      return -(n+1)*np.log(fp)
  
  obs = pymc.Gamma("obs", alpha=alpha, beta=beta, value=sample, observed=True)
  </code></pre>
  
  <p>The physical background of this model is the <a href="http://en.wikipedia.org/wiki/Luminosity_function_%28astronomy%29" rel="nofollow">luminosity function of galaxies</a> (LF), i.e., the probability of a galaxy having luminosity L. For some types of galaxies, the LF is just a gamma function. The potential accounts for data truncation, as galaxy surveys usually miss a substantial fraction of the targets, particularly those of low luminosity. In this model I miss everything below <code>lmin</code></p>
  
  <p>Details of this method can be found in <a href="http://arxiv.org/abs/0805.2946" rel="nofollow">this paper by Kelly et al</a>.</p>
  
  <p>This model works: I run <code>MAP</code> and <code>MCMC</code> on the model and I can recover the parameters <code>alpha</code> and <code>scale</code> from my simulated data <code>sample</code>, with increased uncertainty as <code>lmin</code> grows.</p>
  
  <p>Now I would like to insert gaussian measurement errors. For simplicity all the data has the same precision. I'm not modifying the potential to include the errors also.</p>
  
  <pre><code>alpha = pymc.Uniform('alpha', 0.01, 2.0)
  scale = pymc.Uniform('scale',1.0, 4.0)
  sig = 0.1
  tau = math.pow(sig, -2.0)  
  
  @pymc.deterministic(plot=False)
  def beta(scale=scale):
      return 1.0 / scale
  
  @pymc.potential
  def p_factor(alpha=alpha, scale=scale, lmin=lmin, n=len(sample)):
      dist = gamma(alpha, loc=0., scale=scale)
      fp = 1.0 - dist.cdf(lmin)
      return -(n+1) * np.log(fp)
  
  dist = pymc.Gamma("dist", alpha=alpha, beta=beta)
  obs = pymc.Normal("obs", mu=dist, tau=tau, value=sample, observed=True)
  </code></pre>
  
  <p>But surely I'm doing something wrong here because this model does not work.
  When I run <code>pymc.MAP</code>on this model I recover the initial values of <code>alpha</code> and <code>scale</code></p>
  
  <pre><code>vals = {'alpha': alpha, 'scale': scale, 'beta': beta, 
     'p_factor': p_factor, 'obs': obs, 'dist': dist}
  M2 = pymc.MAP(vals)
  M2.fit()
  print M2.alpha.value, M2.scale.value
  &gt;&gt;&gt; (array(0.010000000006018368), array(1.000000000833973))
  </code></pre>
  
  <p>When I run <code>pymc.MCMC</code>, <code>alpha</code> and <code>beta</code> are no traced at all.</p>
  
  <pre><code>M = pymc.MCMC(vals)
  M.sample(10000, burn=5000)
  ...
  M.stats()['alpha']
  &gt;&gt;&gt; {'95% HPD interval': array([ 0.01000001,  0.01000502]),
  'mc error': 2.1442678276712383e-07,
  'mean': 0.010001588137798096,
  'n': 5000,
  'quantiles': {2.5: 0.0100000088679046,
  25: 0.010000382359859467,
  50: 0.010001100377476166,
  75: 0.010001668672799679,
  97.5: 0.0100050194240779},
  'standard deviation': 2.189828287191421e-06}
  </code></pre>
  
  <p>again initial values. In fact if I change <code>alpha</code> to start in, say, 0.02, the recovered values of <code>alpha</code> is 0.02.</p>
  
  <p>This is <a href="http://nbviewer.ipython.org/url/guaix.fis.ucm.es/~spr/notebooks/lumfunc/recover%20N1.ipynb?create=1" rel="nofollow">a notebook with the working model plus simulated data</a>.</p>
  
  <p>This is <a href="http://nbviewer.ipython.org/url/guaix.fis.ucm.es/~spr/notebooks/lumfunc/Recover%20N%20with%20errors.ipynb?create=1" rel="nofollow">a notebook with the error model plus simulated data</a>.</p>
  
  <p>Any guidance on making this work would be really appreciated.</p>