Python Minimizer.prepare_fit方法代码示例

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
def NIST_Test(DataSet, start='start2', plot=True):

    NISTdata = ReadNistData(DataSet)
    resid, npar, dimx = Models[DataSet]
    y = NISTdata['y']
    x = NISTdata['x']

    params = Parameters()
    for i in range(npar):
        pname = 'b%i' % (i+1)
        cval  = NISTdata['cert_values'][i]
        cerr  = NISTdata['cert_stderr'][i]
        pval1 = NISTdata[start][i]
        params.add(pname, value=pval1)


    myfit = Minimizer(resid, params, fcn_args=(x,), fcn_kws={'y':y},
                      scale_covar=True)

    myfit.prepare_fit()
    myfit.leastsq()

    digs = Compare_NIST_Results(DataSet, myfit, params, NISTdata)

    if plot and HASPYLAB:
        fit = -resid(params, x, )
        pylab.plot(x, y, 'r+')
        pylab.plot(x, fit, 'ko--')
        pylab.show()

    return digs > 2

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
def test_constraints1():
    def residual(pars, x, sigma=None, data=None):
        yg = gaussian(x, pars['amp_g'], pars['cen_g'], pars['wid_g'])
        yl = lorentzian(x, pars['amp_l'], pars['cen_l'], pars['wid_l'])

        model =  yg +  yl + pars['line_off'] + x * pars['line_slope']
        if data is None:
            return model
        if sigma is None:
            return (model - data)
        return (model - data)/sigma


    n = 601
    xmin = 0.
    xmax = 20.0
    x = linspace(xmin, xmax, n)

    data = (gaussian(x, 21, 8.1, 1.2) +
            lorentzian(x, 10, 9.6, 2.4) +
            random.normal(scale=0.23,  size=n) +
            x*0.5)


    pfit = Parameters()
    pfit.add(name='amp_g',  value=10)
    pfit.add(name='cen_g',  value=9)
    pfit.add(name='wid_g',  value=1)

    pfit.add(name='amp_tot',  value=20)
    pfit.add(name='amp_l',  expr='amp_tot - amp_g')
    pfit.add(name='cen_l',  expr='1.5+cen_g')
    pfit.add(name='wid_l',  expr='2*wid_g')

    pfit.add(name='line_slope', value=0.0)
    pfit.add(name='line_off', value=0.0)

    sigma = 0.021  # estimate of data error (for all data points)

    myfit = Minimizer(residual, pfit,
                      fcn_args=(x,), fcn_kws={'sigma':sigma, 'data':data},
                      scale_covar=True)

    myfit.prepare_fit()
    init = residual(myfit.params, x)

    result = myfit.leastsq()

    print(' Nfev = ', result.nfev)
    print( result.chisqr, result.redchi, result.nfree)

    report_fit(result.params)
    pfit= result.params
    fit = residual(result.params, x)
    assert(pfit['cen_l'].value == 1.5 + pfit['cen_g'].value)
    assert(pfit['amp_l'].value == pfit['amp_tot'].value - pfit['amp_g'].value)
    assert(pfit['wid_l'].value == 2 * pfit['wid_g'].value)

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
def test_peakfit():
    from lmfit.utilfuncs import gaussian
    def residual(pars, x, data=None):
        g1 = gaussian(x, pars['a1'].value, pars['c1'].value, pars['w1'].value)
        g2 = gaussian(x, pars['a2'].value, pars['c2'].value, pars['w2'].value)
        model = g1 + g2
        if data is None:
            return model
        return (model - data)

    n    = 601
    xmin = 0.
    xmax = 15.0
    noise = np.random.normal(scale=.65, size=n)
    x = np.linspace(xmin, xmax, n)

    org_params = Parameters()
    org_params.add_many(('a1', 12.0, True, None, None, None),
                        ('c1',  5.3, True, None, None, None),
                        ('w1',  1.0, True, None, None, None),
                        ('a2',  9.1, True, None, None, None),
                        ('c2',  8.1, True, None, None, None),
                        ('w2',  2.5, True, None, None, None))

    data  = residual(org_params, x) + noise


    fit_params = Parameters()
    fit_params.add_many(('a1',  8.0, True, None, 14., None),
                        ('c1',  5.0, True, None, None, None),
                        ('w1',  0.7, True, None, None, None),
                        ('a2',  3.1, True, None, None, None),
                        ('c2',  8.8, True, None, None, None))

    fit_params.add('w2', expr='2.5*w1')

    myfit = Minimizer(residual, fit_params,
                      fcn_args=(x,), fcn_kws={'data':data})

    myfit.prepare_fit()

    init = residual(fit_params, x)


    myfit.leastsq()

    print(' N fev = ', myfit.nfev)
    print(myfit.chisqr, myfit.redchi, myfit.nfree)

    report_fit(fit_params)

    fit = residual(fit_params, x)
    check_paras(fit_params, org_params)

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
	def fitevent(self, edat, initguess):
		try:
			dt = 1000./self.Fs 	# time-step in ms.

			# control numpy error reporting
			np.seterr(invalid='ignore', over='ignore', under='ignore')

			ts = np.array([ t*dt for t in range(0,len(edat)) ], dtype='float64')

			self.nStates=len(initguess)
			initRCConst=dt*5.

			# setup fit params
			params=Parameters()

			for i in range(0, len(initguess)):
				params.add('a'+str(i), value=initguess[i][0]) 
				params.add('mu'+str(i), value=initguess[i][1]) 
				if self.LinkRCConst:				
					if i==0:
						params.add('tau'+str(i), value=initRCConst)
					else:
						params.add('tau'+str(i), value=initRCConst, expr='tau0')
				else:
					params.add('tau'+str(i), value=initRCConst)

			params.add('b', value=self.baseMean )
			

			igdict=params.valuesdict()

			optfit=Minimizer(self._objfunc, params, fcn_args=(ts,edat,))
			optfit.prepare_fit()
			result=optfit.leastsq(xtol=self.FitTol,ftol=self.FitTol,maxfev=self.FitIters)

			if result.success:
				tt=[init[0] for init, final in zip(igdict.items(), (result.params.valuesdict()).items()) if init==final]
				if len(tt) > 0:
					self.flagEvent('wInitGuessUnchanged')

				self._recordevent(result)
			else:
				#print optfit.message, optfit.lmdif_message
				self.rejectEvent('eFitConvergence')
		except KeyboardInterrupt:
			self.rejectEvent('eFitUserStop')
			raise
		except InvalidEvent:
			self.rejectEvent('eInvalidEvent')
		except:
	 		self.rejectEvent('eFitFailure')

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
	def __FitEvent(self):
		try:
			dt = 1000./self.Fs 	# time-step in ms.
			# edat=np.asarray( np.abs(self.eventData),  dtype='float64' )
			edat=self.dataPolarity*np.asarray( self.eventData,  dtype='float64' )

			# control numpy error reporting
			np.seterr(invalid='ignore', over='ignore', under='ignore')

			ts = np.array([ t*dt for t in range(0,len(edat)) ], dtype='float64')

			# estimate initial guess for events
			initguess=self._characterizeevent(edat, np.abs(util.avg(edat[:10])), self.baseSD, self.InitThreshold, 6.)
			self.nStates=len(initguess)-1

			# setup fit params
			params=Parameters()

			for i in range(1, len(initguess)):
				params.add('a'+str(i-1), value=initguess[i][0]-initguess[i-1][0]) 
				params.add('mu'+str(i-1), value=initguess[i][1]*dt) 
				params.add('tau'+str(i-1), value=dt*7.5)

			params.add('b', value=initguess[0][0])
			

			optfit=Minimizer(self.__objfunc, params, fcn_args=(ts,edat,))
			optfit.prepare_fit()

	
			optfit.leastsq(xtol=self.FitTol,ftol=self.FitTol,maxfev=self.FitIters)

			if optfit.success:
				self.__recordevent(optfit)
			else:
				#print optfit.message, optfit.lmdif_message
				self.rejectEvent('eFitConvergence')
		except KeyboardInterrupt:
			self.rejectEvent('eFitUserStop')
			raise
		except InvalidEvent:
			self.rejectEvent('eInvalidEvent')
		except:
	 		self.rejectEvent('eFitFailure')
	 		raise

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
def test_peakfit():
    def residual(pars, x, data=None):
        g1 = gaussian(x, pars['a1'], pars['c1'], pars['w1'])
        g2 = gaussian(x, pars['a2'], pars['c2'], pars['w2'])
        model = g1 + g2
        if data is None:
            return model
        return (model - data)

    n = 601
    xmin = 0.
    xmax = 15.0
    noise = np.random.normal(scale=.65, size=n)
    x = np.linspace(xmin, xmax, n)

    org_params = Parameters()
    org_params.add_many(('a1', 12.0, True, None, None, None),
                        ('c1', 5.3, True, None, None, None),
                        ('w1', 1.0, True, None, None, None),
                        ('a2', 9.1, True, None, None, None),
                        ('c2', 8.1, True, None, None, None),
                        ('w2', 2.5, True, None, None, None))

    data = residual(org_params, x) + noise

    fit_params = Parameters()
    fit_params.add_many(('a1', 8.0, True, None, 14., None),
                        ('c1', 5.0, True, None, None, None),
                        ('w1', 0.7, True, None, None, None),
                        ('a2', 3.1, True, None, None, None),
                        ('c2', 8.8, True, None, None, None))

    fit_params.add('w2', expr='2.5*w1')

    myfit = Minimizer(residual, fit_params, fcn_args=(x,),
                      fcn_kws={'data': data})

    myfit.prepare_fit()
    out = myfit.leastsq()
    check_paras(out.params, org_params)

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
        Parameter(name='cen_g', value=9),
        Parameter(name='wid_g', value=1),
        Parameter(name='frac', value=0.50),
        Parameter(name='amp_l', expr='amp_g'),
        Parameter(name='cen_l', expr='cen_g'),
        Parameter(name='wid_l', expr='wid_g'),
        Parameter(name='line_slope', value=0.0),
        Parameter(name='line_off', value=0.0)]

