Many sites have the solution to the case where one wants to fit an exponential of the form A e ^ {B x}. In general one can just take the log of the function and fit the straight line that results. What about data where you don’t know the y-intersect, i.e. of the form: A + B e ^ {C x}. We are trying to guess A, so we can’t just subtract it.
Skipping the detail, wrote a matlab script to solve this, some people might find it useful.
function [numresult,CIinfos] = fit(xdata,ydata,stdvar) %result = getfit(xdata,ydata,variance) % xdata = xcoordinates % ydata = ycoordinates % stdvar = variance of data points % ALL COLUMN VECTORS PLEASE! % % numresult gives in the first column: estimates of the parameters (e.g. intercept and gradient) % second and third column: estimates of the parameters' standard error confidence intervals % CIinfos gives the parameter values at the limits of the confidence intervals, for each parameter. % % There's plenty of things that need tweaking in the code, it's not plug-and-play. % % e.g. mess with rangesize and rangeres to decide on what values to look at when deciding on the CIs. % % three functions are available, comment out the one you want: % MEAN - one parameter (the mean) % LINEAR - two parameters (intersept & gradient) % EXPONENT - three parameters, A, B and C in A+Be^{Cx} options = optimset('TolX',0.05); if (size(xdata,1)<size(xdata,2)) error('Please use column vectors, not row vectors'); end %MEAN %Nparams = 1; %func = @(a)(-sum( log( normpdf((ydata-a(1))./sqrt(stdvar),0,1)./sqrt(stdvar) ) )); %[numresult,res]=fminsearch(func,[1]); %LINEAR %Nparams = 2; %func = @(a)(-sum( log( normpdf((ydata-(a(1)+a(2)*xdata))./sqrt(stdvar),0,1)./sqrt(stdvar) ) )); %[numresult,res]=fminsearch(func,[1,1]); %EXPONENT Nparams = 3; %more likely to give Inf (as log(norm(X)) depends on norm giving non-zero results). %func = @(a)(-sum( log( normpdf((ydata-( a(1)+a(2)*exp(a(3)*xdata)) )./sqrt(stdvar),0,1)./sqrt(stdvar) ) )); %this will be more robust... temp = @(a)( ydata-( a(1)+a(2)*exp(a(3)*xdata)) ); func = @(a)(sum( temp(a).^2 ./ (2*stdvar) ) ); [result,res]=fminsearch(func,[max(ydata),-max(ydata),-.5]); result = result'; numresult = result; numresult for param=1:Nparams fprintf('Parameter: %i\n',param); testa = result; %set testa to the numerical best-fit result ps = []; rangesize = 8+2*abs(result(param)); %assume that the range we should look over is roughly of the order of the parameter... rangeres = 3000; paramvalrange = result(param)-rangesize:rangesize/rangeres:result(param)+rangesize; %range of values we'll try counter = 0; for paramval = paramvalrange counter=counter+1; if (mod(counter,10)==0); fprintf('%i%% ',round(100*counter/rangeres)); end; testa(param,1) = paramval; %find the best result while fixing this parameter... [temp,x] = fminsearch(@(a) func([a(1:param-1);paramval;a(param+1:end)]),testa,options); %not necessarily the best result: the sample distribution may not be symmetric %ie, this assumes the two parameters are independent! % x = func(testa); ps(end+1) = -x; end ps = exp(ps); if (sum(ps)==0) subplot(2,1,1); plot(ps); subplot(2,1,2); plot(xdata,ydata,'x'); fprintf('Best guess: '); error('Probabilities all zero'); end ps = ps./sum(ps); cumps = cumsum(ps); calcCI = [paramvalrange(find(diff(cumps>0.1587))),paramvalrange(find(diff(cumps>0.8413)))]; %find the parameters that fit this CI limit the best (just run the minsearch again) for ci_i=1:2 [CIinfo] = fminsearch(@(a) func([a(1:param-1);calcCI(ci_i);a(param+1:end)]),testa,options); CIinfos{param,ci_i} = CIinfo; end numresult(param,2) = calcCI(1); numresult(param,3) = calcCI(2); if (1) subplot(Nparams,1,param); plot(paramvalrange,ps); grid; title(sprintf('Parameter %i',param)); hold on; plot([calcCI(1),calcCI(1)],[0,max(ps)],'k-'); plot([calcCI(2),calcCI(2)],[0,max(ps)],'k-'); end end figure; errorbar(xdata,ydata,sqrt(stdvar),'k+-'); hold on; xs=[min(xdata):0.1:max(xdata)]; ys=result(1)+result(2)*exp(result(3)*xs); plot(xs,ys,'b-'); for param=1:3 for ci_i=1:2 CIinfos{param,ci_i}(param,1) = numresult(param,ci_i+1); nr = CIinfos{param,ci_i}; ys=nr(1,1)+nr(2,1)*exp(nr(3,1)*xs); plot(xs,ys,'b:'); end end
There is also a test function, which tests that 69% of the test data is within the CIs the fitting function returns:
%make some test data parameters xdata=[1,2,3,4,5,6,7,8,9,10]'; a=[10,-10,-.4]; %parameters for curve stdvar = abs(rand(size(xdata)))+1; %noisiness of data points stdvar = stdvar./10; guesses = []; rs = []; for it=1:1000 %Linear: % ydata = a(1)+a(2)*xdata + randn(size(xdata)).*sqrt(stdvar); %Exponential: ydata = a(1)+a(2)*exp(a(3)*xdata) + randn(size(xdata)).*sqrt(stdvar); r = getfit_numerical(xdata,ydata,stdvar); rs(it,:) = r(:,1)'; for param = 1:size(r,1) guesses(it,param) = (r(param,2)<a(param)) & (r(param,3)>a(param)); end end mean(guesses)
For the mean and straight-line solutions the probability distributions of the parameters are
Parameter confidence intervals: parameter 1 is the asymptote (if it is an exponential decay), parameter 2 is the height of the curve (from x axis to asymptote), parameter 3 is the exponential curve’s parameter.
(as we already know) normal. For the exponential fit however, they become very non-normal with very long tails. This seemed like a bug at first. But if one considers the extreme of the first parameter: Very large values can be ‘cancelled out’ by very negative values of the 2nd parameter (and a smaller value of the 3rd parameter). This would show up as something closer to a straight line through the points. Plotting the 6 limits to the 3 parameters (and selecting the maximum-likelihood for each other pair of parameters) reveals this. In the figure below the dotted lines illustrate these extremes: Note that a reasonable fit is achieved by lines extending from an almost straight line to highly curved (almost-vertical-then-horizontal) solutions.
Hope this comes in useful!
Fitting exponential curve to a set of points (each with an associated error value). The CIs indicate 1 standard error, the 6 dotted lines are the 3 pairs of CIs for the 3 parameters.
lionfish 