def fit_mcmc(self,nwalkers=<NUM_LIT:200>,nburn=<NUM_LIT:100>,niter=<NUM_LIT:200>,<EOL>p0=None,initial_burn=None,<EOL>ninitial=<NUM_LIT:100>, loglike_kwargs=None,<EOL>**kwargs):
<EOL>if self._samples is not None:<EOL><INDENT>self._samples = None<EOL><DEDENT>if self.fit_for_distance:<EOL><INDENT>npars = <NUM_LIT:6><EOL>if initial_burn is None:<EOL><INDENT>initial_burn = True<EOL><DEDENT><DEDENT>else:<EOL><INDENT>if initial_burn is None:<EOL><INDENT>initial_burn = False<EOL><DEDENT>npars = <NUM_LIT:4><EOL><DEDENT>if p0 is None:<EOL><INDENT>mA_0,age0,feh0 = self.ic.random_points(nwalkers)<EOL>mB_0,foo1,foo2 = self.ic.random_points(nwalkers)<EOL>mA_fixed = np.maximum(mA_0,mB_0)<EOL>mB_fixed = np.minimum(mA_0,mB_0)<EOL>mA_0, mB_0 = (mA_fixed, mB_fixed)<EOL>d0 = <NUM_LIT:10>**(rand.uniform(<NUM_LIT:0>,np.log10(self.max_distance),size=nwalkers))<EOL>AV0 = rand.uniform(<NUM_LIT:0>,self.maxAV,size=nwalkers)<EOL>if self.fit_for_distance:<EOL><INDENT>p0 = np.array([mA_0,mB_0,age0,feh0,d0,AV0]).T<EOL><DEDENT>else:<EOL><INDENT>p0 = np.array([mA_0,mB_0,age0,feh0]).T<EOL><DEDENT>if initial_burn:<EOL><INDENT>sampler = emcee.EnsembleSampler(nwalkers,npars,self.lnpost,<EOL>**kwargs)<EOL>pos, prob, state = sampler.run_mcmc(p0, ninitial) <EOL>wokinds = np.where((sampler.naccepted/ninitial > <NUM_LIT>) &<EOL>(sampler.naccepted/ninitial < <NUM_LIT>))[<NUM_LIT:0>]<EOL>i=<NUM_LIT:1><EOL>while len(wokinds)==<NUM_LIT:0>:<EOL><INDENT>thresh = <NUM_LIT> - i*<NUM_LIT><EOL>if thresh < <NUM_LIT:0>:<EOL><INDENT>raise RuntimeError('<STR_LIT>')<EOL><DEDENT>wokinds = np.where((sampler.naccepted/ninitial > thresh) &<EOL>(sampler.naccepted/ninitial < <NUM_LIT>))[<NUM_LIT:0>]<EOL>i += <NUM_LIT:1><EOL><DEDENT>inds = rand.randint(len(wokinds),size=nwalkers)<EOL>p0 = sampler.chain[wokinds[inds],:,:].mean(axis=<NUM_LIT:1>) <EOL>p0 *= (<NUM_LIT:1> + rand.normal(size=p0.shape)*<NUM_LIT>)<EOL><DEDENT><DEDENT>else:<EOL><INDENT>p0 = np.array(p0)<EOL>p0 = rand.normal(size=(nwalkers,npars))*<NUM_LIT> + p0.T[None,:]<EOL>if self.fit_for_distance:<EOL><INDENT>p0[:,<NUM_LIT:4>] *= (<NUM_LIT:1> + rand.normal(size=nwalkers)*<NUM_LIT:0.5>)<EOL><DEDENT><DEDENT>sampler = emcee.EnsembleSampler(nwalkers,npars,self.lnpost)<EOL>pos, prob, state = sampler.run_mcmc(p0, nburn)<EOL>sampler.reset()<EOL>sampler.run_mcmc(pos, niter, rstate0=state)<EOL>self._sampler = sampler<EOL>return sampler<EOL>
Fits stellar model using MCMC.

        See :func:`StarModel.fit_mcmc`