function [ outStruct ] = ExponentialFit(x,y,sigma_y,tol) %ExponentialFit: Performs a nonlinear weighted fit on data and calculates % uncertainties on the fit parameters [A;xi] % y = A*exp(-(x-x0)/xi) % % % Note: Be wary of offset exponentials. If the tail of the distribution % does not decay to zero, the fit will not work properly. In this % case it is recommended that the user either modify the code to % add parameters for the background (e.g. an exponential plus a % constant) or that the user first subtract the background noise % and then fit the data near the peak. % % Note: This program uses Gauss-Newton iteration to determine the fit % parameters, which for some cases might not converge. In this % case, it is recommended that the tail of the distribution be % cut off, with most of the data kept near the peak. % % % INPUTS: x, y, (sigma_y), (tol) % x: Independant variable (for an exponential distribution) % y: Dependant variable (counts or intensity of x) % sigma_y: Uncertainties on the dependant variable % tol: Relative tolerance for the numerical iteration % % The inputs x, y, and possibly sigma_y must be the same size and 1xN % or Nx1. If sigma_y is not provided, the fit will assume equal % weighting and will solve for the sigma_y that makes the reduced % chi-squared equal to 1. If tol is not provided, the iteration will % use a preset tolerance for its stop criterion. % % % OUTPUTS: A structure array containing the following fields: % x0: Smallest value of the independent variable % A: Amplitude of the exponential % xi: Decay constant of the exponential % sig_[var]: Uncertainties of the fit parameters % corr_[vars]: Correlation coefficients of the fit parameters % fit_x, fit_y: The x and y coordinates of the fit % chi2: The chi-squared estimate of goodness-of-fit % chi2red: Reduced chi-squared % sigma_y: If not an input, estimates the uncertainty of y % % For an equal-weighted fit, a sigma_y will be found such that the % reduced chi-squared will be 1. x = x(:); y = y(:); %Resizes as column arrays %Initial guess of fit parameters x0 = min(x); xi0 = sum(x.*y)/sum(y) - x0; A0 = max(y); %Shifts exponential to improve numerics xn = (x-x0)/xi0; yn = y/A0; P = [1;1]; %Estimates of shifted [Amplitude; Decay Constant] if (~exist('sigma_y','var')) %No uncertainties given (unweighted fit) u = 1; sigma_yn = ones(size(y)); %Temporarily sets all uncertainties to 1 %This is rescaled once sigma_y is estimated else u = 0; sigma_y(sigma_y<=0) = min(sigma_y(sigma_y>0)); %Resets nonpositive %uncertainties sigma_yn = sigma_y/A0; end if (~exist('tol','var')) %No relative tolerance given tol = 1e-12; end %Weight matrix W = diag(sigma_yn.^-2); %Exponential and Jacobian derivative with respect to fit parameters f = @(x,x0,A,xi)(A*exp(-(x-x0)/xi)); Df = @(x,x0,A,xi)([f(x,x0,A,xi)/A,f(x,x0,A,xi).*(x-x0)/xi^2]); %Gauss-Newton iteration err = Inf; while err > tol J = Df(xn,0,P(1),P(2)); dyn = yn-f(xn,0,P(1),P(2)); N = J'*W*J; c = J'*W*dyn; dP = N\c; P = P+dP; err = sqrt((c'*dP)/(dyn'*W*dyn)); end %Unscales A = P(1)*A0; xi = P(2)*xi0; df = length(x)-length(P); %Degrees of freedom dy = y-f(x,x0,A,xi); J = Df(x,x0,A,xi); N = J'*W*J; if u==1 %Unweighted fit; W (and hence N) is scaled to match uncertainties chi2red = 1; chi2 = df; s_y = sqrt((dy'*dy)/df); %Estimates the uncertainty in y V = N\eye(size(N))*s_y^2; %Matrix of variances and covariances else %Weighted fit; W (and hence N) is unscaled here chi2 = dy'*W*dy/A0^2; chi2red = chi2/df; V = N\eye(size(N))*A0^2; %Matrix of variances and covariances end outStruct = struct(); outStruct.x0 = x0; outStruct.A = A; outStruct.sig_A = sqrt(V(1,1)); outStruct.xi = xi; outStruct.sig_xi = sqrt(V(2,2)); outStruct.corr_A_xi = V(1,2)/sqrt(V(1,1)*V(2,2)); outStruct.fit_x = x; outStruct.fit_y = f(x,x0,A,xi); outStruct.chi2 = chi2; outStruct.chi2red = chi2red; if u==1 outStruct.sigma_y = s_y; end end %Fine structure at high energies