Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fitting of stellar mass function with only two, uncertain data points
import matplotlib.pyplot as plt
import numpy
from numpy import log, exp
import scipy.misc
import scipy.optimize
# see Buchner et al (2015), Appendix A, for details on the method
# http://adsabs.harvard.edu/abs/2015ApJ...802...89B
Nsamples = 1000
data = [numpy.random.normal(9, 0.3, size=Nsamples),
numpy.random.normal(10, 0.3, size=Nsamples)]
data = numpy.array(data)
def smf(logM, lognorm, logM0, alpha):
return lognorm + (logM-logM0) * (alpha + 1) - exp(logM-logM0)
Ngrid = 1000
logMgrid = None
def loglike((lognorm, logM0, alpha)):
prob = smf(data, lognorm, logM0, alpha)
#print scipy.misc.logsumexp(prob, axis=1).shape
dataterm = (scipy.misc.logsumexp(prob, axis=1) - log(Nsamples)).mean()
detterm = exp(smf(logMgrid, lognorm, logM0, alpha)).mean()
like = dataterm - detterm
return like
def chi2(params):
like = loglike(params)
#print '%.3f' % like, params
return -2 * like
for lo in 6, 7, 8:
for hi in 12, 13, 14:
logMgrid = numpy.linspace(lo, hi, Ngrid)
params = scipy.optimize.fmin(chi2, [-4, 12, 0], disp=False)
print 'lo=%d hi=%d' % (lo, hi), 'lognorm=%.2e logM0=%.1f alpha=%.2f ' % tuple(params)
plt.plot(logMgrid, smf(logMgrid, *params), label='lo=%d hi=%d' % (lo, hi))
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.