CPlotView:Plotting a Least Squares Fit

The following plot shows the result of a least squares fit to 10 points. Each point is the local group average of several points and the standard deviation among the points is shown as its error bar. The least square fit used optimal weighting in which the weight of each point is equal to the inverse of its variance. This plot shows 2 series in Overplot mode: Series 1 plots the square markers with error bars. Series 2 plots the polynomial fit evaluated at 100 points over the x range of the data. The script that produced this plot is listed below.

Sample Script

The script first enters the data to be plotted. In this case, the values were typed-in, but the data might also have been loaded form a text file or from a Mira table. Each of the 10 entries represents the group average and standard deviation of the grouped data. Three entries are given for x, y, and s. The value of s is the y direction standard deviation for the group of points at average coordinates x, y. The script defines a simple helper function to convert the standard deviation s into a weight using the inverse variance method and then uses it to compute the cubic polynomial fit. The data are then plotted in Series 1 and and the fit is evaluated at 100 points over the data range to create Series 2.

L = new_lsqfit()	-- create the least squares fitting object
L:SetNumCoefs( 4 )	-- fit a 4 term (cubic) polynomial
L:AddPtWt( -4.5, 2, 3.5 )	-- add 10 data (x,y) points with y weights
L:AddPtWt( 5.2, 20, 1.4 )
L:AddPtWt( 10.1, 35.5, 3.8 )
L:AddPtWt( 13, 50, 1.0 )
L:AddPtWt( 18.7, 65, 3.0 )
L:AddPtWt( -3, 7, 1.0 )
L:AddPtWt( 4, 25, 2.5 )
L:AddPtWt( 9.7, 40, 5.0 )
L:AddPtWt( 11, 55, 0.5 )
L:AddPtWt( 14.5, 70, 6.5 )
L:Fit()	-- compute the least squares fit
-- now plot the data points and their fit
P = new_plotview()	-- create a new CPlotView instance
-- Plot the data points as markers
for i = 1, L:GetNumPts() do	-- add each point as a marker
P:Add( L:GetPtX(i)[1], L:GetPtY(i), 0, L:GetPtWt(i), 0 )
end
-- make a caption from the fit sigma and create the plot of data points:
sCap = Sprintf("Fit StdDev = %lg", L:GetSigmaFit() )
P:PlotPoints( "X Value", "Observations", "Fit to Data", sCap )
P:Empty()	-- clear out existing plot points
-- Plot the fit using 100 line segments over the range of data:
nMinX, nMaxX = L:GetRangeX()
nDeltaX = (nMaxX - nMinX) / 100
for x = nMinX, nMaxX, nDeltaX do	-- run through the x range in 100 steps
P:Add( x, L:Eval(x), 0, 0, 2 )	-- plot a line between points
end
P:AddSeries("Cubic Fit")	-- add evaluation points as a new series
P:SetSeriesMode(4)	-- show both series using overplot mode

Sample Script

Related Topics