Fitting of stellar mass function with only two, uncertain data points

smf-fit.py

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()