In [1]:
import lmfit
In [2]:
def gaussian(x, amp, cen, wid):
return amp * exp(-(x-cen)**2 /wid)
In [3]:
from scipy.optimize import curve_fit
In [16]:
x = linspace(-5., 5., 100)
y = gaussian(x, 1., 0., 1.)
yran = y+randn(y.size)*0.05
plot(x,yran, 'o')
Out[16]:
In [21]:
init_vals = [2,1,2]
best_vals, covar = curve_fit(gaussian, x, yran, p0=init_vals)
print(best_vals, best_vals-[1,0,1])
plot(x, yran, 'o')
plot(x, gaussian(x, best_vals[0], best_vals[1], best_vals[2]))
Out[21]:
In [22]:
gmod = lmfit.Model(gaussian)
In [23]:
print(gmod.param_names)
In [25]:
print(gmod.independent_vars)
In [26]:
params = gmod.make_params()
In [28]:
gmod.eval(x=x, amp=10, cen=6.2, wid=0.75)
Out[28]:
In [1]:
from lmfit.models import GaussianModel, ExponentialModel
data = loadtxt('/Users/Benjamin/PythonLib/lmfit/examples/NIST_Gauss2.dat')
x = data[:, 1]
y = data[:, 0]
In [2]:
exp_mod = ExponentialModel(prefix='exp_')
pars = exp_mod.guess(y, x=x)
gauss1 = GaussianModel(prefix='g1_')
pars.update( gauss1.make_params())
pars['g1_center'].set(105, min=75, max=125)
pars['g1_sigma'].set(15, min=3)
pars['g1_amplitude'].set(2000, min=10)
gauss2 = GaussianModel(prefix='g2_')
pars.update(gauss2.make_params())
pars['g2_center'].set(155, min=125, max=175)
pars['g2_sigma'].set(15, min=3)
pars['g2_amplitude'].set(2000, min=10)
mod = gauss1 + gauss2 + exp_mod
init = mod.eval(pars, x=x)
plt.plot(x, y)
plt.plot(x, init, 'k--')
out = mod.fit(y, pars, x=x)
print(out.fit_report(min_correl=0.5))
plot(x, out.best_fit, 'r-')
Out[2]:
In [37]:
model = lmfit.models.QuadraticModel()
In [38]:
pars = model.make_params()
In [44]:
data = loadtxt('temp.txt')
wave = data[:,0]
flux = data[:,1]
In [224]:
nlines = 2
lines = zeros(nlines, dtype=[('SIGN', 'i'),('ELEMENT','S20'),('WAVE','f4'),('EW','f4'), ('WAVELEFT', 'f4')])
lines[0] = (-1, 'MgII', 2796.35, 2., 2789.)
lines[1] = (-1, 'MgII', 2803.53, 2., 2798.)
model, pars = ebossspec.make_model(lines)
In [208]:
reload(ebossspec)
Out[208]:
In [132]:
In [154]:
init = model.eval(pars, x=log10(wave))
out = model.fit(flux,pars,x=log10(wave), weights=1./power(0.05*ones(wave.size),2))
In [157]:
out.fit_report()
Out[157]:
In [205]:
figsize(10,7)
plot(wave, flux, 'o')
#plot(log10(wave), model.eval(pars,x=log10(wave)))
plot(wave, out.best_fit, lw=2)
xlim(2780, 2810)
plot([2796,2796],[0.5,1.2])
Out[205]:
In [100]:
model
Out[100]:
In [101]:
pars
Out[101]:
In [229]:
print((log10(2801.5)-log10(2803.53))/1E-4*69.)
In [225]:
ew_profile = ebossspec.line_property(log10(wave), flux-1., lines)
In [228]:
x = ravel(ew_profile[1]['WAVE'])
y0 = -ravel(ew_profile[1]['EW'])
y = y0/max(y0)*100.
plot(x, y, '-o')
plot([2798, 2808], [50,50])
plot([2803.53, 2803.53], [0,100])
plot([2801.5, 2801.5], [0,100])
Out[228]:
In [230]:
x = ravel(ew_profile[1]['VEL'])
y0 = -ravel(ew_profile[1]['EW'])
y = y0/max(y0)*100.
plot(x, y, '-o')
plot([-400, 600], [50,50])
plot([-217, -217], [0,100])
plot([0, 0], [0,100])
Out[230]:
In [231]:
import fitsio
In [239]:
data = fitsio.read('/Users/Benjamin/AstroData/DEEP2/deep2.ppxf.specdata.dr4_v1.0.fits.gz')
In [269]:
oii1_wave = data['OII_3727_1_WAVE']
oii2_wave = data['OII_3727_2_WAVE']
oii1_flux = data['OII_3727_1_AMP']
oii2_flux = data['OII_3727_2_AMP']
oii1_EW = data['OII_3727_1_EW']
oii2_EW = data['OII_3727_2_EW']
oii_EW = data['OII_3727_EW']
In [273]:
ii = (where(np.logical_and(oii1_EW>0, oii2_EW>0)))[0]
sort_index = argsort(oii_EW[ii])
print(oii1_EW.size, ii.size, oii_EW.size)
In [276]:
ew_bin = arange(0, 34000, 1000)
ew = zeros(ew_bin.size-1)
ratio = zeros(ew_bin.size-1)
for i in arange(ew.size-1):
ew[i] = median(oii_EW[ii[sort_index[ew_bin[i]:ew_bin[i+1]+1]]])
ratio[i] = median(oii2_EW[ii[sort_index[ew_bin[i]:ew_bin[i+1]+1]]]/oii1_EW[ii[sort_index[ew_bin[i]:ew_bin[i+1]+1]]])
#print(ew, ratio)
print(ew, ratio)
fig = figure(figsize=(8,8))
ax = fig.add_subplot(111)
ax.plot(oii_EW[ii], oii2_EW[ii]/oii1_EW[ii], 'o', ms=1)
#ax.plot(ew, ratio, 'o', ms=20)
ax.set_xlim(-5,40)
ax.set_ylim(0,3)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(2.0)
ax.tick_params(axis='both', which='both', length=7, width=2, labelsize=15)
In [277]:
import quicklook
In [283]:
jj = (where(np.logical_and(oii_EW<50, oii_EW>5)))[0]
quicklook.plot_mean(oii_EW, oii2_EW/oii1_EW, xbin=[0,50], ybin=[0,3], statistic='median')
In [14]:
continuum = lmfit.models.PowerLawModel()
powerlaw = lambda x, amp, index: amp * (x**index)
In [115]:
num_points = 40
amp = 10.0
index = -2.0
#xdata = logspace(log10(1.), log10(10.1), num_points)
xdata = linspace(1.1, 10.1, num_points)
ydata0 = powerlaw(xdata, 10.0, -0.5)
yerr = 0.5*ydata0
ydata = ydata0+randn(num_points)*yerr
ydata
Out[115]:
In [113]:
pars = continuum.guess(ydata, x=xdata)
pars['exponent'].set(min=-3, max=-0.1)
pars['amplitude'].set(min=9, max=11.)
out = continuum.fit(ydata, pars, x=xdata, weights=1./yerr)
#out = continuum.fit(ydata, pars, x=xdata)
print(out.fit_report())
In [114]:
figsize(15,10)
plot(xdata, out.init_fit, 'k')
plot(xdata, out.best_fit, 'g')
errorbar(xdata, ydata, yerr=yerr, fmt='k^')
xlim(1E0, 1E1)
Out[114]:
In [101]:
sum(pow(ydata-ydata0, 2)/pow(yerr,2))/40.
Out[101]:
In [ ]:
def residual(params, x, data, eps_data):
amp = params['amp'].value
index = params['index'].value
model = powerlaw(x, amp, index)
return (data-model)/eps_data
params = Parameters()
params.add('amp', value=10)
params.add('decay', value=0.007)
params.add('phase', value=0.2)
params.add('frequency', value=3.0)
out = minimize(residual, params, args=(x, data, eps_data))