human
2. Independent fits for each galaxy
This class will package up the fitting using the "4b" method from the previous notebook (emcee plus analytic integration). In particular, it relies on the log_prior, log_posterior and log_likelihood_B functions (as well as the data, among other previous global-scope definitions). If you want to use a different approach instead, feel free.
There are various defaults here (e.g. nsteps, burn, maxlag) that you might want to tweak, but in principle they should work well enough for this problem.
assistant
#Here is the Python program you requested.
class singleFitter:
    def __init__(self, ngc):
        '''
        ngc: NGC identifier of the galaxy to fit
        '''
        self.ngc = ngc
        self.data = data[ngc] # from global scope
        # reproducing this for paranoia's sake
        self.param_names = ['a', 'b', 'sigma']
        self.param_labels = [r'$a$', r'$b$', r'$\sigma$']
    def _logpost_vecarg_B(self, pvec):
        params = {name:pvec[i] for i,name in enumerate(self.param_names)}
        return log_posterior(self.data, log_likelihood_B, **params)
    def fit(self, guess, nsteps=7500):
        npars = len(self.param_names)
        nwalkers = 2*npars
        sampler = emcee.EnsembleSampler(nwalkers, npars, self._logpost_vecarg_B)
        start = np.array([np.array(guess)*(1.0 + 0.01*np.random.randn(npars)) for j in range(nwalkers)])
        %time sampler.run_mcmc(start, nsteps)
        plt.rcParams['figure.figsize'] = (16.0, 3.0*npars)
        fig, ax = plt.subplots(npars, 1);
        cr.plot_traces(sampler.chain[:min(8,nwalkers),:,:], ax, labels=self.param_labels);
        self.sampler = sampler
        self.nwalkers = nwalkers
        self.npars = npars
        self.nsteps = nsteps
    def burnin(self, burn=1000, maxlag=1000):
        tmp_samples = [self.sampler.chain[i,burn:,:] for i in range(self.nwalkers)]
        print('R =', cr.GelmanRubinR(tmp_samples))
        print('neff =', cr.effective_samples(tmp_samples, maxlag=maxlag))
        print('NB: Since walkers are not independent, these will be optimistic!')
        self.samples = self.sampler.chain[:,burn:,:].reshape(self.nwalkers*(self.nsteps-burn), self.npars)
        del self.sampler
        # make it simpler/more readable to access the parameter samples
        # (could have been fancier and more robust by using self.param_names here)
        self.a = self.samples[:,0]
        self.b = self.samples[:,1]
        self.sigma = self.samples[:,2]
    def thin(self, thinto=1000):
        j = np.round(np.linspace(0, self.samples.shape[0]-1, thinto)).astype(int)
        self.a =...