Search code examples
pythonstandardsmeanpymcdeviation

PyMC MCMC analysis - .summary() "


I'm new to PyMC and am having a little trouble getting out parameters related to my prior, for example, the mean and standard deviation.

I describe my model in a file called 'model.py' like so:

    import pymc
    import numpy

    #constants
    r_div=numpy.loadtxt("r_div", comments="#", unpack=False)
    map=numpy.loadtxt("map", comments="#", delimiter=",", unpack=False)
    M_star=3*10^6;

    #data
    n=numpy.loadtxt("n")

    #priors
    alpha_0=pymc.Uniform('alpha_0end',-10,10, value=0)
    logA_0=pymc.Uniform('logA_0end',-10,10,value=-6.1246)

    #model

    @pymc.deterministic(plot=False)

    def r(logA_0=logA_0,alpha_0=alpha_0,M_star=M_star,r_div=r_div):

        r=r_div*numpy.exp(logA_0)*((numpy.exp(map[:,1])/M_star)**(alpha_0))

        return r

            #likelihood
            Distribution=pymc.Poisson('Distribution',mu=r,value=n,observed=True)

And then I use the following script in ipython to run the MCMC chain:

            import pymc
            import model
            M=pymc.MCMC(model)
            M.sample(100000, burn=10000)
            M.summary()

Everything seems to work until the final command M.summary() where I get the error:


AttributeError Traceback (most recent call last) in () ----> 1 M.summary()

AttributeError: 'MCMC' object has no attribute 'summary'

I am confident the chain has run successfully as the command M.trace('alpha_0end')[:] shows there are chain elements there but I can't get out any information out about the prior, such as a mean or standard deviation. I've tried different permutations of the summary command. For example: M.alpha_0end.summary()

It would be helpful to know if there is an easy way to get out the standard deviation and means of the priors.


Solution

  • I can't run your code, but M.summary() works for me in this minimal example:

    In [1]: import pymc as pm
    
    In [2]: pm.__version__
    Out[2]: '2.3.2'
    
    In [3]: X = pm.Normal('X', 0, 1)
    
    In [4]: M = pm.MCMC(dict(X=X))
    
    In [5]: M.sample(100000, burn=10000)
     [-----------------100%-----------------] 100000 of 100000 complete in 5.0 sec
    In [6]: M.summary()
    
        X:
    
            Mean             SD               MC Error        95% HPD interval
            ------------------------------------------------------------------
            -0.0             1.003            0.003            [-1.897  2.026]
    
    
            Posterior quantiles:
    
            2.5             25              50              75             97.5
             |---------------|===============|===============|---------------|
            -1.961           -0.673          0.0            0.675         1.964