Home › Code Snippets

Chi2 map - calculate chi2 as a parameter varies

Posted January 29th, 2009 by andyfaff
in
Average rating
(0 votes)

Often one wants to see how chi2 varies for a given set of data, and a given fit function, for a specific parameter. The following snippet takes the name of a user defined fitfunction, a set of x,y data, a coefficient wave, the parameter of interest and a given region to examine (expressed as a percentage of the parameter you are interested in).
You can optionally specify which x,y pairs you are interested in using lhs and rhs and you can also weight the data if you have a wave containing the standard deviation.

Function chi2generator(fn,ydata,xdata,coefs,whichParam,searchArea, [edata, lhs, rhs])
	FUNCREF chi2fitfuntemplate  fn
	Wave ydata,xdata,coefs      //y data, x data, and coefficient wave
	variable whichParam, searchArea  //which of the parameters to you want to vary, how far from the original value do you want to search (%)
	Wave edata   //wave containing standard deviation for the x,y pair
	variable lhs, rhs    //specify a region of interest in the data using a left hand side and right hand side.
 
	variable originalvalue = coefs[whichparam]
	variable range = originalValue * searchArea/100
	variable ii
 
	duplicate/o ydata, :theoretical_data, :chi2_data
	Wave theoretical_data, chi2_data
 
	make/o/n=200/d chi2_map
	setscale/I x,  originalvalue - range/2, originalvalue + range/2, chi2_map
 
	if(paramisdefault(lhs))
		lhs = 0
	endif
	if(paramisdefault(rhs))
		rhs = numpnts(ydata)
	endif
 
	for(ii=0 ; ii < numpnts(chi2_map) ; ii+=1)
		coefs[whichparam] = pnt2x(chi2_map, ii)
 
		theoretical_data = fn(coefs, xdata)
		chi2_data = (ydata-theoretical_data)^2
 
		if(!paramisdefault(edata))
			chi2_data /= edata^2
		endif
		Wavestats/q/R=[lhs, rhs] chi2_data
 
		chi2_map[ii] = V_avg * V_npnts
	endfor
 
	coefs[whichparam] = originalvalue
 
	display/K=1 chi2_map
 
	killwaves/z theoretical_data, chi2_data
End
 
Function chi2fitfuntemplate(w, x)
Wave w
variable x
End

Back to top