<jupyter_start><jupyter_text>Fitting MDF modelsI use `dynesty` https://arxiv.org/abs/1904.02180 to do dynamic nested sampling to infer the posterior by default.This is much slower than MCMC like emcee, but pretty foolproof. It also calculates the Bayesian evidence.It's a pretty neat algorithm. Basically, you start with random points from the prior (using a "prior transform"), calculate the likelihood for every point, then iteratively throw out the lowest likelihood point (and add it to your samples, with a weight based on the likelihood) and sample one with strictly higher likelihood from the prior.Yuan-Sen Ting describes this as "ruling out prior volume" rather than "finding the peak and sampling" as a usual MCMC algorithm would work. It costs about 10x as much as MCMC, but is good at finding multiple modes.It's probably overkill for these simple two-parameter models, but I just learned how it works and like it, so you get to suffer/benefit.The code `mdfmodels.fit` has the loglikelihoods for the models if you want to write your own (faster?) fitter.The upshot is that unlike MCMC, if you set a very large prior, it will take a long time to find the answer.  Simple leaky box example Here we will give a simple example of fitting a leaky box model to randomly generated data.Note: data errors are *not* accounted for right now.<jupyter_code># Generate fake data
np.random.seed(230589)
N = 100
p = 0.05
data = utils.sample_leaky_box(N, p)
plt.hist(data, bins='auto');
plt.xlabel("[Fe/H]")
# this step takes a minute or two because dynesty is kinda slow. There's ways to optimize that I haven't found yet
with MultiPool(4) as pool:
    output = fit.fit_leaky_box(data, logpmin=-2, logpmax=-1,
                              pool=pool)
# The output is a DynamicNestedSampler object, the results are here
res = output.results<jupyter_output><empty_output><jupyter_text>Importantly, dynesty returns *weighted* samples of the posterior. If you want to get the equally-weighted samples like a usual MCMC, you need to resample using the weights. This is easy in dynesty:<jupyter_code>res.keys()
# dyplot.resample_equal?
samples = dyplot.resample_equal(res["samples"], np.exp(res["logwt"] - res["logz"][-1]))

## Note: in dynesty/plotting.py they do this:
# try:
#     weights = np.exp(results['logwt'] - results['logz'][-1])
# except:
#     weights = results['weights']
p05, p50, p95 = np.percentile(samples, [5,50,95])
plt.hist(samples, bins='auto', density=True);
plt.axvline(p05,color='r',ls=':')
plt.axvline(p50,color='r',ls='-')
plt.axvline(p95,color='r',ls=':')<jupyter_output><empty_output><jupyter_text>You can then project the fitted parameters into the observed space.<jupyter_code>h,x,p = plt.hist(data, bins='auto', density=True)
plt.xlabel("[Fe/H]")
xplot = np.linspace(-4,0)
plt.plot(xplot, mdfmodels.leaky_box(xplot, p50), 'r-')
plt.plot(xplot, mdfmodels.leaky_box(xplot, p05), 'r:')
plt.plot(xplot, mdfmodels.leaky_box(xplot, p95), 'r:')<jupyter_output><empty_output><jupyter_text>The traceplot shows the 1D posteriors and the time evolution of the values.Note that samples are weighted and progress linearly from low to high -lnX.The coloring shows where the typical set (i.e. most of the volume of the posterior) is.<jupyter_code>fig1 = dyplot.traceplot(res)<jupyter_output><empty_output><jupyter_text>This is a standard corner plot that automatically accounts for the weights.<jupyter_code>fig2 = dyplot.cornerplot(res)<jupyter_output><empty_output><jupyter_text>Fitting pre-enriched boxLet's use the same leaky box data and fit it using the pre enriched model.<jupyter_code># this step takes a minute or two because dynesty is kinda slow. There's ways to optimize that I haven't found yet
with MultiPool(4) as pool:
    output2 = fit.fit_pre_enriched_box(data, logpmin=-2, logpmax=-1,
                                       feh0min=-4, feh0max=-2,
                                       pool=pool)
dyplot.traceplot(output2.results);
dyplot.cornerplot(output2.results);<jupyter_output><empty_output><jupyter_text>You can see in the cornerplot that x2, which is feh0, is essentially an upper limit to be less than -3.3. This makes sense, because the pre-enriched box reduces to the leaky box when feh0 goes to -infinity (or at least less than the observations)  Fitting Extra Gas modelRight now this is very slow. The issue is that there is no analytic form of the normalizing constant, and so you have to numerically integrate it every time.I am going to try to precompute and interpolate to speed this up, but here's the code for now.<jupyter_code># this step took 4 hours, need to optimize this for sure.
with MultiPool(6) as pool:
    output3 = fit.fit_extra_gas(data, logpmin=-2, logpmax=-1,
                                Mmin=1, Mmax=3,
                                pool=pool)

# import pickle
# with open("fitexample_extra_gas_dresults.pkl","wb") as fp:
#     pickle.dump(output3.results,fp)<jupyter_output><empty_output><jupyter_text>M=1 reduces the extra gas model back to the leaky box model, so this is all great!<jupyter_code>dyplot.traceplot(output3.results);
dyplot.cornerplot(output3.results);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09423374613003097, 'fraction_numerical': 0.01586687306501548, 'mean_word_length': 4.528342245989305, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '437793', 'license': 'MIT', 'n_tokens_mistral': 1671, 'n_tokens_neox': 1586, 'n_words': 600, 'path': 'FitExample.ipynb', 'repo_name': 'alexji/mdfmodels', 'size': '5168'}