Source code for fits

''' Collection of fitting functions '''
"""

Author : Thomas Haslwanter
Date : Oct-2013
Version: 1.1
"""

from skimage import filter
import numpy as np
import os
import matplotlib.pyplot as plt
import scipy.stats as stats
from collections import namedtuple
import math
 

[docs]def demo_ransac(): '''Find the best-fit circle in an image, using the RANSAC algorithm ''' debug_flag = 1 # Since this function is only used in demo_ransac, define it here def drawCircle(center,r): '''Draw a circle''' nPts = 100 phi = np.linspace(0,2*np.pi,nPts) x = center[0] + r*np.cos(phi) y = center[1] + r*np.sin(phi) plt.hold(True) plt.plot(y,x,'r') # Get the data os.chdir(r'C:\Users\p20529\Coding\Matlab\imProc\ransac_fitCircle') data = plt.imread('0_5.bmp') # Eliminate edge artefacts rim = 20 data = data[rim:-rim, rim:-rim] imSize = data.shape # Find edges edges = filter.sobel(data) edgePnts = edges>0.15 np.sum(edgePnts) [x,y] = np.where(edgePnts==True) # set RANSAC parameters par={'eps': 3, 'ransac_threshold': 0.1, 'nIter': 500, 'lowerRadius': 5, 'upperRadius': 200 } # Allocate memory, for center(2D), radius, numPoints (structured array) fitted = np.zeros(par['nIter'], \ dtype={'names':['center', 'radius', 'nPts'], 'formats':['2f8', 'f8', 'i4']}) for ii in range(par['nIter']): # Takes 3 random points, and find the corresponding circle randEdges = np.random.permutation(len(x))[:3] (center, radius) = fitCircle(x[randEdges], y[randEdges]) # Eliminate very small and very large circles, and centers outside the image if not (par['lowerRadius'] < radius < par['upperRadius'] and \ 0 <= center[0] < imSize[0] and \ 0 <= center[1] < imSize[1]): continue # Make sure a reasonable number of points lies near that circle centerDistance = np.sqrt((x-center[0])**2 + (y-center[1])**2) inCircle = np.where(np.abs(centerDistance-radius)<par['eps'])[0] inPts = len(inCircle) if inPts < par['ransac_threshold'] *4*np.pi*radius*par['eps'] or inPts < 3: continue # Fit a circle to all good points, and save the corresponding parameters (center, radius) = fitCircle(x[inCircle], y[inCircle]) fitted[ii] = (center, radius, inPts) # If you want to see these points: if debug_flag == 1: plt.plot(y,x,'.') plt.hold(True) plt.plot(y[inCircle], x[inCircle],'r.') plt.plot(y[randEdges], x[randEdges], 'g+', ms=15) plt.axis('equal') plt.show() # Sort the circles, according to number of points included fitted = np.sort(fitted,order='nPts') # Show the best-fitting circle plt.imshow(data, cmap='gray', origin='lower') drawCircle(fitted[-1]['center'], fitted[-1]['radius']) plt.show()
[docs]def fit_circle(x,y): ''' Determine the best-fit circle to given datapoints. Parameters ---------- x : array (N,) x-values. y : array (N,) corresponding y-values. Returns ------- center : array (2,) x/y coordinates of center of the circle radius : float Circle radius. ''' M = np.vstack((2*x,2*y,np.ones(len(x)))).T (par,_,_,_) = np.linalg.lstsq(M,x**2+y**2) center = par[:2] radius = np.sqrt(par[2]+np.sum(center**2)) return(center, radius)
[docs]def fit_exp(tFit, yFit): ''' Calculates best-fit parameters for the exponential decay to an offset. Parameters ---------- tFit : array (N,) Time values. yFit : array (N,) Function values Returns ------- offset : float Function offset/bias. amp : float Amplitude of exponential function tau : float Decay time. ''' from scipy import optimize # Define the fit-function and the error-function fitfunc = lambda p, x: p[0] + p[1]*np.exp(-x/p[2]) errfunc = lambda p,x,y: fitfunc(p,x) - y pInit = np.r_[0, 1, 1] # Initial values # Make the fit pFit, success = optimize.leastsq(errfunc, pInit, args=(tFit, yFit)) # Plot the data and the fit plt.plot(tFit, yFit, label='rawdata') plt.hold(True) plt.plot(tFit, fitfunc(pFit,tFit), label='fit') plt.legend() plt.show() ExpFit = namedtuple('ExpFit', ['offset', 'amplitude', 'tau']) return ExpFit(pFit)
[docs]def fit_line(x, y, alpha=0.05, newx=[], plotFlag=1): """ Linear regression fit. Parameters ---------- x : ndarray Input / Predictor. y : ndarray Input / Estimator. alpha : float Confidence limit [default=0.05] newx : float or ndarray Values for which the fit and the prediction limits are calculated (optional) plotFlag: int, optional 1 = plot, 0 = no_plot [default] Returns ------- a : float Intercept b : float Slope ci : ndarray Lower and upper confidence interval for the slope info : dictionary contains return information on - residuals - var_res - sd_res - alpha - tval - df newy : list(ndarray) Predictions for (newx, newx-ciPrediction, newx+ciPrediction) Examples -------- >>> import numpy as np >>> from fitLine import fitLine >>> x = np.r_[0:10:11j] >>> y = x**2 >>> (a,b,(ci_a, ci_b),_)=fitLine(x,y) Notes ----- Example data and formulas are taken from D. Altman, "Practical Statistics for Medicine" """ # Summary data n = len(x) # number of samples Sxx = np.sum(x**2) - np.sum(x)**2/n # Syy = np.sum(y**2) - np.sum(y)**2/n # not needed here Sxy = np.sum(x*y) - np.sum(x)*np.sum(y)/n mean_x = np.mean(x) mean_y = np.mean(y) # Linefit b = Sxy/Sxx a = mean_y - b*mean_x # Residuals fit = lambda xx: a + b*xx residuals = y - fit(x) var_res = np.sum(residuals**2)/(n-2) sd_res = np.sqrt(var_res) # Confidence intervals se_b = sd_res/np.sqrt(Sxx) se_a = sd_res*np.sqrt(np.sum(x**2)/(n*Sxx)) df = n-2 # degrees of freedom tval = stats.t.isf(alpha/2., df) # appropriate t value ci_a = a + tval*se_a*np.array([-1,1]) ci_b = b + tval*se_b*np.array([-1,1]) # create series of new test x-values to predict for npts = 100 px = np.linspace(np.min(x),np.max(x),num=npts) se_fit = lambda x: sd_res * np.sqrt( 1./n + (x-mean_x)**2/Sxx) se_predict = lambda x: sd_res * np.sqrt(1+1./n + (x-mean_x)**2/Sxx) print('Summary: a={0:5.4f}+/-{1:5.4f}, b={2:5.4f}+/-{3:5.4f}'.format(a,tval*se_a,b,tval*se_b)) print('Confidence intervals: ci_a=({0:5.4f} - {1:5.4f}), ci_b=({2:5.4f} - {3:5.4f})'.format(ci_a[0], ci_a[1], ci_b[0], ci_b[1])) print('Residuals: variance = {0:5.4f}, standard deviation = {1:5.4f}'.format(var_res, sd_res)) print('alpha = {0:.3f}, tval = {1:5.4f}, df={2:d}'.format(alpha, tval, df)) # Return info ri = {'residuals': residuals, 'var_res': var_res, 'sd_res': sd_res, 'alpha': alpha, 'tval': tval, 'df': df} if plotFlag == 1: # Plot the data plt.figure() plt.plot(px, fit(px),'k', label='Regression line') #plt.plot(x,y,'k.', label='Sample observations', ms=10) plt.plot(x,y,'k.') x.sort() limit = (1-alpha)*100 plt.plot(x, fit(x)+tval*se_fit(x), 'r--', lw=2, label='Confidence limit ({0:.1f}%)'.format(limit)) plt.plot(x, fit(x)-tval*se_fit(x), 'r--', lw=2 ) plt.plot(x, fit(x)+tval*se_predict(x), '--', lw=2, color=(0.2,1,0.2), label='Prediction limit ({0:.1f}%)'.format(limit)) plt.plot(x, fit(x)-tval*se_predict(x), '--', lw=2, color=(0.2,1,0.2)) plt.xlabel('X values') plt.ylabel('Y values') plt.title('Linear regression and confidence limits') # configure legend plt.legend(loc=0) leg = plt.gca().get_legend() ltext = leg.get_texts() plt.setp(ltext, fontsize=10) # show the plot plt.show() if newx != []: try: newx.size except AttributeError: newx = np.array([newx]) print('Example: x = {0}+/-{1} => se_fit = {2:5.4f}, se_predict = {3:6.5f}'\ .format(newx[0], tval*se_predict(newx[0]), se_fit(newx[0]), se_predict(newx[0]))) newy = (fit(newx), fit(newx)-se_predict(newx), fit(newx)+se_predict(newx)) return (a,b,(ci_a, ci_b), ri, newy) else: return (a,b,(ci_a, ci_b), ri)
[docs]def fit_sin(tList,yList,freq): ''' Fit a sine wave with a known frequency to a given set of data. y = amplitude * sin(2*pi*freq * tList + phase*pi/180) + bias Parameters ---------- yList : array datapoints tList : float time base, in sec freq : float in Hz Returns ------- phase : float in degrees amplitude : float bias : float Examples -------- >>> (phase, amp, offset) = thLib.fitSin.fitSine(t, data, freq) ''' # Get input data b = yList rows = [ [np.sin(freq*2*np.pi*t), np.cos(freq*2*np.pi*t), 1] for t in tList] A = np.array(rows) # Make the fit (w,residuals,rank,sing_vals) = np.linalg.lstsq(A,b) # Extract desired parameters ... phase = math.atan2(w[1],w[0])*180/np.pi amplitude = np.linalg.norm([w[0],w[1]],2) bias = w[2] # ... and return them SinFit = namedtuple('SinFit', ['phase', 'amp', 'offset']) return SinFit(phase,amplitude,bias)
[docs]def regress(x, y, alpha=0.05): ''' Linear regression and confidence intervals, for a linear regression y = k * x + d (kd, ci) = regress(x,y,[alpha=0.05]) Parameters ---------- x : ndarray (N,) predictors y : ndarray(N,) responses alpha: float Defines the 100*(1-alpha)% confidence level in ci. Returns ------- k : float estimated slope d : float estimated intercept ci : ndarray (2,) confidence intervals for the slope Note ---- Similar to the MATLAB command "regress" Examples -------- >>> x = arange(100) >>> y = 0.3*x + 10 + randn(len(x)) >>> (kd, ci) = thLib.fits.regress(x,y) See also -------- fit_line ''' regressed = stats.linregress(x,y) kd = regressed[0:2] se_k = regressed[4] level = (1.-alpha/2.) tVal = stats.t.ppf(level,len(x)-2) ci = kd[0] + se_k*tVal*np.array([-1, 1]) RegressionLine = namedtuple('Regression', ['slope', 'intercept', 'ci_slope']) return RegressionLine(kd[0], kd[1], ci)
if __name__=='__main__': print ('Done')