Best way to fit functionals with error bars

Hi,

I have the following expression for a rate coefficient . For the purposes of this question, ρ is a smooth continous increasing function of E. I have values of k (with an error bar associated with them) for different values of E. I need to find the best fitting that will fit the values i have (keeping into account the error bars). ρ does not have an analytical form, i need to find numerically compute it. So basically i have an array of ρ pre computed.

What is the best way to do this?

Thanks!

[Matt J]

keeping into account the error bars

Meaning what? You want to constrain the final fit to fall within the corridor between the error bars? With only 2 degrees of freedom, that might be difficult to enforce rigidly.

[Arun]

I am not sure how to even begin fitting. Lets say I have rho = [array of size 1*1e6] where every data point is rho(e) for some e = [array of size 1*16].

How do i do fit k(e) which is this ratio of these rhos at different energies. It is not an analytical function, so i can not use fit().

If i figure this out, I could then start thinking about accounting for the spread in k(e)

Answer 1

[Star Strider]

That depends on what the error bars indicate, and the information you have about them.

A common way of weighting a regression (linear or nonlinear) is to use inverse variance weighting.. If you can extract the vairances from the error bars, then you can use this method. (This may involve inverting the t-distribution or normal distribution, or the distribution that describes them. A common criterion is to define them as the 95% confidence interval (2.5% to 97.5%), so use that as the default.)

Otherwise, an alternate (and to the best of my knowledge, unproven) method would be to use a function of the inverses of the error bar lengths.

[Arun]

I usually do inverse variance weighing if I have the analytical form of the function to fit the data with. But here I dont, as \rho does not have an analytical form that I can define as a fittype to pass into the fit function. I am more interested to know how to do the fitting itself>

[Star Strider]

I am not certain what you are asking.

If ρ does not have an analytic form although is something you can compute (I have no idea what its independnet variable would be), and knowing nothing more about it, one option sould be to interpolate it using the interp1 function, since I assume that it is a function of one (independent) variable. If it has more than one variable (for example one or more parameters), what you would do would depend on the ρ function itself. I have no idea what that is. You might still be able to interpolate it.

As a general rule, the least square estimation would go something like this --

where is the known value, is the regression estimate and a function of the independent variable and parameters, and is the weighting value at that point.

Everything now depends on what ρ is, whether it has any parameters you might want to fit, and what you want to do with it.

I generally use the Statistics and Machine Learning or Optimization or Global Optimization Toolboxes for these sorts of problems.

.

[Arun]

ρ is a function of E (E is the only independent variable).

So if i want to find k at E, then k is α times the ratio of ρ at E - Ea and ρ at E. α and Ea are the two parameters I have to find by fitting the functional form i have in the question to the values of k i have at some energies, E (These values themselves have uncertainities, how to handle them is my second question, but i understand now)

Is there a function in these toolboxes that does the minimization of the least square estimate that you propose?

[Matt J]

Is there a function in these toolboxes that does the minimization of the least square estimate that you propose?

It depend on what constraints (if any) you want to apply. If you only intend to apply simple bounds on α and , then fit() with a custom fittype would work. However, if you need more complicated constraints, you will need something like lsqcurvefit in the Optimization Toolbox.

Either way, the representation of is not the issue. You would just use interpolation as mentioned but @Star Strider

rho=@(E) interp1(E_samp,rho_samp,E,'linear','extrap') %continuous rho()

[Star Strider]

'I need to find the best fitting and α that will fit the values i have (keeping into account the error bars).'

This defines and α as parameters to fit.

The sort of function I would use for this would be something like:

k_fcn = @(b,E) b(1).*rhofcn(E-b(2))./rhofcn(E);

with 'rhofcn' defined as @Matt J defined it above. (I call it 'rhofcn' for clarity, avoiding problems since 'rho' could be used for something else.) Here, 'b(1)' is α, and 'b(2)' is . The independent variable is E.

You would then use 'k_fcn' as the objective function with lsqcurvefit or fitnlm. Use the weighting vector derived from the error bar values with it.

.

[Arun]

This approach worked. Thank you!

[Star Strider]

As always, my pleasure!