curve fitting with numeric array instead of analytic form for function of x

Hello,

Sorry if this is a silly question, but I've been searching through the help documentation and I've been unable to find an answer. Does MATLAB have the ability to do linear least squares fitting of the type y = a*f1(x) + b*f2(x) + ..., where fi(x) are numerically defined functions of the independent variable? I have f1, f2, f3, and f4 defined in vectors as the normalized (to maximum value) results of source signal measurements. I have another set of signal measurements that I want to 'decompose' into a linear sum of the four sources, with amplitudes determined by least squares.

I've found that I can define functions of 'x' and use the matlab function 'fit' with no problem, but I'm unable to write down an analytic functional form for this type of model, so I don't think the curve fitting toolbox will be of much use. Does anyone have any suggestions?

Thanks a lot,

Chris

[Star Strider]

Do your functions f1(x) ... fi(x) produce vectors equal to length(y) (assuming y is a vector and not a matrix), or are they nonlinear functions that have to be evaluated during the parameter estimation process?

[Chris]

Hello,

Yes, the fi(x) are vectors of the same length y, the data I would like to decompose.

Answer 1

[Star Strider]

‘Yes, the fi(x) are vectors of the same length y, the data I would like to decompose.’

In that situation, this seems to me a multiple linear regression problem. If your matrix of f(*) vectors is F, your vector of parameters a, b, ... is P, then your system becomes (with the appropriate dimensions of F):

y = F*P

and the solution:

P = F\y;

A;though you may need to augment F with a constant vector such as:

F1 = [ones(size(F,1)  F];

and then:

P = F1\y;

A more formal implementation of this — with statistics — is regress in the Statistics Toolbox.

[Chris]

This would work, except that regress does not allow for use of measurement error as weights for fitting. Thanks though! It is looking like I'll be writing a function myself...

[Chris]

I'm inverting a matrix to find the best fit coefficients in a regular least squares sense, with consideration of measurement error in the dependent variable only... with the matrix defined as described in "The Art of Scientific Computing" data modeling chapter. I don't have the book with me right now, so sorry for not citing the specific page.

[Star Strider]

I didn't know you wanted to do weighted regression. I must have missed that.

It took a bit to find the online reference, but here's a demo how you can do Weighted Nonlinear Regression. I've used this technique in the past for inverse-variance weighting in both linear and nonlinear regression. This demo uses nlinfit but you can easily give it a linear system in your model. It will then give you access to nlparci and nlpredci. I've done multivariable regression with nlinfit in the past as well. It may require some experimenting, but barring unusual circumstances, it will work.

I'll keep this thread up for a while so I can find it easily.

[Chris]

The demo you referenced is what brought me to ask this question. In it they define:

modelFun = @(b,x) b(1).*(1-exp(-b(2).*x));

where the model function is written as an analytic function of x. In my case, Fi(x(i)) = Fi_i where Fi_i is the value of Fi at x_i. I have all the Fi tabulated, one for each of the x_i, but can not write down an explicit functional dependence on x.

I wonder if I can get around writing the F's as a function of x by using:

modelFun = @(a,b,c,d,F1,F2,F3,F4) a*F1 + b*F2 + c*F3 + d*F4;

This was my main problem with using 'fit' and why I kept complaining about needing an analytic function of x for the Fi's... I'm going to give it a shot and see how it works.

Thanks a lot for your time and responses!

[Chris]

Well, my idea did not work...

[Star Strider]

I'm feeling a bit chagrined right now that I didn't think of lscov earlier. See if it will do what you want. Its documentation says it can.

If asked, it will return the standard errors of the parameter vector, so you can calculate the parameter confidence intervals from them.

It will probably also let you use nlparci because it can return a covariance matrix for the parameters if asked. You will have to calculate your own vector of residuals, but that shouldn't be a problem.

[Chris]

Thanks a lot this will work... but I was actually able write up something pretty quickly that utilizes Singular Value Decomposition. Works good.

[Star Strider]

I hadn't thought about using SVD for WLS regression. I went online and found lecture slides from the University of Southampton and a paper from the Technical University of Denmark that discusses it. I doubt I'd have even looked had you not mentioned it, so thanks!

And thank you for accepting my answer.

[Chris]

Basically, I followed the description given in section 15.4 of Art of Scientific Computing and used MATLAB's svd function. I have to do some more research to bound the fit parameters tho...Here's what I have so far:

For a model of form a*F1(x) + b*F2(x) + ...

modarray is the tabulated response values for which the sum of coeff(i)*column(i) makes up the best fit model. Column(i) = Fi(x). indep is the independent variable, though it is not used. datatofit is a column vector of the data to fit. Assumes no error in independent variable.

function [ coeffs redchisq ] = nvar_lls( modarray, indep, datatofit, err )

%build the design matrix A and build vector b

%define vector b of length rows (number of data points)

[rows columns] = size(modarray);

b = zeros(rows,1);

A = zeros(size(modarray));

X = modarray; y = datatofit; sigma = err;

for i = 1:rows

    for j = 1:columns    

       A(i,j) = X(i,j)/sigma(i,1); 

    end

    b(i,1) = y(i,1)/sigma(i,1);

end

%denote the vector of best fit coefficients a

a = zeros(length(alphabet),1);

%get singular value decomposition of A

[U,S,V] = svd(A);

%fill the values of a

[srows scolumns] = size(S);

for i = 1:scolumns

   if S(i,i) > 0 

   a = a + (dot(U(:,i),b)/S(i,i))*V(:,i);

   else

   a = a + 0;

   end

end

coeffs = a;

%now to calculate the reduced chi-square value for the fit

numdof = rows - columns - 1; calcmodel = zeros(rows,1);

%build the model:

for i = 1:columns

    calcmodel = calcmodel + a(i)*X(:,i);

end

chisq = sum(((y-calcmodel)./sigma).^2);

redchisq = chisq/numdof;

[Star Strider]

Interesting! Thank you. The paper I found from the Technical University of Denmark is ‘Least Squares Adjustment: Linear and Nonlinear Weighted Regression Analysis’. I'm sending the link along in case you might find it of interest. It even has MATLAB code, but circa 2003. (I haven't had time to read it yet in any detail, but it appears to be a clear and cogent presentation.)

