Prev: Open a BMP image
Next: Changing Default Stream within Simulink for generating random numbers
From: Pat on 28 Apr 2010 18:15
I'm reading in data that is 10,000 points long. I then select a certain section of the data...say 8500 points. I fit an exponential to this data using fminsearch which works great. Finally, I need to calulate the error of the fit which brings us to the Jacobian, the Covariance matrix and the standard error.

COV3=Standard_Error^2*inv(jac'*jac)

This gives an "Out of Memory" response and I was wondering if there was a decent way to compute the inverse without getting the "Out of memory" error.
Fit Function: F(time)=A+B.*exp(-c.*time)
jac will always be less than 10,000 points.

Thanks in advance.
From: Bruno Luong on 29 Apr 2010 01:02
"Pat " <heisenberg12501-shopping(a)yahoo.com> wrote in message <hrac1p$akg$1(a)fred.mathworks.com>...
> I'm reading in data that is 10,000 points long. I then select a certain section of the data...say 8500 points. I fit an exponential to this data using fminsearch which works great. Finally, I need to calulate the error of the fit which brings us to the Jacobian, the Covariance matrix and the standard error.
>
> COV3=Standard_Error^2*inv(jac'*jac)
>
> This gives an "Out of Memory" response and I was wondering if there was a decent way to compute the inverse without getting the "Out of memory" error.
> Fit Function: F(time)=A+B.*exp(-c.*time)
> jac will always be less than 10,000 points.
>

It does not make sense: I have a doubt whereas the dimension of your jacobian is even correct. The dimension of jac should be 8500 x 3. The memory requirement for inv(jac'*jac) is small.

In any case please attach the code together with verbose description.

Bruno
From: Pat on 30 Apr 2010 17:12
Bruno my friend! You are absolutely 100% correct!!
I messed up the Jacobian.

I think I've got it right now, but lets check:
In reality function is the following:
f(t)=A+Bexp(-t/Tau) interested in Tau (lifetime).
In program fitting function is the following:
f(t)=A+B.*exp(-lambda.*time4)
A=Estimates(1)
B=Estimates(2)
lambda=Estimates(3)

-----------------------------------------------------------------------------------------
Estimates=fminsearch(@myfit2,Starting,options,time4,data4);
tau_1=1/abs(Estimates(3)); %This is lifetime
y_exp=Estimates(1)+Estimates(2).*exp(-Estimates(3).*time4);
corr_exp=corrcoef(data4,y_exp);
exp_r_square=abs(corr_exp(1,2));

%Calculating Error of Fit
%Jacobian
jac=[ones(length(time4),1),exp(-Estimates(3).*time4'),(-1*Estimates(3)^2)*(Estimates(2).*exp(-Estimates(3).*time4')).*-time4'];
%Compute Residuals
res=y_exp-data4;
%Compute standard error dividing by the minus 3 due to 3 Estimate parameters
ser=sqrt(sum(res.^2)/(length(time4)-3));
COV3=ser^2*inv(jac'*jac);
StdDev3=sqrt(diag(COV3));%Standard Deviation of Fitting Parameters
__________________________________________________________

StdDev3(3) is the number I'm interested in and I've had some success calculating it.
However I've noted that as my R-squared fit value goes down, I'm getting
Warning Matrix is close to Singular Results May be Inaccurate

I agree I think they are inaccurate.

Thanks in Advance!
From: us on 30 Apr 2010 17:31
"Pat " <heisenberg12501-shopping(a)yahoo.com> wrote in message <hrfh34$56b$1(a)fred.mathworks.com>...
> Bruno my friend! You are absolutely 100% correct!!
> I messed up the Jacobian.
>
> I think I've got it right now, but lets check:
> In reality function is the following:
> f(t)=A+Bexp(-t/Tau) interested in Tau (lifetime).

just a note:
this problem can be dealt with greatly with NLINFIT available in the stats tbx, which you ...hopefully... own...
it returns all you're looking for and more...

us
From: Pat on 3 May 2010 12:25
Hopefully this doesn't double post...had some issues earlier.
Created function:
function Fitted_Curve=myfit5(params,Input)
A=params(1);
B=params(2);
lambda=params(3);

Fitted_Curve=A+(B.*exp((-1/lambda).*Input));

Used in Main Program:
[Estimates,R,J,COVB,MSE]=nlinfit(time4,data4,@myfit5,Starting,options);
tau_1=Estimates(3);
y_exp=Estimates(1)+Estimates(2).*exp((-1/Estimates(3)).*time4);

In my opinion the fitting using the fminsearch is better. I get higher r-squared values when using that method. I do like the fact that nlinfit gives me what I'm looking for quickly and in one command, but it doesn't seem to work quite as well as fminsearch. Also, I still get warnings such as:
The Jacobian at the solution is ill-conditioned, and some
model parameters may not be estimated well (they are not identifiable).

The matrix errors I was receiving before were only on ill-behaving data sets.
 | 
Pages: 1
Prev: Open a BMP image
Next: Changing Default Stream within Simulink for generating random numbers