possible pymc.Gamma issue?

gistfile1.txt

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Gamma Model"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pymc as mc"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# simple gamma model\n",
      "with mc.Model() as gamma_model:\n",
      "    tau = mc.Gamma('tau', alpha=1.0, beta=1.0)\n",
      "    start = mc.find_MAP()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "ValueError",
       "evalue": "Optimization error: max, logp or dlogp at max have bad values. max: array([-0.00048828]) logp: array(-inf) dlogp: array([ 0.])Check that 1) you don't have hierarchical parameters, these will lead to points with infinite density. 2) your distribution logp's are properly specified.",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-2-b98bea2619ec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mmc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mModel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mgamma_model\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0mtau\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGamma\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'tau'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0malpha\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbeta\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m     \u001b[0mstart\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfind_MAP\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/pymc/tuning/starting.pyc\u001b[0m in \u001b[0;36mfind_MAP\u001b[0;34m(start, vars, fmin, return_raw, disp, model, *args, **kwargs)\u001b[0m\n\u001b[1;32m     63\u001b[0m             not allfinite(dlogp(mx))):\n\u001b[1;32m     64\u001b[0m             raise ValueError(\"Optimization error: max, logp or dlogp at max have bad values. max: \" + repr(mx) + \" logp: \" + repr(logp(mx)) + \" dlogp: \" + repr(dlogp(mx)) +\n\u001b[0;32m---> 65\u001b[0;31m                              \"Check that 1) you don't have hierarchical parameters, these will lead to points with infinite density. 2) your distribution logp's are properly specified.\")\n\u001b[0m\u001b[1;32m     66\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     67\u001b[0m     \u001b[0mmx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbij\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;31mValueError\u001b[0m: Optimization error: max, logp or dlogp at max have bad values. max: array([-0.00048828]) logp: array(-inf) dlogp: array([ 0.])Check that 1) you don't have hierarchical parameters, these will lead to points with infinite density. 2) your distribution logp's are properly specified."
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Find Starting Values Step by Step"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Use code from https://github.com/pymc-devs/pymc/blob/pymc3/pymc/tuning/starting.py to inspect what's going on in fmin_bfgs"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# just get variables setup\n",
      "model = gamma_model\n",
      "start = { 'tau': 1.0 } # set a valid starting value for tau\n",
      "start = mc.Point(start, model=gamma_model)\n",
      "vars = model.vars"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# setup logp/dlogp\n",
      "bij = mc.DictToArrayBijection(mc.ArrayOrdering(vars), start)\n",
      "logp = bij.mapf(model.logpc)\n",
      "dlogp = bij.mapf(model.dlogpc(vars))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# some helpfer functions\n",
      "def nan_to_high(x):\n",
      "    return mc.core.np.where(mc.core.np.isfinite(x), x, 1.0e100)\n",
      "def logp_o(point):\n",
      "    return nan_to_high(-logp(point))\n",
      "def grad_logp_o(point):        \n",
      "    return mc.core.np.nan_to_num(-dlogp(point))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# run the optimization\n",
      "from scipy.optimize import fmin_bfgs\n",
      "r = fmin_bfgs(logp_o, bij.map(start), fprime=grad_logp_o)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Optimization terminated successfully.\n",
        "         Current function value: 10000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815104.000000\n",
        "         Iterations: 1\n",
        "         Function evaluations: 17\n",
        "         Gradient evaluations: 6\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# inspect values\n",
      "print('tau: {tau}\\nlogp: {logp}\\ndlogp: {dlogp}'.format(tau=r, logp=logp(r), dlogp=dlogp(r)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "tau: [-0.00048828]\n",
        "logp: -inf\n",
        "dlogp: [ 0.]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Different Priors"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If we use a Gamma(10, 1) it works..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "with mc.Model() as gamma_model_2:\n",
      "    tau = mc.Gamma('tau', alpha=10.0, beta=1.0)\n",
      "    start = mc.find_MAP()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    }
   ],
   "metadata": {}
  }
 ]
}