least squares fit for a histogram

Hi,

I've created a histogram and need to least-squares fit it to a function with the form N(t)=Bexp(-t/τ)+A . I can't use the histfit function since I don't have the statistics toolbox but I don't know what else to try or even where to start.

Thanks in advance!

Answer 1

[Matt J]

FMINSPLEAS ( available here ) should work well on this problem, and doesn't require any toolboxes. It should work well because N(t) has a linear dependence on all but one of your parameters, and you can tell it to take advantage of this.

Answer 2

[Star Strider]

It’s easy enough to do with fminsearch and two other lines of code, one being your objective function and the second a sum-squared-error function:

t = linspace(0,1,10);                       % Create Bin Data
H = 3.*exp(-t/5) + 7 + 0.01*rand(size(t));  % Create Histogram Data

N = @(b,t) b(1).*exp(-t./b(2)) + b(3);      % Objective Function
SSE = @(b) sum((H - N(b,t)).^2);            % Cost Function (Sum Squared Error)

[b, fv] = fminsearch(SSE, [1; 1; 1]);

fprintf(1, '\n\tB   = %.5f\n\ttau = %.5f\n\tA   = %.5f\n\n',b)

I named the histogram data ‘H’ rather than ‘N’ here because it is advisable to have different names for functions and variables to avoid confusion. Change the variable name in the ‘SSE’ function to fit your needs. You may also have to change the initial parameter estimates (here the [1; 1; 1] vector) if the routine has problems fitting your data. Nonlinear parameter estimation sometimes requires some experimentation with the initial estimates.

[Matt J]

FMINSPLEAS will be more reliable than direct application of FMINSEARCH in this problem. FMINSEARCH isn't guaranteed to converge when there is more than 1 unknown (even though it does pretty well empirically for small numbers of unknowns). Conversely FMINSPLEAS will reduce the problem to 1 unknown variable and apply FMINSEARCH to that, guaranteeing convergence.

[Star Strider]

The fminsearch function worked well here when I tested it.

Simple, two lines of code, and problem solved.

[Matt J]

The performance will be more sensitive to your initial guess, however. The following test shows that you get more robust behavior from fminspleas, when initializing with 2x the true parameters, even with 1/10th the noise of your original test,

    t = linspace(0,1,10);                       % Create Bin Data
    H = 3.*exp(-t/5) + 7 + 0.001*rand(size(t));  % Create Histogram Data

    [tau,BA] = fminspleas({@(tau,t)exp(-t./tau),1},10,t,H);

    disp 'FMISNPLEAS'
    fprintf(1, '\n\tB   = %.5f\n\ttau = %.5f\n\tA   = %.5f\n\n',[BA(1),tau, BA(2)])

    N = @(b,t) b(1).*exp(-t./b(2)) + b(3);      % Objective Function
    SSE = @(b) sum((H - N(b,t)).^2);            % Cost Function (Sum Squared Error)
    [b, fv] = fminsearch(SSE, [6 10 14]);
    disp 'FMINSEARCH'
    fprintf(1, '\n\tB   = %.5f\n\ttau = %.5f\n\tA   = %.5f\n\n',b)

Results:

    FMISNPLEAS

      B   = 2.98051
      tau = 4.96771
      A   = 7.01970

    FMINSEARCH

      B   = 6.81130
      tau = 12.09818
      A   = 3.18261

[Star Strider]

I still favor the simpler approach with fminsearch. I don’t understand why this has become such a contentious issue. There could be many appropriate parameter sets providing an equivalent fit.

Let Katelyn decide what most closely meets her needs. Neither of us know.

[Matt J]

I don’t understand why this has become such a contentious issue.

Has it? I think we want to be able to have robust discussions here about the strengths and pitfalls of the different methods proposed. Part of the point of the forum is to dig into things and pull out subtleties.

I'm not sure why you view fminsearch as simpler. They both require comparable lines of code.