sigma = 0.021  # estimate of data error (for all data points)

myfit = Minimizer(residual, pfit, # iter_cb=per_iteration,
                  fcn_args=(x,), fcn_kws={'sigma':sigma, 'data':data},
                  scale_covar=True)

myfit.prepare_fit()
init = residual(myfit.params, x)

if HASPYLAB:
    pylab.plot(x, init, 'b--')

# fit with Nelder-Mead simplex method
supported_methods = ('BFGS', 'COBYLA', 'SLSQP', 'Powell', 'Nelder-Mead')
myfit.scalar_minimize(method='Nelder-Mead')


print(' Nfev = ', myfit.nfev)
# print( myfit.chisqr, myfit.redchi, myfit.nfree)
# report_errors(myfit.params, modelpars=p_true)

fit = residual(myfit.params, x)

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
def test_constraints(with_plot=True):
    with_plot = with_plot and WITHPLOT

    def residual(pars, x, sigma=None, data=None):
        yg = gaussian(x, pars['amp_g'], pars['cen_g'], pars['wid_g'])
        yl = lorentzian(x, pars['amp_l'], pars['cen_l'], pars['wid_l'])

        model =  yg +  yl + pars['line_off'] + x * pars['line_slope']
        if data is None:
            return model
        if sigma is None:
            return (model - data)
        return (model - data) / sigma


    n = 201
    xmin = 0.
    xmax = 20.0
    x = linspace(xmin, xmax, n)

    data = (gaussian(x, 21, 8.1, 1.2) +
            lorentzian(x, 10, 9.6, 2.4) +
            random.normal(scale=0.23,  size=n) +
            x*0.5)

    if with_plot:
        pylab.plot(x, data, 'r+')

    pfit = Parameters()
    pfit.add(name='amp_g',  value=10)
    pfit.add(name='cen_g',  value=9)
    pfit.add(name='wid_g',  value=1)

    pfit.add(name='amp_tot',  value=20)
    pfit.add(name='amp_l',  expr='amp_tot - amp_g')
    pfit.add(name='cen_l',  expr='1.5+cen_g')
    pfit.add(name='wid_l',  expr='2*wid_g')

    pfit.add(name='line_slope', value=0.0)
    pfit.add(name='line_off', value=0.0)

    sigma = 0.021  # estimate of data error (for all data points)

    myfit = Minimizer(residual, pfit,
                      fcn_args=(x,), fcn_kws={'sigma':sigma, 'data':data},
                      scale_covar=True)

    myfit.prepare_fit()
    init = residual(myfit.params, x)

    result = myfit.leastsq()

    print(' Nfev = ', result.nfev)
    print( result.chisqr, result.redchi, result.nfree)

    report_fit(result.params, min_correl=0.3)

    fit = residual(result.params, x)
    if with_plot:
        pylab.plot(x, fit, 'b-')
    assert(result.params['cen_l'].value == 1.5 + result.params['cen_g'].value)
    assert(result.params['amp_l'].value == result.params['amp_tot'].value - result.params['amp_g'].value)
    assert(result.params['wid_l'].value == 2 * result.params['wid_g'].value)

    # now, change fit slightly and re-run
    myfit.params['wid_l'].expr = '1.25*wid_g'
    result = myfit.leastsq()
    report_fit(result.params, min_correl=0.4)
    fit2 = residual(result.params, x)
    if with_plot:
        pylab.plot(x, fit2, 'k')
        pylab.show()

    assert(result.params['cen_l'].value == 1.5 + result.params['cen_g'].value)
    assert(result.params['amp_l'].value == result.params['amp_tot'].value - result.params['amp_g'].value)
    assert(result.params['wid_l'].value == 1.25 * result.params['wid_g'].value)

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
	def __FitEvent(self):
		try:
			varyBlockedCurrent=True

			i0=np.abs(self.baseMean)
			i0sig=self.baseSD
			dt = 1000./self.Fs 	# time-step in ms.
			# edat=np.asarray( np.abs(self.eventData),  dtype='float64' )
			edat=self.dataPolarity*np.asarray( self.eventData,  dtype='float64' )

			blockedCurrent=min(edat)
			tauVal=dt

			estart 	= self.__eventStartIndex( self.__threadList( edat, range(0,len(edat)) ), i0, i0sig ) - 1
			eend 	= self.__eventEndIndex( self.__threadList( edat, range(0,len(edat)) ), i0, i0sig ) - 2

			# For long events, fix the blocked current to speed up the fit
			#if (eend-estart) > 1000:
			#	blockedCurrent=np.mean(edat[estart+50:eend-50])

			# control numpy error reporting
			np.seterr(invalid='ignore', over='ignore', under='ignore')

			ts = np.array([ t*dt for t in range(0,len(edat)) ], dtype='float64')

			#pl.plot(ts,edat)
			#pl.show()

			params=Parameters()

			# print self.absDataStartIndex

			params.add('mu1', value=estart * dt)
			params.add('mu2', value=eend * dt)
			params.add('a', value=(i0-blockedCurrent), vary=varyBlockedCurrent)
			params.add('b', value = i0)
			params.add('tau1', value = tauVal)

			if self.LinkRCConst:
				params.add('tau2', value = tauVal, expr='tau1')
			else:
				params.add('tau2', value = tauVal)


			optfit=Minimizer(self.__objfunc, params, fcn_args=(ts,edat,))
			optfit.prepare_fit()

			result=optfit.leastsq(xtol=self.FitTol,ftol=self.FitTol,maxfev=self.FitIters)

			# print optfit.params['b'].value, optfit.params['b'].value - optfit.params['a'].value, optfit.params['mu1'].value, optfit.params['mu2'].value
			if result.success:
				if result.params['mu1'].value < 0.0 or result.params['mu2'].value < 0.0:
					# print 'eInvalidFitParams1', optfit.params['b'].value, optfit.params['b'].value - optfit.params['a'].value, optfit.params['mu1'].value, optfit.params['mu2'].value
					self.rejectEvent('eInvalidResTime')
				# The start of the event is set past the length of the data
				elif result.params['mu1'].value > ts[-1]:
					# print 'eInvalidFitParams2', optfit.params['b'].value, optfit.params['b'].value - optfit.params['a'].value, optfit.params['mu1'].value, optfit.params['mu2'].value
					self.rejectEvent('eInvalidEventStart')
				else:
					self.mdOpenChCurrent 	= result.params['b'].value 
					self.mdBlockedCurrent	= result.params['b'].value - result.params['a'].value
					self.mdEventStart		= result.params['mu1'].value 
					self.mdEventEnd			= result.params['mu2'].value
					self.mdRCConst1			= result.params['tau1'].value
					self.mdRCConst2			= result.params['tau2'].value
					self.mdAbsEventStart	= self.mdEventStart + self.absDataStartIndex * dt

					self.mdBlockDepth		= self.mdBlockedCurrent/self.mdOpenChCurrent
					self.mdResTime			= self.mdEventEnd - self.mdEventStart
					
					self.mdRedChiSq			= result.chisqr/( np.var(result.residual) * (len(self.eventData) - result.nvarys -1) )

					# if (eend-estart) > 1000:
					# 	print blockedCurrent, self.mdBlockedCurrent, self.mdOpenChCurrent, self.mdResTime, self.mdRiseTime, self.mdRedChiSq, optfit.chisqr
					# if self.mdBlockDepth > self.BlockRejectRatio:
					# 	# print 'eBlockDepthHigh', optfit.params['b'].value, optfit.params['b'].value - optfit.params['a'].value, optfit.params['mu1'].value, optfit.params['mu2'].value
					# 	self.rejectEvent('eBlockDepthHigh')
						
					if math.isnan(self.mdRedChiSq):
						self.rejectEvent('eInvalidChiSq')
					if self.mdBlockDepth < 0 or self.mdBlockDepth > 1:
						self.rejectEvent('eInvalidBlockDepth')
					if self.mdRCConst1 <= 0 or self.mdRCConst2 <= 0:
						self.rejectEvent('eInvalidRCConstant')

					#print i0, i0sig, [optfit.params['a'].value, optfit.params['b'].value, optfit.params['mu1'].value, optfit.params['mu2'].value, optfit.params['tau'].value]
			else:
				# print optfit.message, optfit.lmdif_message
				self.rejectEvent('eFitConvergence')

		except KeyboardInterrupt:
			self.rejectEvent('eFitUserStop')
			raise
		except:
			# print optfit.message, optfit.lmdif_message
	 		self.rejectEvent('eFitFailure')

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
    def featureMap(self, DICOMImages, img_features, time_points, featuresKeys, caseLabeloutput, path_outputFolder):
        """Extracts feature maps per pixel based on request from vector of keywords featuresKeys """
        ## Retrive image data
        VOIshape = img_features['VOI0'].shape
        print VOIshape
        self.init_features(img_features, featuresKeys)
        data_deltaS=[]  
        self.allvar_F_r_i=[]
        
        # append So and to
        data_deltaS.append( 0 )  
        
        # Based on the course of signal intensity within the lesion
        So = array(img_features['VOI0']).astype(float)
        Crk = {'Cr0': mean(So)}  
        C = {}
        Carray = []
        
        # iterate point-by-point to extract feature map       
        for i in range(VOIshape[0]):
            for j in range(VOIshape[1]):
                for k in range(VOIshape[2]):
                    for timep in range(1, len(DICOMImages)):
                        pix_deltaS = (img_features['VOI'+str(timep)][i,j,k].astype(float) - img_features['VOI0'][i,j,k].astype(float))/img_features['VOI0'][i,j,k].astype(float)
                        if pix_deltaS<0: pix_deltaS=0 
                        data_deltaS.append( pix_deltaS )
                        
                        F_r_i =  array(img_features['VOI'+str(timep)]).astype(float)
                        n_F_r_i, min_max_F_r_i, mean_F_r_i, var_F_r_i, skew_F_r_i, kurt_F_r_i = stats.describe(F_r_i)
                        self.allvar_F_r_i.append(var_F_r_i)
                    
                    data = array(data_deltaS)
                    #print data                        
                    
                    # create a set of Parameters
                    params = Parameters()
                    params.add('amp',   value= 10,  min=0)
                    params.add('alpha', value= 1, min=0) 
                    params.add('beta', value= 0.05, min=0.0001, max=0.9)
                    
                    # do fit, here with leastsq self.model
                    myfit = Minimizer(self.fcn2min,  params, fcn_args=(time_points,), fcn_kws={'data':data})
                    myfit.prepare_fit()
                    myfit.leastsq()
 
                    ####################################                              
                    # Calculate R-square: R_square = sum( y_fitted - y_mean)/ sum(y_data - y_mean)
                    R_square = sum( (self.model - mean(data))**2 )/ sum( (data - mean(data))**2 )
                    #print "R^2:"
                    #print R_square
                    self.R_square_map[i,j,k] = R_square
                    
                    if 'amp' in featuresKeys:
                        amp = params['amp'].value
                        print "amp:"
                        print amp
                        self.amp_map[i,j,k] = amp
                    
                    if 'beta' in featuresKeys:
                        beta = params['beta'].value
                        self.beta_map[i,j,k] = beta
                        
                    if 'alpha' in featuresKeys:
                        alpha = params['alpha'].value
                        print "alpha:"
                        print alpha
                        self.alpha_map[i,j,k] = alpha
                        
                    if 'iAUC1' in featuresKeys:
                        iAUC1 = params['amp'].value *( ((1-exp(-params['beta'].value*t[1]))/params['beta'].value) + (exp((-params['alpha'].value+params['beta'].value)*t[1])-1)/(params['alpha'].value+params['beta'].value) )
                        print "iAUC1"
                        print iAUC1
                        self.iAUC1_map[i,j,k] = iAUC1
                        
                    if 'Slope_ini' in featuresKeys:
                        Slope_ini = params['amp'].value*params['alpha'].value
                        print "Slope_ini"
                        print Slope_ini
                        self.Slope_ini_map[i,j,k] = Slope_ini
                    
                    if 'Tpeak' in featuresKeys:
                        Tpeak = (1/params['alpha'].value)*log(1+(params['alpha'].value/params['beta'].value))
                        self.Tpeak_map[i,j,k] = Tpeak
                    
                    if 'Kpeak' in featuresKeys:
                        Kpeak = -params['amp'].value * params['alpha'].value * params['beta'].value
                        self.Kpeak_map[i,j,k] = Kpeak
                    
                    if 'SER' in featuresKeys:
                        SER = exp( (t[4]-t[1])*params['beta'].value) * ( (1-exp(-params['alpha'].value*t[1]))/(1-exp(-params['alpha'].value*t[4])) )
                        print "SER"
                        print SER
                        self.SER_map[i,j,k] = SER
                        
                    if 'maxCr' in featuresKeys:
                        print "Maximum Upate (Fii_1) = %d " %  self.maxCr
                        self.maxC_map[i,j,k] = self.maxCr
                        
                    if 'peakCr' in featuresKeys:
                        print "Peak Cr (Fii_2) = %d " %  self.peakCr
