我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用scipy.optimize.curve_fit()。
def fit_rabi(xdata, ydata): """Analyze Rabi amplitude data to find pi-pulse amplitude and phase offset. Arguments: xdata: ndarray of calibration amplitudes. length should be even. ydata: measurement amplitudes Returns: pi_amp: Fitted amplitude of pi pulsed offset: Fitted mixer offset fit_pts: Fitted points.""" #seed Rabi frequency from largest FFT component N = len(ydata) yfft = fft(ydata) f_max_ind = np.argmax(np.abs(yfft[1:N//2])) f_0 = 0.5 * max([1, f_max_ind]) / xdata[-1] amp_0 = 0.5*(ydata.max() - ydata.min()) offset_0 = np.mean(ydata) phase_0 = 0 if ydata[N//2 - 1] > offset_0: amp_0 = -amp_0 popt, _ = curve_fit(rabi_model, xdata, ydata, [offset_0, amp_0, f_0, phase_0]) f_rabi = np.abs(popt[2]) pi_amp = 0.5/f_rabi offset = popt[3] return pi_amp, offset, popt
def fit_Murn(V,E,guess=[0.0,0.0,900/RY_KBAR,1.15],lm_pars={}): """ This is the function for fitting with the Murnaghan EOS as a function of volume only. The input variable *V* is an 1D array of volumes, *E* are the corresponding energies (or other analogous quantity to be fitted with the Murnaghan EOS. *a* Note: volumes must be in a.u.^3 and energies in Rydberg. """ # reasonable initial guesses for EOS parameters if guess[0]==0.0: guess[0] = E[len(E) / 2] if guess[1]==0.0: guess[1] = V[len(V) / 2] a, pcov = curve_fit(E_MurnV, V, E, p0=guess, **lm_pars) chi = 0 for i in range(0,len(V)): chi += (E[i]-E_Murn(V[i],a))**2 return a, pcov, chi
def fit_photon_number(xdata, ydata, params): ''' Fit number of measurement photons before a Ramsey. See McClure et al., Phys. Rev. App. 2016 input params: 1 - cavity decay rate kappa (MHz) 2 - detuning Delta (MHz) 3 - dispersive shift 2Chi (MHz) 4 - Ramsey decay time T2* (us) 5 - exp(-t_meas/T1) (us), only if starting from |1> (to include relaxation during the 1st msm't) 6 - initial qubit state (0/1) ''' params = [2*np.pi*p for p in params[:3]] + params[3:] # convert to angular frequencies def model_0(t, pa, pb): return (-np.imag(np.exp(-(1/params[3]+params[1]*1j)*t + (pa-pb*params[2]*(1-np.exp(-((params[0] + params[2]*1j)*t)))/(params[0]+params[2]*1j))*1j))) def model(t, pa, pb): return params[4]*model_0(t, pa, pb) + (1-params[4])*model_0(t, pa+np.pi, pb) if params[5] == 1 else model_0(t, pa, pb) popt, pcov = curve_fit(model, xdata, ydata, p0 = [0, 1]) perr = np.sqrt(np.diag(pcov)) finer_delays = np.linspace(np.min(xdata), np.max(xdata), 4*len(xdata)) fit_curve = model(finer_delays, *popt) return popt[1], perr[1], (finer_delays, fit_curve)
def find_null_offset(xpts, powers, default=0.0): """Finds the offset corresponding to the minimum power using a fit to the measured data""" def model(x, a, b, c): return a*(x - b)**2 + c powers = np.power(10, powers/10.) min_idx = np.argmin(powers) try: fit = curve_fit(model, xpts, powers, p0=[1, xpts[min_idx], powers[min_idx]]) except RuntimeError: logger.warning("Mixer null offset fit failed.") return default, np.zeros(len(powers)) best_offset = np.real(fit[0][1]) best_offset = np.minimum(best_offset, xpts[-1]) best_offset = np.maximum(best_offset, xpts[0]) xpts_fine = np.linspace(xpts[0],xpts[-1],101) fit_pts = np.array([np.real(model(x, *fit[0])) for x in xpts_fine]) if min(fit_pts)<0: fit_pts-=min(fit_pts)-1e-10 #prevent log of a negative number return best_offset, xpts_fine, 10*np.log10(fit_pts)
def activate(self): if self.pMethod.value() == 'Poisson': v = self.pPlot.value() index = None if v != 'All': # index here if specific plot chosen index = self.pPlot.opts['limits'].index(v) - 1 plots = self.display.widget.getData(index) for plot in plots: # fit with curve_fit x, y = plot.x, plot.y if not self.guess: self.guess = _Poisson.initGuess(x, y) parameters, cov_matrix = curve_fit( _Poisson.function, x, y, self.guess) # plot poisson-deviation with fitted parameter x_vals = self._getXVals() y_vals = _Poisson.function(x_vals, *parameters) self.display.addLayer(data=(x_vals, y_vals), filename='X - Poission fitted')
def make_dirichlet(block_len, carrier_len, width=6): def _fit_model(xdata, amplitude, time_offset): xdata = np.array(xdata, dtype=np.float64) dirichlet = _dirichlet_kernel(xdata-time_offset, block_len, carrier_len) return amplitude * np.abs(dirichlet) def _interpolator(fft_mag, peak): xdata = np.array(np.arange(-(width//2), width//2+1)) ydata = fft_mag[peak + xdata] initial_guess = (fft_mag[peak], 0) popt, _ = curve_fit(_fit_model, xdata, ydata, p0=initial_guess) _, fit_offset = popt return fit_offset return _interpolator
def identify_bifurcation_point(omgp, n_splits=30): ''' Linear breakpoint model to infer drastic likelihood decrease ''' mix_m = OMGP(omgp.X, omgp.Y, K=omgp.K, kernels=omgp.kern) mix_m.variance = omgp.variance phi = omgp.phi log_liks = [] t_splits = np.linspace(mix_m.X.min(), mix_m.X.max(), n_splits) for t_s in tqdm(t_splits): mask = mix_m.X > t_s phi_mix = np.ones_like(phi) * 0.5 phi_mix[mask[:, 0]] = phi[mask[:, 0]] mix_m.phi = phi_mix log_liks.append(mix_m.log_likelihood()) x = t_splits y = np.array(log_liks) p, e = optimize.curve_fit(breakpoint_linear, x, y) return p[0]
def optimize(self): """ Applies the chosen optimization technique to the model's evaluate method and the appropriate data structures. """ # use helper method to possibly narrow data by user-entered window x, y = self.data.get_data_to_fit() # optimize based on model-specific evaluate() optimal_params, covar = sp_optimize.curve_fit( self.model.evaluate, x, y, self.model.parameters ) # store the final fit parameters and covariance in the model object for JIT calculations self.model.parameters = optimal_params self.model.covariance = covar return optimal_params, covar
def findStar(self, x, y): x = int(x) y = int(y) roi_size = Configuration.gaussian_roi_size roi = self.image[y - roi_size:y + roi_size + 1, x - roi_size:x + roi_size + 1] X, Y = np.meshgrid(range(x - roi_size, x + roi_size + 1), range(y - roi_size, y + roi_size + 1)) p0 = (1, x, y, 1, 0, 1) try: popt, _ = optimize.curve_fit(self.gaussBivarFit, (X, Y), roi.ravel(), p0=p0, maxfev=10000) except Exception as e: return 0, (0, 0), np.matrix([[0, 0], [0, 0]]), roi A, x0, y0, v1, v2, v3 = popt cov = np.matrix([[v1, v2 / 2], [v2 / 2, v3]]).I mu = (x0, y0) return A, mu, cov, roi
def dephasing(f): """ Computes the dephasing time of a given function using optical response formalisms: S. Mukamel, Principles of Nonlinear Optical Spectroscopy, 1995 About the implementation we use the 2nd order cumulant expansion. See also eq. (2) in : Kilina et al. Phys. Rev. Lett., 110, 180404, (2013) To calculate the dephasing time tau we fit the dephasing function to a gaussian of the type : exp(-0.5 * (-x / tau) ** 2) """ ts = np.arange(f.shape[0]) cumu_ii = np.stack(np.sum(f[0:i]) for i in range(ts.size)) / hbar cumu_i = np.stack(np.sum(cumu_ii[0:i]) for i in range(ts.size)) / hbar deph = np.exp(-cumu_i) np.seterr(over='ignore') popt = curve_fit(gauss_function, ts, deph)[0] xs = np.exp(-0.5 * (-ts / popt[0]) ** 2) deph = np.column_stack((deph, xs)) rate = popt[0] return deph, rate
def fit_SaddleConnGoldenL_pdf(hist_data, xdata, ydata, hrange): [sl, su] = hrange [fl, fu] = [2.3528, 8.1976] b = (fu - fl) / (su - sl) c = fu - b * su #if len(xdata) > 2 * MAX_FIT_POINTS: # trim to the middle 9000 x values # idxl, idxu = len(xdata) / 2 - MAX_FIT_POINTS/4, len(xdata) / 2 + MAX_FIT_POINTS/4 # endu = len(xdata) - MAX_FIT_POINTS / 4 - 1 # xdata = list(xdata)[idxl:idxu+1] + list(xdata)[0:MAX_FIT_POINTS/4+1] + list(xdata)[endu:] # ydata = list(ydata)[idxl:idxu+1] + list(ydata)[0:MAX_FIT_POINTS/4+1] + list(ydata)[endu:] ## # it appears that scipy has a bug, let's do an adjustment to make it go away: if len(ydata) - len(xdata) == 1: ydata = ydata[:-1] ## fit_func_v1 = lambda x, a, b, c, d: a * sc_fitfunc(b*x+c) + d fit_func = lambda x, a, d: fit_func_v1(x, a, b, c, d) p, pcov = curve_fit(fit_func, xdata, ydata)#, bounds = ([0,0], [10,10])) print p, pcov return lambda x: fit_func(x, p[0], p[1]), (1 - c) / b, [p[0], b, c, p[1]] ##
def fun(): T = [1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968] S = [29.72, 30.61, 31.51, 32.13, 32.34, 32.85, 33.56, 34.20, 34.83] xdata = np.array(T) ydata = np.log(np.array(S)) def func(x, a, b): return a + b * x # ?????????????? popt, pcov = curve_fit(func, xdata, ydata) plt.figure() plt.plot(xdata, ydata, 'ko', label="Original Noised Data") plt.plot(xdata, func(xdata, *popt), 'r', label="Fitted Curve") plt.show() plt.savefig('test.png')
def psf_zplane(stack, dz, w0, de = 1): ''' determine the position of the real focal plane. Don't mistake with psf_slice! ''' nz, ny, nx = stack.shape cy, cx = np.unravel_index(np.argmax(gf(stack,2)), (nz,ny,nx))[1:] zrange = (nz-1)*dz*0.5 zz = np.linspace(-zrange, zrange, nz) center_z = stack[:,cy-de:cy+de+1,cx-de:cx+de+1] im_z = center_z.mean(axis=2).mean(axis=1) b = np.mean((im_z[0],im_z[-1])) a = im_z.max() - b p0 = (a,0,w0,b) popt = optimize.curve_fit(gaussian, zz, im_z, p0)[0] z_offset = popt[1] # The original version is wrong return z_offset, zz
def __init__(self, A, h, alpha, CL, CD, x0=0.95): rho = 1025 super().__init__(A, h, alpha, CL, CD, rho) self.x0 = x0 # Create force functions based on input data fitParams, fitCovariances = curve_fit(liftFunc, self.alpha, self.CL) self.CL_a1 = fitParams[0] fitParams, fitCovariances = curve_fit(dragFunc, self.alpha, self.CD) self.CD_a0 = fitParams[0] self.CD_a2 = fitParams[1] # Max and min alpha self.alpha_min = 0 self.alpha_max = 30*np.pi/180
def __init__(self, A, h, alpha, CL, CD, x0=0.95): rho = 1025 self.removeBaseDrag = True super().__init__(A, h, alpha, CL, CD, rho) self.x0 = x0 # Create force functions based on input data fitParams, fitCovariances = curve_fit(liftFunc, self.alpha, self.CL) self.CL_a1 = fitParams[0] fitParams, fitCovariances = curve_fit(dragFunc, self.alpha, self.CD) self.CD_a0 = fitParams[0] self.CD_a2 = fitParams[1] # Max and min alpha self.alpha_min = 0 self.alpha_max = 40*np.pi/180 self.alpha_max = 30*np.pi/180
def setDriftData(self, alpha, CL, CDi, CM): self.alpha = alpha self.CL = CL self.CDi = CDi self.CM = CM # Calculate coefficients for lift fitParams, fitCovariances = curve_fit(liftFunc, self.alpha, self.CL) self.CL_a1 = fitParams[0] self.CL_a2 = fitParams[1] # Calculate coefficients for induced drag fitParams, fitCovariances = curve_fit(inducedDragFunc, self.alpha, self.CDi) self.CDi_a2 = fitParams[0] # Spline interpolation for moment self.CMspl = interpolate.splrep(self.alpha, self.CM, s=1e-9)
def ratingCurve(discharge, stage): """Computes rating curve based on discharge measurements coupled with stage readings. discharge = array of measured discharges; stage = array of corresponding stage readings; Returns coefficients a, b for the rating curve in the form y = a * x**b https://github.com/hydrogeog/hydro/blob/master/hydro/core.py """ from scipy.optimize import curve_fit exp_curve = lambda x, a, b: (a * x ** b) popt, pcov = curve_fit(exp_curve, stage, discharge) a = 0.0 b = 0.0 for i, j in zip(discharge, stage): a += (i - exp_curve(j, popt[0], popt[1]))**2 b += (i - np.mean(discharge))**2 r_squ = 1 - a / b return popt, r_squ
def ratingCurve(discharge, stage): """Computes rating curve based on discharge measurements coupled with stage readings. discharge = array of measured discharges; stage = array of corresponding stage readings; Returns coefficients a, b for the rating curve in the form y = a * x**b """ from scipy.optimize import curve_fit popt, pcov = curve_fit(exp_curve, stage, discharge) def r_squ(): a = 0.0 b = 0.0 for i, j in zip(discharge, stage): a += (i - exp_curve(j, popt[0], popt[1]))**2 b += (i - np.mean(discharge))**2 return 1 - a / b return popt, r_squ()
def fit(df, cell, plot=False): big_time = df.delta_t.quantile(0.75) p0 = [ 1.3, -0.38, df.adc_counts[df.delta_t >= big_time].mean(), ] try: (a, b, c), cov = curve_fit( f, df['delta_t'], df['adc_counts'], p0=p0, ) except RuntimeError: logging.error('Could not fit cell {}'.format(cell)) return np.full(4, np.nan) ndf = len(df.index) - 3 residuals = df['adc_counts'] - f(df['delta_t'], a, b, c) chisquare = np.sum(residuals**2) / ndf return a, b, c, chisquare
def fit_optimal(blocks, keys): """Fit the cumulative distribution, by optimization.""" values = list() for k in keys: values.extend([max_val for (_, _, max_val) in blocks[k]]) hist, bin_edges = np.histogram(values, bins=100, normed=True) bin_centers = np.array( [(bin_edges[i - 1] + bin_edges[i]) / 2 for i in range(1, len(bin_edges))]) cumulative = np.cumsum(hist) / np.sum(hist) popt, _ = optimize.curve_fit(func, bin_centers, cumulative) return popt
def fit_model(func, xdata, ydata, yerrdata, p0): """ Fit the provided data with the beta model. Arguments: p0: initial values for the parameters of beta model Return: (popt, infodict) popt: optimal values for the parameters infodict: * fvec: the function evaluated at the output parameters * dof: degree of freedom * chisq: chi-squared * perr: one standard deviation errors on the parameters """ popt, pcov = curve_fit(func, xdata, ydata, p0=p0, sigma=yerrdata) # the function evaluated at the output parameters fvec = lambda x: func(x, *popt) # degree of freedom dof = len(xdata) - len(popt) - 1 # chi squared chisq = np.sum(((ydata - fvec(xdata)) / yerrdata) ** 2) # one standard deviation errors on the parameters perr = np.sqrt(np.diag(pcov)) infodict = { 'fvec': fvec, 'dof': dof, 'chisq': chisq, 'perr': perr } return (popt, infodict)
def fit_drag(data, DRAG_vec, pulse_vec): """Fit calibration curves vs DRAG parameter, for variable number of pulses""" num_DRAG = len(DRAG_vec) num_seqs = len(pulse_vec) xopt_vec = np.zeros(num_seqs) perr_vec = np.zeros(num_seqs) popt_mat = np.zeros((4, num_seqs)) data = data.reshape(len(data)//num_DRAG, num_DRAG) #first fit sine to lowest n, for the full range data_n = data[1, :] T0 = 2*(DRAG_vec[np.argmax(data_n)] - DRAG_vec[np.argmin(data_n)]) #rough estimate of period p0 = [0, 1, T0, 0] popt, pcov = curve_fit(sinf, DRAG_vec, data_n, p0 = p0) perr_vec[0] = np.sqrt(np.diag(pcov))[0] x_fine = np.linspace(min(DRAG_vec), max(DRAG_vec), 1001) xopt_vec[0] = x_fine[np.argmin(sinf(x_fine, *popt))] popt_mat[:,0] = popt for ct in range(1, len(pulse_vec)): #quadratic fit for subsequent steps, narrower range data_n = data[ct, :] p0 = [1, xopt_vec[ct-1], 0] #recenter for next fit closest_ind =np.argmin(abs(DRAG_vec - xopt_vec[ct-1])) fit_range = int(np.round(0.5*num_DRAG*pulse_vec[0]/pulse_vec[ct])) curr_DRAG_vec = DRAG_vec[max(0, closest_ind - fit_range) : min(num_DRAG-1, closest_ind + fit_range)] reduced_data_n = data_n[max(0, closest_ind - fit_range) : min(num_DRAG-1, closest_ind + fit_range)] #quadratic fit popt, pcov = curve_fit(quadf, curr_DRAG_vec, reduced_data_n, p0 = p0) perr_vec[ct] = np.sqrt(np.diag(pcov))[0] x_fine = np.linspace(min(curr_DRAG_vec), max(curr_DRAG_vec), 1001) xopt_vec[ct] = x_fine[np.argmin(quadf(x_fine, *popt))] #why not x0? popt_mat[:3,ct] = popt return xopt_vec, perr_vec, popt_mat
def fit_quad(xdata, ydata): popt, pcov = curve_fit(quadf, xdata, ydata, p0 = [1, min(ydata), 0]) perr = np.sqrt(np.diag(pcov)) return popt, perr
def fit_gaussian(x, y, z_2d, save_fits=False): z = z_2d max_idx = np.unravel_index(z.argmax(), z.shape) max_row = max_idx[0] - 1 max_col = max_idx[1] - 1 z_max_row = z[max_row, :] z_max_col = z[:, max_col] A = z[max_row, max_col] p_guess_x = (A, x[max_col], 0.1*(x[-1] - x[0])) p_guess_y = (A, y[max_row], 0.1*(y[-1] - y[0])) coeffs_x, var_matrix_x = sciopt.curve_fit(gaussian, x, z_max_row, p_guess_x) coeffs_y, var_matrix_y = sciopt.curve_fit(gaussian, y, z_max_col, p_guess_y) c_x = (x[-1]-x[0])*(max_col+1)/x.size + x[0] c_y = (y[-1]-y[0])*(y.size-(max_row+1))/y.size + y[0] centre = (c_x, c_y) sigma = np.array([coeffs_x[2], coeffs_y[2]]) fwhm = 2.355 * sigma sigma_2 = 1.699 * fwhm if save_fits: with open('x_fit.dat', 'w') as fs: for c in np.c_[x, z_max_row, gaussian(x, *coeffs_x)]: s = ','.join([str(v) for v in c]) fs.write(s+'\n') with open('y_fit.dat', 'w') as fs: for c in np.c_[y, z_max_col, gaussian(y, *coeffs_y)]: s = ','.join([str(v) for v in c]) fs.write(s+'\n') return A, centre, sigma_2
def fit_Pbath(T_pts, P_pts): ''' find best fit to the P_bath function ''' # make sure T_pts and P_pts are 1d arrays npts=len(T_pts) T=np.array(T_pts).reshape(npts) P=np.array(P_pts).reshape(npts) try: ret=curve_fit(P_bath_function,T,P) except: print('insufficient data for TES %i' % TES) ret=None return ret
def myfit(x, y, betagamma, qL, orientation, T, yerr=None, krange=[]): if len(krange) != 0: # trim to given range booleanrange = [(x >= krange[0]) & (x <= krange[1])] x = x[booleanrange] y = y[booleanrange] if yerr is None: popt, pcov = curve_fit(fun, x, y, sigma=None) else: if len(krange) != 0: yerr = yerr[booleanrange] popt, pcov = curve_fit(fun, x, y, sigma=yerr, absolute_sigma=True) x = linspace(x[0], x[-1], 1e3) fit = fun(x, *popt) perr = sqrt(diag(pcov)) A = array([popt[0], perr[0]])/qL**2 B = array([popt[1], perr[1]])/qL*orientation C = array([popt[2], perr[2]]) eps, bet, alp, gam = twissfromparabola(A, B, C, T) epsn = betagamma*eps*1e6 eps *= 1e9 string = (r'$\beta=%.2f \pm %.2f$''\n' r'$\alpha=%.2f \pm %.2f$''\n' r'$\gamma=%.2f \pm %.2f$''\n' r'$\epsilon=(%.2f \pm %.2f)\pi\, nm\, rad$''\n' r'$\epsilon^*=(%.2f \pm %.2f)\pi\, \mu m\, rad$' % (bet[0], bet[1], alp[0], alp[1], gam[0], gam[1], eps[0], eps[1], epsn[0], epsn[1])) return x, fit, string
def dofit(xdata, ydata, betagamma, T, dE, D, QL, orientation, yerr=None, krange=[]): # make sure that k is positive k = abs(g)/Brho_0 xdata = abs(xdata) if len(krange) != 0: # trim to given range xdata = xdata[krange] ydata = ydata[krange] yerr = yerr[krange] if yerr is None: popt, pcov = curve_fit(fun, xdata, ydata, sigma=None) else: popt, pcov = curve_fit(fun, xdata, ydata, sigma=yerr, absolute_sigma=True) fit = fun(xdata, *popt) perr = np.sqrt(np.diag(pcov)) A = np.array([popt[0], perr[0]]) B = np.array([popt[1], perr[1]]) C = np.array([popt[2], perr[2]]) eps, bet, alp, gam = twissfromparabola(A, B, C, T, dE, D, QL, orientation) epsn = betagamma*eps*1e6 eps *= 1e9 A *= 1e9 B *= 1e9 C *= 1e9 string = (r'$function:\;y=ax^2+bx+c$''\n' r'$a=(%.2f \pm %.2f)\cdot10^{-9}$''\n' r'$b=(%.2f \pm %.2f)\cdot10^{-9}$''\n' r'$c=(%.2f \pm %.2f)\cdot10^{-9}$''\n' r'$\beta=%.2f \pm %.2f$''\n' r'$\alpha=%.2f \pm %.2f$''\n' r'$\gamma=%.2f \pm %.2f$''\n' r'$\epsilon=(%.2f \pm %.2f)\pi\, nm\, rad$''\n' r'$\epsilon^*=(%.2f \pm %.2f)\pi\, \mu m\, rad$' % (A[0], A[1], B[0], B[1], C[0], C[1], bet[0], bet[1], alp[0], alp[1], gam[0], gam[1], eps[0], eps[1], epsn[0], epsn[1])) return string, fit, xdata
def parabolafit(xdata, ydata, init, sigma): popt, pcov = curve_fit(fun, xdata, ydata, p0=init, sigma=sigma) perr = np.sqrt(np.diag(pcov)) x = np.linspace(min(xdata), max(xdata), 1e3) fit = fun(x, *popt) A = np.array([popt[0], perr[0]]) B = np.array([popt[1], perr[1]]) C = np.array([popt[2], perr[2]]) return x, fit, A, B, C
def bpm_fit(a, p): if a is None: return None x = a['channel']['data'][p][:, 1] y = a['channel']['data'][p][:, 2] ar = np.trapz(y / np.sum(y) * len(y), x) mean = np.mean(x * y / np.sum(y) * len(y)) rms = np.std(x * y / np.sum(y) * len(y)) popt, pcov = curve_fit(gaussian, x, y, p0=[ar, mean, rms]) return [popt, np.sqrt(pcov.diagonal())]
def variquad_fit(x, y): popt, pcov = curve_fit(quadratic, x, y, p0=[0, 0, 1]) return [popt, np.sqrt(pcov.diagonal())]
def _peakdetect_parabole_fitter(raw_peaks, x_axis, y_axis, points): """ Performs the actual parabole fitting for the peakdetect_parabole function. keyword arguments: raw_peaks -- A list of either the maximium or the minimum peaks, as given by the peakdetect_zero_crossing function, with index used as x-axis x_axis -- A numpy list of all the x values y_axis -- A numpy list of all the y values points -- How many points around the peak should be used during curve fitting, must be odd. return -- A list giving all the peaks and the fitted waveform, format: [[x, y, [fitted_x, fitted_y]]] """ func = lambda x, k, tau, m: k * ((x - tau) ** 2) + m fitted_peaks = [] for peak in raw_peaks: index = peak[0] x_data = x_axis[index - points // 2: index + points // 2 + 1] y_data = y_axis[index - points // 2: index + points // 2 + 1] # get a first approximation of tau (peak position in time) tau = x_axis[index] # get a first approximation of peak amplitude m = peak[1] # build list of approximations # k = -m as first approximation? p0 = (-m, tau, m) popt, pcov = curve_fit(func, x_data, y_data, p0) # retrieve tau and m i.e x and y value of peak x, y = popt[1:3] # create a high resolution data set for the fitted waveform x2 = np.linspace(x_data[0], x_data[-1], points * 10) y2 = func(x2, *popt) fitted_peaks.append([x, y, [x2, y2]]) return fitted_peaks
def _fit_func_factory(a=None, b=None, c=None, d=None, e=None): # Specific with correction if a and b and c and d and e: def func(x): """Specific Function.""" return a * ((b * (x + c))**d) + e # Specific without correction elif a and b and c and d: def func(x): """Specific Function.""" return a * ((b * (x + c))**d) # Generic, which we have scipy.optimize.curve_fit run on # scipy.optimize.curve_fit will then give us a, b, c, d, # however, scipy.optimize has touble with e. We correct e "by hand" # at the end. else: def func(x, a, b, c, d): """Unparameterized Function.""" return a * ((b * (x + c))**d) return func
def _generate_func_fit(func_factory, x, y): popt, pcov = optimize.curve_fit(func_factory(), x, y) logging.debug(popt) logging.debug(pcov) # Static shift to have the end condition be nice unshifted_func = func_factory(*popt) end_diff = y.iloc[-1] - unshifted_func(x.iloc[-1]) logging.debug("end diff: %s", end_diff) popt_shifted = np.append(popt, end_diff) return func_factory(*popt_shifted)
def fit(self, lambdav, fluxv, sigmav): if np.all(np.asarray(fluxv)>0): popt, pcov = curve_fit(self._bound_greybody, lambdav, fluxv, p0=(17,1e7), sigma=sigmav, absolute_sigma=False, maxfev=5000) return popt else: return (0,0) # ----------------------------------------------------------------------
def initialGuessNoE(lower, upper, JDp, RVp): ''' Make a guess of the values of the elements, assuming it is circular. Parameters ---------- lower : list(ndim) Lower bounds of the orbital elements. upper : list(ndim) Upper bounds of the orbital elements. JDp : list Times of observation of the primary component. RVp : list Radial velocities of the primary component. Returns ------- (K, T, P, y) : list Output from curve_fit. ''' from scipy.optimize import curve_fit # This is a simple harmonic function. def alteredNoERV(x, K, T, P, y): return K*cos((2*pi/P)*(x-T))+y return curve_fit(alteredNoERV, JDp, np.asarray(RVp), bounds=(lower, upper))[0]
def scaleSignal(img, fitParams=None, backgroundToZero=False, reference=None): ''' scale the image between... backgroundToZero=True -> 0 (average background) and 1 (maximum signal) backgroundToZero=False -> signal+-3std reference -> reference image -- scale image to fit this one returns: scaled image ''' img = imread(img) if reference is not None: # def fn(ii, m,n): # return ii*m+n # curve_fit(fn, img[::10,::10], ref[::10,::10]) low, high = signalRange(img, fitParams) low2, high2 = signalRange(reference) img = np.asfarray(img) ampl = (high2 - low2) / (high - low) img -= low img *= ampl img += low2 return img else: offs, div = scaleParams(img, fitParams, backgroundToZero) img = np.asfarray(img) - offs img /= div print('offset: %s, divident: %s' % (offs, div)) return img
def scaleParamsFromReference(img, reference): # saving startup time: from scipy.optimize import curve_fit def ff(arr): arr = imread(arr, 'gray') if arr.size > 300000: arr = arr[::10, ::10] m = np.nanmean(arr) s = np.nanstd(arr) r = m - 3 * s, m + 3 * s b = (r[1] - r[0]) / 5 return arr, r, b img, imgr, imgb = ff(img) reference, refr, refb = ff(reference) nbins = np.clip(15, max(imgb, refb), 50) refh = np.histogram(reference, bins=nbins, range=refr)[ 0].astype(np.float32) imgh = np.histogram(img, bins=nbins, range=imgr)[0].astype(np.float32) import pylab as plt plt.figure(1) plt.plot(refh) plt.figure(2) plt.plot(imgh) plt.show() def fn(x, offs, div): return (x - offs) / div params, fitCovariances = curve_fit(fn, refh, imgh, p0=(0, 1)) perr = np.sqrt(np.diag(fitCovariances)) print('error scaling to reference image: %s' % perr[0]) # if perr[0] < 0.1: return params[0], params[1]
def _fit(x, y): popt, _ = curve_fit(function, x, y, check_finite=False) return popt
def g(x,a,b): return a * np.exp(b * () popt, pcov = opt.curve_fit(f, xdata, w, p0 = (0.1, 10)); plt.figure(1); plt.clf(); plt.plot(xdata, w) plt.plot(xdata, f(xdata, *popt), 'r') plt.title(str(popt))
def sigmoid(x, x0, k): import numpy as np import pylab from scipy.optimize import curve_fit y = 1 / (1 + np.exp(-k*(x-x0))) return y
def __init__(self, data, mleValue, fitParameters=True, a=None, b=None): super(uniformModel, self).__init__(data) self.MLE = mleValue if(None in [a,b]): fitParameters = True if(fitParameters): self.b = max(self.getDataSet()) self.a = min(self.getDataSet()) def uniDist(x, a, b): return scipy.where((a<=x) & (x<=b), 1.0/float(b-a), 0.0) try: self.n, self.bins, patches = plt.hist(self.getDataSet(), self.getDatasetSize()/10, normed=1, facecolor='blue', alpha = 0.55) popt,pcov = curve_fit(uniDist,self.bins[:-1], self.n, p0=[self.a, self.b]) ##plt.plot(bins[:-1], gaus(bins[:-1],*popt),'c-',label="Gaussian Curve with params\na=%f\nx0=%f\nsigma=%f" % (popt[0], popt[1], popt[2]), alpha=0.5) print "Fitted uniform distribution pdf curve to data with params a %f, b %f" % (popt[0], popt[1]) self.a = popt[0] self.b = popt[1] ##self.sigma = popt[2] self.fitted = True except RuntimeError: print "Unable to fit data to uniform distribution pdf curve" raise except Warning: raise RuntimeError else: self.a = a self.b = b
def __init__(self, data, mleValue, fitParameters=True, mean=None, sigma=None): super(normalModel, self).__init__(data) self.MLE = mleValue if(None in [mean, sigma]): fitParameters=True if(fitParameters): mean = Mean(self.getDataSet()) sigma = standardDeviation(self.getDataSet()) try: def normDist(x, x0, sigma): output = 1.0/sqrt(2*np.pi*(sigma**2)) output *= exp(-0.5*((x - x0)/sigma)**2) return output self.n, self.bins, patches = plt.hist(self.getDataSet(), self.getDatasetSize()/10, normed=1, facecolor='blue', alpha = 0.55) popt,pcov = curve_fit(normDist,self.bins[:-1], self.n, p0=[mean, sigma]) ##plt.plot(bins[:-1], gaus(bins[:-1],*popt),'c-',label="Gaussian Curve with params\na=%f\nx0=%f\nsigma=%f" % (popt[0], popt[1], popt[2]), alpha=0.5) print "Fitted gaussian curve to data with params x0 %f, sigma %f" % (popt[0], popt[1]) self.x0 = popt[0] self.sigma = popt[1] self.fitted = True except RuntimeError: print "Unable to fit data to normal curve, runtime error" raise except Warning: raise RuntimeError else: self.x0 = mean self.sigma = sigma
def __init__(self, data, mleValue, fitParameters=True, mean=None): super(exponentialModel, self).__init__(data) self.MLE = mleValue if(None in [mean]): fitParameters = True if(fitParameters): mean = Mean(self.getDataSet()) try: def expDist(x, x0): return (exp(-(x/x0))/x0) self.n, self.bins, patches = plt.hist(self.getDataSet(), self.getDatasetSize()/10, normed=1, facecolor='blue', alpha = 0.55) popt,pcov = curve_fit(expDist,self.bins[:-1], self.n, p0=[mean]) ##plt.plot(bins[:-1], gaus(bins[:-1],*popt),'c-',label="Gaussian Curve with params\na=%f\nx0=%f\nsigma=%f" % (popt[0], popt[1], popt[2]), alpha=0.5) print "Fitted gaussian curve to data with params x0 %f" % (popt[0]) self.x0 = popt[0] ##self.sigma = popt[2] self.fitted = True except RuntimeError: print "Unable to fit data to exponential curve" raise except Warning: raise RuntimeError else: self.x0 = mean
def calcslope(self, fmin=1./1000.0, fmax=1./2.0): """ Measures the slope of the PS, between fmin and fmax. All info about this slope is sored into the dictionnary self.slope. Should this depend on self.flux ?? No, only the constructor depends on this ! This is just fitting the powerspectrum as given by the constructor. """ if fmin == None: fmin = self.f[1] if fmax == None: fmax = self.f[-1] reg = np.logical_and(self.f <= fmax, self.f >= fmin) fitx = np.log10(self.f[reg]).flatten() fity = np.log10(self.p[reg]).flatten() if not np.all(np.isfinite(fity)): print "Skipping calcsclope for flat function !" return self.slope = {} def func(x, m, h): return m*x + h popt = spopt.curve_fit(func, fitx, fity, p0=[0.0, 0.0])[0] sepfit = func(fitx, popt[0], popt[1]) # we evaluate the linear fit. self.slope["slope"] = popt[0] self.slope["f"] = 10.0**fitx self.slope["p"] = 10.0**sepfit self.slope["fmin"] = fmin self.slope["fmax"] = fmax
def _create_initial_peak_position(self, frequencies, sf, prec=1e-12): position = frequencies[np.argmax(sf)] # "curve_fit" does not work well for extremely small initial guess. # To avoid this problem, "position" is rounded. # See also "http://stackoverflow.com/questions/15624070" if abs(position) < prec: position = 0.0 return position
def fit_psychometric(xdata, ydata, func=None, p0=None): """ Fit a psychometric function. """ if func is None: func = 'cdf_gaussian' if p0 is None: if func == 'cdf_gaussian': p0 = [np.mean(xdata), np.std(xdata)] elif func == 'cdf_gaussian_with_guessing': p0 = [np.mean(xdata), np.std(xdata), 0.1] else: raise ValueError("[ pycog.fittools.fit_psychometric ] Need initial guess p0.") if isinstance(func, str): func = fit_functions[func] popt_list, pcov_list = curve_fit(func, xdata, ydata, p0=p0) # Return parameters with names args = inspect.getargspec(func).args popt = OrderedDict() for name, value in zip(args[1:], popt_list): popt[name] = value return popt, func
def expfit(x, p1, p2): return p1 * np.exp(p2*x) # popt, pcov = curve_fit(func, xdata, ydata) # popt = optimized parameters matrix (in the order they ar in your func) # pcov = covariance stats # p0 = (p1,p2,...) - initial guesses at the coeff, optional
def normalized_cut_values(self, N=10): f_res, f_asph, _estimated = self.inert_fractions() cuts = list(self.culled_cuts()) if len(cuts) == 0: if self.record.api is not None: oil_api = self.record.api else: oil_rho = self.density_at_temp(288.15) oil_api = est.api_from_density(oil_rho) BP_i = est.cut_temps_from_api(oil_api) fevap_i = np.cumsum(est.fmasses_flat_dist(f_res, f_asph)) else: BP_i, fevap_i = zip(*[(c.vapor_temp_k, c.fraction) for c in cuts]) popt, _pcov = curve_fit(_linear_curve, BP_i, fevap_i) f_cutoff = _linear_curve(732.0, *popt) # center of asymptote (< 739) popt = popt.tolist() + [f_cutoff] fevap_i = np.linspace(0.0, 1.0 - f_res - f_asph, (N * 2) + 1)[1:] T_i = _inverse_linear_curve(fevap_i, *popt) fevap_i = fevap_i.reshape(-1, 2)[:, 1] T_i = T_i.reshape(-1, 2)[:, 0] above_zero = T_i > 0.0 T_i = T_i[above_zero] fevap_i = fevap_i[above_zero] return T_i, fevap_i
