curvefitting comparison in Kaleidagraph and Matlab

Hello

At the moment I'm curve fitting longitudinal growth data into a non-parametric model containing 8 parameters.

At first I used kaleidagraph. A curve fitting program using the levenberg-marquardt algorithm. (as to vary between the inverse Hessian method and the steepest descent method). Starting values are putting in and minimization proceeds iteratively.

Trying to reproduce this using matlab i've tried following code opts = optimoptions(@lsqcurvefit, 'Algorithm', 'levenberg-marquardt','DerivativeCheck','on','FinDiffType', 'central','scaleproblem','Jacobian'); [personalParams,personalRes,personalResidual] = lsqcurvefit(heightModel,initialValues,personalData(:,1),personalData(:,2),[],[],opts);

In the opts command only the algoritm has stayed the same, for the rest I've tried to make a lot of the combinations I could make.

But if i compare the function parameters obtained from matlab or kaleidagraph, results vary dramastically. Curve fitting is however comparable.

How can I make sure that both curve fitting procedures will give an exact match. I've looked into some data as how the curve fitting occurs but maybe someone can give me pointers?

I'm a biologist so my background in mathematics is not this big.

Thanks in advance

Niels

Answer 1

[Alan Weiss]

I doubt that you can get two different implementations, the lsqcurvefit one and the Kaleidagraph one, to take the same iterations. Sorry.

You can examine the lsqcurvefit code to see exactly what it is doing:

 edit lsqcurvefit

When you come to a subroutine, such as lsqncommon, click the subroutine name and hit control-D to see that code, too. But I doubt that Kaleidagraph makes its code available to its users, so this might not help you in any case.

I am not really sure what sort of assurance you are looking for. Perhaps if you tell us what you really want (provably optimal parameters? fit guaranteed to be within 1% of optimum? ISO certification of underlying algorithm?) then we might be able to suggest something.

Alan Weiss

MATLAB mathematical toolbox documentation

[Niels]

thanks for answering

What I would want is that my functional parameters more or less have a normal gaussian distribution (by saphiro wilks) if taken enough samples.

What it seems like is that in some cases the distribution of 1 parameter just centers extremely around one specific value. Kinda like it is more or less treated as a fixed value (again I'm not a mathematic so I could be totally wrong here). This is in most of the cases the latest parameter.

Is there a possibility that lsqcurvefit has problems with handeling up to 8 initial values?

Greetings

[Alan Weiss]

lsqcurvefit treats all parameters symmetrically. So I do not think that it is the case that the "latest parameter" has any special significance.

You might want to check this section of the documentation for suggestions on how to improve your fit. In particular, the starting point is very important for many nonlinear fitting problems, and it is possible that your problem would benefit from trying a variety of starting points. If you have Global Optimization Toolbox, the MultiStart solver can try these points automatically. Be sure to give sensible bounds on all components.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Niels]

Thanks

I've found the code written in C I think used by the other program to perform the fitting. This is written below, for me this is a bit to complex but this is how i want to perform my fittings:

I'm hoping that you'll say

• oh thats this code in matlab that does the trick for you :-)
• that i can just copy paste this code and use it

thanks in advance

void mrqmin(float x[], float y[], float sig[], int ndata, float a[], int ia[],
int ma, float **covar, float **alpha, float *chisq,
void (*funcs)(float, float [], float *, float [], int), float *alamda)

{
void covsrt(float **covar, int ma, int ia[], int mfit);
void gaussj(float **a, int n, float **b, int m);
void mrqcof(float x[], float y[], float sig[], int ndata, float a[],
int ia[], int ma, float **alpha, float beta[], float *chisq,
void (*funcs)(float, float [], float *, float [], int));
int j,k,l;
static int mfit;
static float ochisq,*atry,*beta,*da,**oneda;
if (*alamda < 0.0) { Initialization.
atry=vector(1,ma);
beta=vector(1,ma);
da=vector(1,ma);
for (mfit=0,j=1;j<=ma;j++)
if (ia[j]) mfit++;
oneda=matrix(1,mfit,1,1);
*alamda=0.001;
mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,chisq,funcs);
ochisq=(*chisq);
for (j=1;j<=ma;j++) atry[j]=a[j];
}
for (j=1;j<=mfit;j++) { Alter linearized fitting matrix, by augmenting difor (k=1;k<=mfit;k++) covar[j][k]=alpha[j][k]; agonal elements.
covar[j][j]=alpha[j][j]*(1.0+(*alamda));
oneda[j][1]=beta[j];
}
gaussj(covar,mfit,oneda,1); Matrixsolution.
for (j=1;j<=mfit;j++) da[j]=oneda[j][1];
if (*alamda == 0.0) { Once converged, evaluate covariance matrix.
covsrt(covar,ma,ia,mfit);
covsrt(alpha,ma,ia,mfit); Spread out alpha to its full size too.
free_matrix(oneda,1,mfit,1,1);
free_vector(da,1,ma);
free_vector(beta,1,ma);
free_vector(atry,1,ma);
return;
}
for (j=0,l=1;l<=ma;l++) Did the trial succeed?
if (ia[l]) atry[l]=a[l]+da[++j];
mrqcof(x,y,sig,ndata,atry,ia,ma,covar,da,chisq,funcs);
if (*chisq < ochisq) { Success, accept the new solution.
*alamda *= 0.1;
ochisq=(*chisq);
for (j=1;j<=mfit;j++) {
for (k=1;k<=mfit;k++) alpha[j][k]=covar[j][k];
beta[j]=da[j];
}
for (l=1;l<=ma;l++) a[l]=atry[l];
} else { Failure, increase alamda and return.
*alamda *= 10.0;
*chisq=ochisq;
}
}

{
int i,j,k,l,m,mfit=0;
float ymod,wt,sig2i,dy,*dyda;
dyda=vector(1,ma);
for (j=1;j<=ma;j++)
if (ia[j]) mfit++;
for (j=1;j<=mfit;j++) { Initialize (symmetric) alpha, beta.
for (k=1;k<=j;k++) alpha[j][k]=0.0;
beta[j]=0.0;
}
*chisq=0.0;
for (i=1;i<=ndata;i++) { Summation loop over all data.
(*funcs)(x[i],a,&ymod,dyda,ma);
sig2i=1.0/(sig[i]*sig[i]);
dy=y[i]-ymod;
for (j=0,l=1;l<=ma;l++) {
if (ia[l]) {
wt=dyda[l]*sig2i;
for (j++,k=0,m=1;m<=l;m++)
if (ia[m]) alpha[j][++k] += wt*dyda[m];
beta[j] += dy*wt;
}
}
*chisq += dy*dy*sig2i; And find χ
2
.
}
for (j=2;j<=mfit;j++) Fill in the symmetric side.
for (k=1;k<j;k++) alpha[k][j]=alpha[j][k];
free_vector(dyda,1,ma);
}