Curve Fitting Problem

[Fredrik Lingvall]

A. Kalten wrote:
> Hello,
>
> I need some advice on solving a curve fitting problem.
>
> Ordinarily, a curve fitting algorithm, such as polyfit
> or wpolyfit in octave, will determine the coefficients of a
> polynomial that best fits the empirical data.  If the model
> is a second degree polynomial, a x^2 + b x + c, the algorithm
> will return the best coefficients a, b, and c.
>
> The trouble I am having is that I need to fit a second
> degree polynomial where some of the coefficients are
> either already known or include another known parameter
> as a factor.  For example, I need to find the best value
> for k in this equation:
>
> k^2/4 x + k p x + p^2
>
> where p is a known quantity.
>
> Although there are three coefficients, there is really
> only one unknown parameter, that is k.  Using polyfit from octave
> will return a, b, and c where a = k^2/4, b = k p, and c = p^2.
>
> Solving the two equations for k gives values that
> are significantly different (they differ in the second
> decimal place).
>
> Is there an algorithm available that allows some of the
> parameters of the model to be predetermined or to include
> some predetermined factor?
>
> I suspect that the only solution would be to write a
> custom least squares program.
>
> AK
>
>   
AK,

Do you have any prior information of the parameter k? Do you, for 
example, know if k is positive, or do you have any idea how much it may 
vary (variance of k),  mean value, if k is bound to some interval, etc?

An idea is then to formulate the problem as a parameter estimation problem, using the model

y = k^2/4 x^2 + k p x + p^2 + e

where e "describes" the misfit of the model. Say, for example, that you know that k is in the interval [a,b] and by assuming Gassian errors, e, 
you can find the most likely k (i.e, the maximum a posteriori estimate) by finding the k which maximizes


\hat{k} = arg max { exp (-1/ (2 sigma_e^2)  ( Y-k^2/4 X.^2 + k p X + 
p^2)^T ( Y-k^2/4 X.^2 + k p X + p^2) ) }
                       k

for k in [a,b] where Y is the data vector Y =[y(1) y(2) ... y(N)]^T.  If 
you instead have knowledge of the mean and variance of k then you can 
assign a Gaussian prior for k which will result in a different 
(unbounded) optimization problem.

/Fredrik



 






-- 
Fredrik Lingvall, PhD                                                              
E-mail: fl at ifi.uio.no, fredrik.lingvall at gmail.com
Web:    http://folk.uio.no/fl/