It seems I need to invest in a copy of The Art of Scientific Computing. I had no idea it had so much information.

[Chris]

Yes, get the book...it's really good! You can't accidentally flip to interesting pages in a pdf, though it is available online free. The only thing that, in my opinion, detracts from the online version is that you have to register with the website. You can usually find pdf versions of chapters and sections online at other web pages though. I'll check out the paper..thanks.

Chris

[Chris]

I have finished calculating the confidence limits on the fit parameters...let me know if you want the full script.

[Star Strider]

I would appreciate it, yes.

You might consider writing it up formally as a function, documenting it internally, and submitting it to the MATLAB File Exchange. (I believe the FEX has some suggestions on how best to do that. I haven't checked, since I haven't submitted any files in a while.) It seems to be needed. Even though lscov exists, not everyone may have the Statistics Toolbox, or want to go through everything you had to do to re-create your function.

Answer 2

[Teja Muppirala]

If I understand you correctly, then I think all you need is the "backslash" operator. Just like this:

x = (0:0.1:10)';

f1 = sin(x);
f2 = exp(x/10);
f3 = 3./(1+x);

y = 3*f1 - 4*f2 + 8*f3 + 1*randn(size(x));

plot(x,y);

estimated_coeffs = [f1 f2 f3]\y

hold on;
plot(x,[f1 f2 f3]*estimated_coeffs,'r');

[Chris]

Hello Teja,

No, that is my problem. I do not have an analytic expression for the fi(x). All I have is the tabulated values of fi(x). The fi(x) are emission spectra of different fluorophores and the y(x), the data I want to analyze, is assumed to be a linear sum of these spectra. I can not write down for example, that f1(x) = cos(x) because the functional form of the dependence on x is unknown.

Chris

[Tom Lane]

Teja's example calculated f1-f3 as expressions. I would think you could calculate them by fetching them out of the table, and still carry out the backslash operation that he recommends.

[Chris]

Yes, it is simple to get the 'best fit' coefficients through the backslash operator. This is basically how I am currently performing the fit. However, MATLAB already has functionality for evaluating goodness of fit (also simple) and confidence intervals for the fitted parameters(not as simple). I was hoping to be able to make use of that functionality, instead of finding the bounds of the hyperdimensional ellipse of the chi-square surface myself, in order to find confidence limits on my fitted parameters.