本文整理汇总了Python中scipy.odr.ODR.set_job方法的典型用法代码示例。如果您正苦于以下问题:Python ODR.set_job方法的具体用法?Python ODR.set_job怎么用?Python ODR.set_job使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类scipy.odr.ODR
的用法示例。
在下文中一共展示了ODR.set_job方法的9个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: test_explicit
# 需要导入模块: from scipy.odr import ODR [as 别名]
# 或者: from scipy.odr.ODR import set_job [as 别名]
def test_explicit(self):
explicit_mod = Model(
self.explicit_fcn,
fjacb=self.explicit_fjb,
fjacd=self.explicit_fjd,
meta=dict(name='Sample Explicit Model',
ref='ODRPACK UG, pg. 39'),
)
explicit_dat = Data([0.,0.,5.,7.,7.5,10.,16.,26.,30.,34.,34.5,100.],
[1265.,1263.6,1258.,1254.,1253.,1249.8,1237.,1218.,1220.6,
1213.8,1215.5,1212.])
explicit_odr = ODR(explicit_dat, explicit_mod, beta0=[1500.0, -50.0, -0.1],
ifixx=[0,0,1,1,1,1,1,1,1,1,1,0])
explicit_odr.set_job(deriv=2)
explicit_odr.set_iprint(init=0, iter=0, final=0)
out = explicit_odr.run()
assert_array_almost_equal(
out.beta,
np.array([1.2646548050648876e+03, -5.4018409956678255e+01,
-8.7849712165253724e-02]),
)
assert_array_almost_equal(
out.sd_beta,
np.array([1.0349270280543437, 1.583997785262061, 0.0063321988657267]),
)
assert_array_almost_equal(
out.cov_beta,
np.array([[4.4949592379003039e-01, -3.7421976890364739e-01,
-8.0978217468468912e-04],
[-3.7421976890364739e-01, 1.0529686462751804e+00,
-1.9453521827942002e-03],
[-8.0978217468468912e-04, -1.9453521827942002e-03,
1.6827336938454476e-05]]),
)
示例2: scipyODR
# 需要导入模块: from scipy.odr import ODR [as 别名]
# 或者: from scipy.odr.ODR import set_job [as 别名]
def scipyODR(recipe, *args, **kwargs):
from scipy.odr import Data, Model, ODR, RealData, odr_stop
# FIXME
# temporarily change _weights to _weights**2 to fit the ODR fits
a = [w ** 2 for w in recipe._weights]
recipe._weights = a
model = Model(recipe.evaluateODR,
# implicit=1,
meta=dict(name='ODR fit'),
)
x = [recipe._contributions.values()[0].profile.x]
y = [recipe._contributions.values()[0].profile.y]
dy = [recipe._contributions.values()[0].profile.dy]
cont = recipe._contributions.values()
for i in range(1, len(cont)):
xplus = x[-1][-1] - x[-1][-2] + x[-1][-1]
x.append(cont[i].profile.x + x[-1][-1] + xplus)
y.append(cont[i].profile.y)
dy.append(cont[i].profile.dy)
x.append(np.arange(len(recipe._restraintlist))
* (x[-1][-1] - x[-1][-2]) + x[-1][-1])
y.append(np.zeros_like(recipe._restraintlist))
dy = np.concatenate(dy)
dy = np.concatenate(
[dy, np.ones_like(recipe._restraintlist) + np.average(dy)])
data = RealData(x=np.concatenate(x), y=np.concatenate(y), sy=dy)
odr_kwargs = {}
if kwargs.has_key('maxiter'):
odr_kwargs['maxit'] = kwargs['maxiter']
odr = ODR(data, model, beta0=recipe.getValues(), **odr_kwargs)
odr.set_job(deriv=1)
out = odr.run()
# out.pprint()
# FIXME
# revert back
a = [np.sqrt(w) for w in recipe._weights]
recipe._weights = a
return {'x': out.beta,
'esd': out.sd_beta,
'cov': out.cov_beta,
'raw': out,
}
示例3: test_multi
# 需要导入模块: from scipy.odr import ODR [as 别名]
# 或者: from scipy.odr.ODR import set_job [as 别名]
def test_multi(self):
multi_mod = Model(
self.multi_fcn,
meta=dict(name='Sample Multi-Response Model',
ref='ODRPACK UG, pg. 56'),
)
multi_x = np.array([30.0, 50.0, 70.0, 100.0, 150.0, 200.0, 300.0, 500.0,
700.0, 1000.0, 1500.0, 2000.0, 3000.0, 5000.0, 7000.0, 10000.0,
15000.0, 20000.0, 30000.0, 50000.0, 70000.0, 100000.0, 150000.0])
multi_y = np.array([
[4.22, 4.167, 4.132, 4.038, 4.019, 3.956, 3.884, 3.784, 3.713,
3.633, 3.54, 3.433, 3.358, 3.258, 3.193, 3.128, 3.059, 2.984,
2.934, 2.876, 2.838, 2.798, 2.759],
[0.136, 0.167, 0.188, 0.212, 0.236, 0.257, 0.276, 0.297, 0.309,
0.311, 0.314, 0.311, 0.305, 0.289, 0.277, 0.255, 0.24, 0.218,
0.202, 0.182, 0.168, 0.153, 0.139],
])
n = len(multi_x)
multi_we = np.zeros((2, 2, n), dtype=float)
multi_ifixx = np.ones(n, dtype=int)
multi_delta = np.zeros(n, dtype=float)
multi_we[0,0,:] = 559.6
multi_we[1,0,:] = multi_we[0,1,:] = -1634.0
multi_we[1,1,:] = 8397.0
for i in range(n):
if multi_x[i] < 100.0:
multi_ifixx[i] = 0
elif multi_x[i] <= 150.0:
pass # defaults are fine
elif multi_x[i] <= 1000.0:
multi_delta[i] = 25.0
elif multi_x[i] <= 10000.0:
multi_delta[i] = 560.0
elif multi_x[i] <= 100000.0:
multi_delta[i] = 9500.0
else:
multi_delta[i] = 144000.0
if multi_x[i] == 100.0 or multi_x[i] == 150.0:
multi_we[:,:,i] = 0.0
multi_dat = Data(multi_x, multi_y, wd=1e-4/np.power(multi_x, 2),
we=multi_we)
multi_odr = ODR(multi_dat, multi_mod, beta0=[4.,2.,7.,.4,.5],
delta0=multi_delta, ifixx=multi_ifixx)
multi_odr.set_job(deriv=1, del_init=1)
out = multi_odr.run()
assert_array_almost_equal(
out.beta,
np.array([4.3799880305938963, 2.4333057577497703, 8.0028845899503978,
0.5101147161764654, 0.5173902330489161]),
)
assert_array_almost_equal(
out.sd_beta,
np.array([0.0130625231081944, 0.0130499785273277, 0.1167085962217757,
0.0132642749596149, 0.0288529201353984]),
)
assert_array_almost_equal(
out.cov_beta,
np.array([[0.0064918418231375, 0.0036159705923791, 0.0438637051470406,
-0.0058700836512467, 0.011281212888768],
[0.0036159705923791, 0.0064793789429006, 0.0517610978353126,
-0.0051181304940204, 0.0130726943624117],
[0.0438637051470406, 0.0517610978353126, 0.5182263323095322,
-0.0563083340093696, 0.1269490939468611],
[-0.0058700836512467, -0.0051181304940204, -0.0563083340093696,
0.0066939246261263, -0.0140184391377962],
[0.011281212888768, 0.0130726943624117, 0.1269490939468611,
-0.0140184391377962, 0.0316733013820852]]),
)
示例4: linear_beta
# 需要导入模块: from scipy.odr import ODR [as 别名]
# 或者: from scipy.odr.ODR import set_job [as 别名]
return y
def linear_beta(x, beta_0, beta_1):
return beta_0 + beta_1 * x
intensity = np.array(As).reshape(9, 5).T[[0, 3, 4, 1, 2]].T # move peak 2 & 3 to the correct position
A_lambda = intensity[:-1].T # entries of each peak put together
peaks = ['558.3', '563.5', '564.7', '570.8', '576.7']
intersects = []
x = np.array(concentration[:-1])
sx = 4
for i, peak, A in enum(peaks, A_lambda):
y = un.nominal_values(A)
sy = un.std_devs(A)
data = RealData(x, y, sx=sx, sy=sy)
model = Model(func)
odr = ODR(data, model, [6, 0])
odr.set_job(fit_type=0)
output = odr.run()
beta = uc.correlated_values(output.beta, output.cov_beta)
fit = linear_beta(x_fit, *beta)
#Intersects
y0 = intensity[-1][i] # peak intensities of the original sample
x0 = (y0 - beta[0]) / beta[1]
intersects.append(x0)
lamb_all = np.array(x0s).reshape(9, 5)
d_lin = np.load(npy_dir + 'ccd_calibration.npy')
lamb_all = linear(lamb_all, *d_lin) # Integrate error on wavelength of CCD
lamb_mean = np.mean(lamb_all, 0)
lamb_laser = np.load(npy_dir + 'ccd_lamb_laser.npy')[0]
dnu = (1 / lamb_mean - 1 / lamb_laser) * 10**7
lamb_mean = np.sort(lamb_mean)
lits = np.array([888, 1028, 1091, 1274, 1464])
示例5: zeroth_fit_2var
# 需要导入模块: from scipy.odr import ODR [as 别名]
# 或者: from scipy.odr.ODR import set_job [as 别名]
def zeroth_fit_2var(Nml,Eth,nu_des_ini,t_start,t_stop):
deca=0
sizeOfFont = 16
fontProperties = {'family':'sans-serif',
'weight' : 'normal', 'size' : sizeOfFont}
#rc('text', usetex=True)
rc('font',**fontProperties)
ticks_font = font_manager.FontProperties(family='cm', style='normal',
size=sizeOfFont, weight='normal', stretch='normal')
filenm='12CO_13CO_1_1_fit.txt'
#if molec =='n2':
# m_mol=30.0e-3
# xtex=0.042
#if molec =='co':
m_mol=28.0e-3
xtex=0.0375
#if molec =='co2':
# m_mol=44.0e-3
#if molec =='h2o':
# m_mol=18.0e-3
errort=0.1
T_tot, N_tot,N_13co = np.loadtxt(filenm,unpack=True,skiprows=1)
N_tot=N_13co
T_ini=T_tot
T_tot=T_tot
T_tot_err=T_tot*0.0+errort
N_tot=np.clip(N_tot,1e-12,1e30)
N_tot_err=1e12+N_tot*0
h_rate=1/60.0 #Kmin-1
T=T_tot[np.where((T_tot>t_start) & (T_tot<t_stop))]
N=N_tot[np.where((T_tot>t_start) & (T_tot<t_stop))]
N_err=1e10+N*0
T_err=T*0.0+errort
x=1.0/T
x_err=T_err/(T**2)
#print x_err
y=np.log(N/(1e15))
y_err=N_err/N
#y=np.log(N)
p=[0.,0.]
#nu_des_ini=7e11#np.sqrt(2.0*1e19*8.31*Eth/(np.pi*np.pi*m_mol))
#slope, intercept, r_value, p_value, std_err=scipy.stats.linregress(x, y)
def func(p, t):
return p[1] +t*p[0]
#def func (p,t):
# return np.log(1e12/h_rate)+t*p[0]
#def func(p, t):
# return np.log(np.sqrt(2.0*1e19*8.31*-1*p[0]/(np.pi*np.pi*m_mol))/h_rate) +t*p[0]
#p[1]=np.log(1e12/h_rate)
#p,pcov=curve_fit(func, x, y)
#err = np.sqrt(np.diag(pcov))
data = RealData(x, y, sx=x_err,sy=y_err)#, sy=np.log(0.001))
model = Model(func)
#odr = ODR(data, model, beta0=[-800,1e12])
odr = ODR(data, model, beta0=[-750,7e11])
odr.set_job(fit_type=2)
output = odr.run()
p=output.beta
err=output.sd_beta
E_des=-p[0]
nu_des=h_rate*np.exp(p[1])
nu_theo=h_rate*np.exp(np.log(np.sqrt(2.0*1e19*8.31*-1*p[0]/(np.pi*np.pi*m_mol))/h_rate))
#p[1]=np.log(np.sqrt(2.0*1e19*8.31*-1*p[0]/(np.pi*np.pi*m_mol))/h_rate)
#err[1]=0.0
#nu_theo=np.sqrt(2.0*1e19*8.31*Eth/(np.pi*np.pi*m_mol))
print 'E_des fit', p[0]*-1, '+-',err[0]
print 'E_des litt', Eth
print 'nu fit e12', 1e-12*h_rate*np.exp(p[1]), '+-',err[1]*1e-12*h_rate*np.exp(p[1])
print 'nu litt e12', 1e-12*nu_des_ini
print 'nu theo e12', 1e-12*nu_theo
#nu_des=2e11
#print "nu e12n", nu_des*1e-12
# def funct(x,t):
# if x[0]>0.0:
# f_ml = -(nu_des/h_rate)*np.exp(-E_des/t)
# f_g= f_ml-x[1]
# if x[0]<0.0:
# f_ml=0
#.........这里部分代码省略.........
示例6: main
# 需要导入模块: from scipy.odr import ODR [as 别名]
# 或者: from scipy.odr.ODR import set_job [as 别名]
def main():
try:
if (len(sys.argv) > 1 and sys.argv[1] == '1'):
data_filename = 'number.txt'
start_time = time.time()
number = np.zeros((41,3))
index = -1
for det in range(1,8):
if (det!=6):
num_n = 6
else:
num_n = 5
for n in range(0,num_n):
index = index + 1
temp = np.zeros(3)
for angle in range(1,4):
filename = 'detector_' + str(det) + '_no_' + str(n) \
+ '_angle_' + str(angle) + '.jpg'
refname = 'detector_' + str(det) + '_no_' + str(n) \
+ '_background.jpg'
elapse_time = (time.time()-start_time)
if(elapse_time >= 1):
remain_time = elapse_time/(index*3+angle-1)*41*3-elapse_time
print 'Processing .. ' + filename \
+ time.strftime(" %H:%M:%S", time.gmtime(elapse_time)) \
+ ' has past. ' + 'Remaining time: ' \
+ time.strftime(" %H:%M:%S", time.gmtime(remain_time))
else:
print 'Processing .. ' + filename \
+ time.strftime(" %H:%M:%S", time.gmtime(elapse_time)) \
+ ' has past'
image = rgb2gray(io.imread(filename))
ref = rgb2gray(io.imread(refname))
temp[angle-1] = ellipse.count_bubble(image,ref)
#temp[angle-1] = ellipse.count_bubble(image)
number[index,1] = np.mean(temp)
number[index,2] = np.std(temp)
manual_count = np.array([1,27,40,79,122,160,1,18,28,42,121,223,0,11,24,46,\
142,173,3,19,23,76,191,197,0,15,24,45,91,152,0,\
16,27,34,88,0,9,12,69,104,123])
number[:,0] = manual_count.T
number.tofile(data_filename,sep=" ")
elif(len(sys.argv) > 1):
data_filename = sys.argv[1]
data = np.loadtxt(data_filename,skiprows=0)
number = np.reshape(data,(41,3))
else:
data_filename = 'number.txt'
data = np.loadtxt(data_filename,skiprows=0)
number = np.reshape(data,(41,3))
#emperical error
xerr = np.sqrt(number[:,0])/3
yerr = number[:,2]
data = Data(number[:,0].T, number[:,1].T, we = 1/(np.power(xerr.T,2)+np.spacing(1)), wd = 1/(np.power(yerr.T,2)+np.spacing(1)))
model = Model(ord_function)
odr = ODR(data, model, beta0=[0, 0])
odr.set_job(fit_type=2)
output = odr.run()
popt = output.beta
perr = output.sd_beta
output.pprint()
# popt, pcov = curve_fit(linear_fit_function, number[:,0], number[:,1], [0, 0], number[:, 2])
# perr = np.sqrt(np.diag(pcov))
fitting_error = np.mean(np.sqrt(np.power(popt[0]*number[:,0]+popt[1] - number[:,1],2)))
labels = np.array([[1,1,1,1,1,1,\
2,2,2,2,2,2,\
3,3,3,3,3,3,\
4,4,4,4,4,4,\
5,5,5,5,5,5,\
6,6,6,6,6,\
7,7,7,7,7,7],
[0,1,2,3,4,5,\
0,1,2,3,4,5,\
0,1,2,3,4,5,\
0,1,2,3,4,5,\
0,1,2,3,4,5,\
0,1,2,3,4,
0,1,2,3,4,5]])
fig, ax = plt.subplots(ncols = 1)
ax.errorbar(number[:,0], number[:,1], xerr = xerr, yerr = yerr, fmt='o')
ax.plot(number[:,0], popt[0]*number[:,0]+popt[1], '-r')
bbox_props = dict(boxstyle="square,pad=0.3", fc="white", ec="black", lw=2)
annotation_text = "function: y = kx + b \n" \
"k = %.2f" % popt[0] + " +/- %.2f" % perr[0] + '\n' \
#.........这里部分代码省略.........
示例7: show
# 需要导入模块: from scipy.odr import ODR [as 别名]
# 或者: from scipy.odr.ODR import set_job [as 别名]
def show():
data_filename = 'number.txt'
data = np.loadtxt(data_filename,skiprows=0)
number = np.reshape(data,(41,3))
#emperical error
xerr = np.sqrt(number[:,0])/3
yerr = number[:,2]
data = Data(number[:,0].T, number[:,1].T, we = 1/(np.power(xerr.T,2)+np.spacing(1)), wd = 1/(np.power(yerr.T,2)+np.spacing(1)))
model = Model(ord_function)
odr = ODR(data, model, beta0=[0, 0])
odr.set_job(fit_type=2)
output = odr.run()
popt = output.beta
perr = output.sd_beta
output.pprint()
fitting_error = np.mean(np.sqrt(np.power(popt[0]*number[:,0]+popt[1] - number[:,1],2)))
labels = np.array([[1,1,1,1,1,1,\
2,2,2,2,2,2,\
3,3,3,3,3,3,\
4,4,4,4,4,4,\
5,5,5,5,5,5,\
6,6,6,6,6,\
7,7,7,7,7,7],
[0,1,2,3,4,5,\
0,1,2,3,4,5,\
0,1,2,3,4,5,\
0,1,2,3,4,5,\
0,1,2,3,4,5,\
0,1,2,3,4,
0,1,2,3,4,5]])
fig, ax = plt.subplots(ncols = 1)
ax.errorbar(number[:,0], number[:,1], xerr = xerr, yerr = yerr, fmt='o')
ax.plot(number[:,0], popt[0]*number[:,0]+popt[1], '-r')
bbox_props = dict(boxstyle="square,pad=0.3", fc="white", ec="black", lw=2)
annotation_text = "function: y = kx + b \n" \
"k = %.2f" % popt[0] + " +/- %.2f" % perr[0] + '\n' \
"b = %.2f" % popt[1] + " +/- %.2f" % perr[1] + '\n' \
"Error: %.2f" % fitting_error
ax.text(10, np.amax(number[:,1])+10, annotation_text, ha="left", va="top", rotation=0,
size=15, bbox=bbox_props)
for i in range(0, len(number[:,0])): # <--
ax.annotate('(%s, %s)' % (labels[0,i], labels[1,i]),\
(number[i,0]+1, number[i,1]-1)) #
ax.set_title('Algorithom Performance')
ax.set_xlabel('Bubble Number Counted Manually')
ax.set_ylabel('Bubbble Number Counted by Algorithom')
plt.grid()
plt.xlim((np.amin(number[:,0])-5,np.amax(number[:,0])+5))
plt.ylim((0,np.amax(number[:,1])+20))
plt.show()
示例8: main
# 需要导入模块: from scipy.odr import ODR [as 别名]
# 或者: from scipy.odr.ODR import set_job [as 别名]
def main():
try:
legend = ['curvature', 'neural-network', 'hough transform', 'SVM-Preliminary']
data_filename = ['number.txt', 'skflow.txt', 'curvature.txt', 'SVM.txt']
colors = ['blue', 'red', 'black', 'green']
annotation_text = "function: y = kx + b";
fig, ax = plt.subplots(ncols = 1)
for i in range(len(data_filename)):
temp = np.loadtxt(data_filename[i],skiprows=0)
data = np.reshape(temp,(41,len(temp)/41))
#print data
data[:,1:3] = np.true_divide(data[:,1:3], np.amax(data[:,1]))
data[:,1] = data[:,1]+i
#print data
#emperical error
xerr = np.sqrt(data[:,0])/3
yerr = data[:,2]
data_fit = Data(data[:,0].T, data[:,1].T, \
we = 1/(np.power(xerr.T,2)+np.spacing(1)),\
wd = 1/(np.power(yerr.T,2)+np.spacing(1)))
model = Model(ord_function)
odr = ODR(data_fit, model, beta0=[0, 0])
odr.set_job(fit_type=2)
output = odr.run()
popt = output.beta
perr= output.sd_beta
popt, pcov = curve_fit(linear_fit_function,
data[:,0], data[:,1], [1, 0], data[:, 2])
perr = np.sqrt(np.diag(pcov))
A = popt[0]/np.sqrt(popt[0]*popt[0]+1)
B = -1/np.sqrt(popt[0]*popt[0]+1)
C = popt[1]/np.sqrt(popt[0]*popt[0]+1)
fitting_error= np.mean(np.abs(A*data[:,0]+B*data[:,1]+C))
ax.errorbar(data[:,0], data[:,1], xerr = xerr,\
yerr = yerr, fmt='o',color=colors[i])
### not using error bar in fitting #########
popt[0],popt[1] = np.polyfit(data[:,0], data[:,1], 1)
###
ax.plot(data[:,0], popt[0]*data[:,0]+popt[1], colors[i], linewidth = 2)
annotation_text = annotation_text + "\n" +\
legend[i] + "(" + \
colors[i] + ") \n" + \
"k = %.2f b = %.2f Error = %.2f"%( \
popt[0], popt[1], fitting_error)
bbox_props = dict(boxstyle="square,pad=0.3", fc="white", ec="black", lw=2)
#ax.text(0, 4.15, annotation_text, ha="left", va="top", \
# rotation=0, size=14, bbox=bbox_props)
ax.set_title('Algorithom Performance')
ax.set_xlabel('Bubble Number Counted Manually')
ax.set_ylabel('Normalized Bubbble Number Counted by Algorithom')
plt.grid()
plt.xlim((np.amin(data[:,0])-5,np.amax(data[:,0])+5))
plt.ylim((0,np.amax(data[:,1])+0.2))
plt.show()
except KeyboardInterrupt:
print "Shutdown requested... exiting"
except Exception:
traceback.print_exc(file=sys.stdout)
sys.exit(0)
示例9: fitDoubleGaussians
# 需要导入模块: from scipy.odr import ODR [as 别名]
# 或者: from scipy.odr.ODR import set_job [as 别名]
def fitDoubleGaussians(charges, charge_err, data, data_err):
# some simple peakfinding to find the first and second peaks
peak1_value = np.max(data)
peak1_bin = [x for x in range(len(data)) if data[x]==peak1_value][0]
peak1_charge = charges[peak1_bin]
# Preliminary fits to get better estimates of parameters...
first_peak_data = RealData(charges[peak1_bin-30:peak1_bin+30],
data[peak1_bin-30:peak1_bin+30],
charge_err[peak1_bin-30:peak1_bin+30],
data_err[peak1_bin-30:peak1_bin+30],)
first_peak_model = Model(gauss_fit)
first_peak_odr = ODR(first_peak_data, first_peak_model,
[peak1_value, peak1_charge, 0.1*peak1_charge])
first_peak_odr.set_job(fit_type=2)
first_peak_output = first_peak_odr.run()
first_peak_params = first_peak_output.beta
#second_peak_params, covariance = curve_fit(gauss_fit2,
# data[peak1_bin-30:peak1_bin+30],
# data[peak1_bin-30:peak1_bin+30],
# p0 = [peak1_value, peak1_charge, 0.1*peak1_charge])
# subtract the largest peak so we can search for the other one
updated_data = data-gauss_fit(first_peak_params, charges)
#updated_data = data[:int(len(data)*2.0/3.0)]
peak2_value = np.max(updated_data)
peak2_bin = [x for x in range(len(updated_data)) if updated_data[x]==peak2_value][0]
peak2_charge = charges[peak2_bin]
#first_peak_params, covariance = curve_fit(gauss_fit2,
# updated_data[peak2_bin-30:peak2_bin+30],
# updated_data[peak2_bin-30:peak2_bin+30],
# p0 = [peak2_value, peak2_charge, 0.1*peak2_charge])
# and the second peak...
second_peak_data = RealData(charges[peak2_bin-30:peak2_bin+30],
data[peak2_bin-30:peak2_bin+30],
charge_err[peak2_bin-30:peak2_bin+30],
data_err[peak2_bin-30:peak2_bin+30],)
second_peak_model = Model(gauss_fit)
second_peak_odr = ODR(second_peak_data, second_peak_model,
[peak2_value, peak2_charge, first_peak_params[2]])
second_peak_odr.set_job(fit_type=2)
second_peak_output = second_peak_odr.run()
second_peak_params = second_peak_output.beta
# Now fit both gaussians simultaneously
double_peak_data = RealData(charges, data,
charge_err, data_err)
double_peak_model = Model(double_gauss_fit)
double_peak_odr = ODR(double_peak_data, double_peak_model,
[first_peak_params[0], first_peak_params[1], first_peak_params[2],
second_peak_params[0], second_peak_params[1], second_peak_params[2]])
double_peak_odr.set_job(fit_type=0)
double_peak_output = double_peak_odr.run()
double_peak_params = double_peak_output.beta
double_peak_sigmas = double_peak_output.sd_beta
#double_peak_params = [first_peak_params[0], first_peak_params[1], first_peak_params[2],
# second_peak_params[0], second_peak_params[1], second_peak_params[2]]
return double_peak_params, double_peak_sigmas