prappleizer / prappleizer.github.io

This repo contains my introductory python textbook for astronomy students, which covers the basics of learning the language with an emphasis on astronomical applications.
http://prappleizer.github.io/textbook.pdf
55 stars 22 forks source link

Fitting a multivariate curve using Bayesian inference #5

Closed haricash closed 2 years ago

haricash commented 2 years ago

Hello! I am learning how to use emcee and I came across your tutorial where you fit a single variable model to the data. I wanted to try out an example with multiple variables and so did this :

import numpy as np

x = np.linspace(0,10,100)
y = np.linspace(0,10,100)
z = np.linspace(0,10,100)
noise = np.random.rand(100)
f = 1.5*x + 3.0*y + 4.5*z + noise + 6 

import sys
sys.path.append("PATH/TO/EMCEE")

import emcee

def model(params, var):
    a, b, c, d = params
    x, y, z = var
    return a*x + b*y + c*z + d

def log_likelihood(params, var, f, ferr):
    return ((f - model(params, var))/ferr)**2

def log_prior(params):
    a, b, c, d = params
    if (0 < a < 5) and (0 < b < 6) and (0 < c < 7) and (0 < d < 8):
        return 0.0
    else:
        return -np.inf

def log_prob(params, var, f, ferr):
    lp = log_prior(params)
    if not np.isfinite(lp):
        return -np.inf
    else:
        return lp + log_likelihood(params, var, f, ferr)

var = x, y, z
data = (var, f, noise) 
nwalkers = 32
niter = 500
initial = np.array([1, 2, 3, 4]) 
ndim = len(initial)
p0 = [np.array(initial) + 1e-7 * np.random.randn(ndim) for i in range(nwalkers)] 

def main(p0,nwalkers,niter,ndim,log_prob,data):
    sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob, args=data)

    print("Running burn-in...")
    p0, _, _ = sampler.run_mcmc(p0, 100)
    sampler.reset()

    print("Running production...")
    pos, prob, state = sampler.run_mcmc(p0, niter)

    print("Run Complete!")

    return sampler, pos, prob, state

sampler, pos, prob, state = main(p0,nwalkers,niter,ndim,log_prob,data)

and it is showing this error :

Running burn-in...

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
~/PATH/TO/lib/python3.9/site-packages/emcee/ensemble.py in _get_lnprob(self, pos)
    384         try:
--> 385             lnprob = np.array([float(l[0]) for l in results])
    386             blob = [l[1] for l in results]

~/PATH/TO/lib/python3.9/site-packages/emcee/ensemble.py in <listcomp>(.0)
    384         try:
--> 385             lnprob = np.array([float(l[0]) for l in results])
    386             blob = [l[1] for l in results]

TypeError: 'float' object is not subscriptable

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
<ipython-input-73-9dd151c3debc> in <module>
     58     return sampler, pos, prob, state
     59 
---> 60 sampler, pos, prob, state = main(p0,nwalkers,niter,ndim,log_prob,data)

<ipython-input-73-9dd151c3debc> in main(p0, nwalkers, niter, ndim, log_prob, data)
     48 
     49     print("Running burn-in...")
---> 50     p0, _, _ = sampler.run_mcmc(p0, 100)
     51     sampler.reset()
     52 

~/PATH/TO/lib/python3.9/site-packages/emcee/sampler.py in run_mcmc(self, pos0, N, rstate0, lnprob0, **kwargs)
    169                 rstate0 = self._last_run_mcmc_result[2]
    170 
--> 171         for results in self.sample(pos0, lnprob0, rstate0, iterations=N,
    172                                    **kwargs):
    173             pass

~/PATH/TO/lib/python3.9/site-packages/emcee/ensemble.py in sample(self, p0, lnprob0, rstate0, blobs0, iterations, thin, storechain, mh_proposal)
    256                 first, second = slice(halfk), slice(halfk, self.k)
    257                 for S0, S1 in [(first, second), (second, first)]:
--> 258                     q, newlnp, acc, blob = self._propose_stretch(p[S0], p[S1],
    259                                                                  lnprob[S0])
    260                     if np.any(acc):

~/PATH/TO/lib/python3.9/site-packages/emcee/ensemble.py in _propose_stretch(self, p0, p1, lnprob0)
    330         # Calculate the proposed positions and the log-probability there.
    331         q = c[rint] - zz[:, np.newaxis] * (c[rint] - s)
--> 332         newlnprob, blob = self._get_lnprob(q)
    333 
    334         # Decide whether or not the proposals should be accepted.

~/PATH/TO/lib/python3.9/site-packages/emcee/ensemble.py in _get_lnprob(self, pos)
    386             blob = [l[1] for l in results]
    387         except (IndexError, TypeError):
--> 388             lnprob = np.array([float(l) for l in results])
    389             blob = None
    390 

~/PATH/TO/lib/python3.9/site-packages/emcee/ensemble.py in <listcomp>(.0)
    386             blob = [l[1] for l in results]
    387         except (IndexError, TypeError):
--> 388             lnprob = np.array([float(l) for l in results])
    389             blob = None
    390 

TypeError: only size-1 arrays can be converted to Python scalars

I think I kinda understand the issue here - there is something wrong with the way I am passing the variables into the data holder (that goes into args). I was hoping you could help me understand what is going on and how to rectify it.