#.........这里部分代码省略.........

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
class FitModel(object):
    """base class for fitting models

    only supports polynomial background (offset, slop, quad)

    """
    invalid_bkg_msg = """Warning: unrecoginzed background option '%s'
expected one of the following:
   %s
"""
    def __init__(self, background=None, **kws):
        self.params = Parameters()
        self.has_initial_guess = False
        self.bkg = None
        self.initialize_background(background=background, **kws)

    def initialize_background(self, background=None,
                              offset=0, slope=0, quad=0):
        """initialize background parameters"""
        if background is None:
            return
        if background not in VALID_BKGS:
            print( self.invalid_bkg_msg % (repr(background),
                                          ', '.join(VALID_BKGS)))

        kwargs = {'offset':offset}
        if background.startswith('line'):
            kwargs['slope'] = slope
        if background.startswith('quad'):
            kwargs['quad'] = quad

        self.bkg = PolyBackground(**kwargs)

        for nam, par in self.bkg.params.items():
            self.params[nam] = par

    def calc_background(self, x):
        if self.bkg is None:
            return 0
        return self.bkg.calculate(x)

    def __objective(self, params, y=None, x=None, dy=None, **kws):
        """fit objective function"""
        bkg = 0
        if x is not None: bkg = self.calc_background(x)
        if y is None:     y   = 0.0
        if dy is None:    dy  = 1.0
        model = self.model(self.params, x=x, dy=dy, **kws)
        return (model + bkg - y)/dy

    def model(self, params, x=None, **kws):
        raise NotImplementedError

    def guess_starting_values(self, params, y, x=None, **kws):
        raise NotImplementedError

    def fit_report(self, params=None, **kws):
        if params is None:
            params = self.params
        return lmfit.fit_report(params, **kws)

    def fit(self, y, x=None, dy=None, **kws):
        fcn_kws={'y':y, 'x':x, 'dy':dy}
        fcn_kws.update(kws)
        if not self.has_initial_guess:
            self.guess_starting_values(y, x=x, **kws)
        self.minimizer = Minimizer(self.__objective, self.params,
                                   fcn_kws=fcn_kws, scale_covar=True)
        self.minimizer.prepare_fit()
        self.init = self.model(self.params, x=x, **kws)
        self.minimizer.leastsq()

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]

