From 89de9a2d5d55782c6b70eab22616b42b158c9ab9 Mon Sep 17 00:00:00 2001 From: e-koch Date: Wed, 13 Aug 2014 14:02:00 -0700 Subject: [PATCH] More testing using scipy fit funcs. --- plnfit.py | 36 +++++++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 9 deletions(-) diff --git a/plnfit.py b/plnfit.py index 2962b63..5891076 100644 --- a/plnfit.py +++ b/plnfit.py @@ -1,23 +1,41 @@ import numpy as np -import scipy.stats as ss -from scipy.stats import norm -from scipy.special import erf, erfc import matplotlib.pyplot as plt import scipy.optimize as opt +from astropy.io.fits import getdata from pln_distrib import pln - def objfunc(x, p): if p[2] <= 0: return np.inf vals = -np.sum(pln.logpdf(x, *p)) return vals -sampledata = np.exp(np.random.randn(1000) + 2) -params = opt.minimize( - lambda p: objfunc(sampledata, p), [0, 0, 1], method='Powell') +# img = getdata("../testingdata/chamaeleonI-250.fits") +# sampledata = img[np.isfinite(img)] +# sampledata = sampledata + np.abs(sampledata.min()) + 0.1 # [sampledata > 0] + + +# sampledata = np.exp(np.random.randn(1000) + 2) +sampledata = pln.rvs(2, 2, 0.1, size=1000) + 0.1 * np.random.randn(1000) + +guess = [1.8, 2.1, 0.3] +params = opt.minimize( + lambda p: objfunc(sampledata, p), guess, method='Powell') +print("Finished minimization.") p = params['x'] -plt.hist(sampledata, normed=1, log=True) +plt.hist(sampledata, bins=100, normed=1, log=True) trialx = np.linspace(sampledata.min(), sampledata.max(), 1000) -plt.plot(trialx, pln.pdf(trialx, *p)) +plt.plot(trialx, pln.logpdf(trialx, *p), "b-", label="Minimize Powell") + +params2 = opt.basinhopping(lambda p: objfunc(sampledata, p), guess, + minimizer_kwargs={"method": "Powell"}, niter=100) +print("Finished basinhopping.") +p2 = params2['x'] +plt.plot(trialx, pln.logpdf(trialx, *p2), "r-", label="basinhopping") + + +p3 = pln.fit(sampledata, floc=0, fscale=1.0) +plt.plot(trialx, pln.logpdf(trialx, *p3), "g-", label="MaxLik - built-in.") +plt.legend() +plt.show()