#.........这里部分代码省略.........
        # make array for data_deltaS
        data = array(data_deltaS)

        print "\n================\nMean and SE (i.e VOI sample data)"
        print mean_deltaS
        print se_deltaS

        # create a set of Parameters
        params = Parameters()
        params.add("amp", value=10, min=0)
        params.add("alpha", value=1, min=0)
        params.add("beta", value=0.05, min=0.0001, max=0.9)

        # do fit, here with leastsq model
        # define objective function: returns the array to be minimized
        def fcn2min(params, t, data):
            global model, model_res, x
            """ model EMM for Bilateral DCE-MRI, subtract data"""
            # unpack parameters:
            #  extract .value attribute for each parameter
            amp = params["amp"].value  # Upper limit of deltaS
            alpha = params["alpha"].value  # rate of signal increase min-1
            beta = params["beta"].value  # rate of signal decrease min-1

            model = amp * (1 - exp(-alpha * t)) * exp(-beta * t)

            x = linspace(0, t[4], 101)
            model_res = amp * (1 - exp(-alpha * x)) * exp(-beta * x)

            return model - data

        #####
        myfit = Minimizer(fcn2min, params, fcn_args=(t,), fcn_kws={"data": data})
        myfit.prepare_fit()
        myfit.leastsq()

        # On a successful fit using the leastsq method, several goodness-of-fit statistics
        # and values related to the uncertainty in the fitted variables will be calculated
        print "myfit.success"
        print myfit.success
        print "myfit.residual"
        print myfit.residual
        print "myfit.chisqr"
        print myfit.chisqr
        print "myfit.redchi"
        print myfit.redchi

        # calculate final result
        final = data + myfit.residual
        # write error report
        report_errors(params)

        # Calculate R-square
        # R_square = sum( y_fitted - y_mean)/ sum(y_data - y_mean)
        R_square = sum((model - mean(data)) ** 2) / sum((data - mean(data)) ** 2)
        print "R^2"
        print R_square

        self.amp = params["amp"].value
        self.alpha = params["alpha"].value
        self.beta = params["beta"].value

        ##################################################
        # Now Calculate Extract parameters from model
        self.iAUC1 = params["amp"].value * (
            ((1 - exp(-params["beta"].value * t[1])) / params["beta"].value)

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]
p_fit.add('cen_g', value=max_x)
p_fit.add('wid_g', value=2.0)
p_fit.add('frac',  value=0.50)
p_fit.add('amp_l', expr='0.5*amp_g')
p_fit.add('cen_l', value=12.5)
p_fit.add('wid_l', expr='2.5*wid_g')
p_fit.add('line_slope', value=0.0)
p_fit.add('line_off', value=0.0)

sigma = 0.041  # estimate of data error (for all data points)

myfit = Minimizer(residual, None,            # iter_cb=per_iteration,
                  fcn_args=(x,), fcn_kws={'sigma':sigma, 'data':data},
                  scale_covar=True)

myfit.prepare_fit(params=p_fit)
init = residual(p_fit, x)

if HASPYLAB:
    pylab.plot(x, init, 'b--')

myfit.leastsq()

print(' Nfev = ', myfit.nfev)
print( myfit.chisqr, myfit.redchi, myfit.nfree)

report_errors(myfit.params, modelpars=p_true)

fit = residual(myfit.params, x)

if HASPYLAB:

# 需要导入模块: from lmfit import Minimizer [as 别名]
# 或者: from lmfit.Minimizer import prepare_fit [as 别名]

#.........这里部分代码省略.........
        # make array for data_deltaS
        data = array(data_deltaS)
        
        print "\n================\nMean and SE (i.e VOI sample data)"
        print mean_deltaS
        print se_deltaS
        
        # create a set of Parameters
        params = Parameters()
        params.add('amp',   value= 10,  min=0)
        params.add('alpha', value= 1, min=0) 
        params.add('beta', value= 0.05, min=0.0001, max=0.9)
        
        # do fit, here with leastsq model
        # define objective function: returns the array to be minimized
        def fcn2min(params, t, data):
            global model, model_res, x
            """ model EMM for Bilateral DCE-MRI, subtract data"""
            # unpack parameters:
            #  extract .value attribute for each parameter
            amp = params['amp'].value    # Upper limit of deltaS
            alpha = params['alpha'].value    # rate of signal increase min-1
            beta = params['beta'].value        # rate of signal decrease min-1
                    
            model = amp * (1- exp(-alpha*t)) * exp(-beta*t)
            
            x = linspace(0, t[4], 101)
            model_res = amp * (1- exp(-alpha*x)) * exp(-beta*x)
        
            return model - data
        
        #####
        myfit = Minimizer(fcn2min,  params, fcn_args=(t,), fcn_kws={'data':data})
        myfit.prepare_fit()
        myfit.leastsq()
            
        # On a successful fit using the leastsq method, several goodness-of-fit statistics
        # and values related to the uncertainty in the fitted variables will be calculated
        print "myfit.success"
        print myfit.success
        print "myfit.residual"
        print myfit.residual
        print "myfit.chisqr"
        print myfit.chisqr
        print "myfit.redchi"
        print myfit.redchi
            
        # calculate final result
        #final = data + myfit.residual
        # write error report
        report_errors(params)
        
        # Calculate R-square
        # R_square = sum( y_fitted - y_mean)/ sum(y_data - y_mean)
        R_square = sum( (model - mean(data))**2 )/ sum( (data - mean(data))**2 )
        print "R^2"
        print R_square
        
        self.amp = params['amp'].value
        self.alpha = params['alpha'].value
        self.beta = params['beta'].value
        
        ##################################################
        # Now Calculate Extract parameters from model
        self.iAUC1 = params['amp'].value *( ((1-exp(-params['beta'].value*t[1]))/params['beta'].value) + (exp((-params['alpha'].value+params['beta'].value)*t[1])-1)/(params['alpha'].value+params['beta'].value) )
        print "iAUC1"