| |
 |
|
|
Science Forum Index » Math - Numerical Analysis Forum » Biexponential curve-fitting
Page 1 of 1
|
| Author |
Message |
| spasmous |
 Posted: Thu Mar 08, 2007 5:55 pm |
|
|
|
Guest
|
I've been thinking about curve-fitting a biexponential model to some
data:
s = A.exp(-B.t) + C.exp(-D.t)
One problem seems that the model has two identical solutions (by
swapping A & C and B & D). When A->C and B->D the Jacobian becomes
singular. One way I thought of to avoid this is to use
s = A.exp(-B.t) + C.exp(-(B+D).t)
as the fitting function.It is identical, but the Jacobian cannot be
singular except when D->0, which can be excluded easily by bound
constraints. Simulation results show some slight improvement in
reliability using MATLAB's lsqcurvefit function. Just wondering if
anyone has any objection to this approach - or any better ideas? |
|
|
| Back to top |
|
| Peter Spellucci |
| Posted: Fri Mar 09, 2007 5:55 am |
|
|
|
Guest
|
In article <1173390944.187218.69330@30g2000cwc.googlegroups.com>,
"spasmous" <spasmous@gmail.com> writes:
Quote: I've been thinking about curve-fitting a biexponential model to some
data:
s = A.exp(-B.t) + C.exp(-D.t)
One problem seems that the model has two identical solutions (by
swapping A & C and B & D). When A->C and B->D the Jacobian becomes
singular. One way I thought of to avoid this is to use
s = A.exp(-B.t) + C.exp(-(B+D).t)
as the fitting function.It is identical, but the Jacobian cannot be
singular except when D->0, which can be excluded easily by bound
constraints. Simulation results show some slight improvement in
reliability using MATLAB's lsqcurvefit function. Just wondering if
anyone has any objection to this approach - or any better ideas?
this is quite o.k. you might do even better by using the
structure of the problem, expressing A and C in terms of B and D
from the optimality conditions which reduces the dimension of the
problem to 2.
hth
peter |
|
|
| Back to top |
|
| Fred Krogh |
| Posted: Fri Mar 09, 2007 5:45 pm |
|
|
|
Guest
|
spasmous wrote:
Quote: On Mar 9, 7:55 am, spellu...@fb04373.mathematik.tu-darmstadt.de (Peter
Spellucci) wrote:
In article <1173390944.187218.69...@30g2000cwc.googlegroups.com>, "spasmous" <spasm...@gmail.com> writes:
I've been thinking about curve-fitting a biexponential model to some
data:
s = A.exp(-B.t) + C.exp(-D.t)
One problem seems that the model has two identical solutions (by
swapping A & C and B & D). When A->C and B->D the Jacobian becomes
singular. One way I thought of to avoid this is to use
s = A.exp(-B.t) + C.exp(-(B+D).t)
as the fitting function.It is identical, but the Jacobian cannot be
singular except when D->0, which can be excluded easily by bound
constraints. Simulation results show some slight improvement in
reliability using MATLAB's lsqcurvefit function. Just wondering if
anyone has any objection to this approach - or any better ideas?
this is quite o.k. you might do even better by using the
structure of the problem, expressing A and C in terms of B and D
from the optimality conditions which reduces the dimension of the
problem to 2.
hth
peter
Hm, I'm not sure I understand what you mean. How do I express A and C
in terms of B and D?
Use google to search on "variable projection algorithm".
Fred |
|
|
| Back to top |
|
| Guest |
| Posted: Fri Mar 09, 2007 6:26 pm |
|
|
|
|
On Mar 8, 2:55 pm, "spasmous" <spasm...@gmail.com> wrote:
Quote: I've been thinking about curve-fitting a biexponential model to some
data:
s = A.exp(-B.t) + C.exp(-D.t)
One problem seems that the model has two identical solutions (by
swapping A & C and B & D). When A->C and B->D the Jacobian becomes
singular. One way I thought of to avoid this is to use
s = A.exp(-B.t) + C.exp(-(B+D).t)
as the fitting function.It is identical, but the Jacobian cannot be
singular except when D->0, which can be excluded easily by bound
constraints. Simulation results show some slight improvement in
reliability using MATLAB's lsqcurvefit function. Just wondering if
anyone has any objection to this approach - or any better ideas?
Is there significant noise in your data? If not, try Prony's method
(google for it) |
|
|
| Back to top |
|
| Alex. Lupas |
| Posted: Sat Mar 10, 2007 12:51 am |
|
|
|
Guest
|
On Mar 8, 11:55 pm, "spasmous" <spasm...@gmail.com> wrote:
Quote: I've been thinking about curve-fitting a biexponential model to some
data:
s = A.exp(-B.t) + C.exp(-D.t)
One problem seems that the model has two identical solutions (by
swapping A & C and B & D). When A->C and B->D the Jacobian becomes
singular. One way I thought of to avoid this is to use
s = A.exp(-B.t) + C.exp(-(B+D).t)
as the fitting function.It is identical, but the Jacobian cannot be
singular except when D->0, which can be excluded easily by bound
constraints. Simulation results show some slight improvement in
reliability using MATLAB's lsqcurvefit function. Just wondering if
anyone has any objection to this approach - or any better ideas?
======================
Next, by biexponential model of some data I mean
a linear combination of {1,exp(-Bt),exp(-Ct)} where
B,C are non-zero, different (given) numbers.
Then s=s(t)= N + M*exp(-Bt)+P*exp(-Ct)
and the 1question is to find N,M,P.
=========
In other words, consider following table:
x_1 , x_2,..., x_n
-------------------------
f(x_1),f(x_2),...,f(x_n)
{a=<x_j=<b for j in {1,2,...,n}}
and let
A(g):=(1/n)*SUM_{1=<k=<n}g(x_k)
(f,g):=A(fg)=(1/n)*SUM_{1=<k=<n}f(x_k)g(x_k)
||f||:=sqrt{(f,f)}
where f,g are real functions defined at x_j.
Further {1,f_1,f_2} are real functions defined on [a,b]
such thta
i) f_1 is sstrict monotonic on [a,b],
ii) for any system {t_1,t_2,t_3} of distinct points from [a,b],
the determinant
| 1 f_1(t_1) f_2(t_1) |
| 1 f_1(t_2) f_2(t_2) |
| 1 f_1(t_3) f_2(t_3) |
is non-zero.
Denote by T_2 the set of all linear combinations of
{1,f_1,f_2}. Of interest is to find in T_2 (if existst)
an element
(1) V*(t)= M + N*f_1(t) +P*f_2(t)
such that
(2) ||f-V*||=min_{V in T_2}||f-V||.
In your case
f_1(t):=e^{-Bt}
f_2(t):=e^{-Ct}
where B,C are non=zero and B=/=C.
Solution:
=========
1) find A(f_1) and then
X:= f_1^{-1}(A(f_1)) .
2) Construct
P_1(t):=f_1(t)-A(f_1)
P_2(t):= C0+C1*f_1(t)+C2*f_2(t)
where
C0:= A(f_1)A(f_1*f_2)-A(f_2)A(f_1^2)
C1:= A(f_1)A(f_2)-A(f_1f_2)
C2:= A(f_1^2)-A^2(f_1).
3) For j in {1,2}, find
w_j:=1/(P_j,P_j) = 1/A(P_j^2).
4) The unique element satisfying (1)-(2), is
V*(t)=A(f)+w_1*A(fP_1)P_1(t)+ w_2*A(fP_2)P_2(t). |
|
|
| Back to top |
|
| Thomas Mautsch |
| Posted: Sun Mar 11, 2007 8:52 pm |
|
|
|
Guest
|
In news:<1173390944.187218.69330@30g2000cwc.googlegroups.com>
schrieb spasmous <spasmous@gmail.com>:
Quote: I've been thinking about curve-fitting a biexponential model to some
data:
s = A.exp(-B.t) + C.exp(-D.t)
One problem seems that the model has two identical solutions (by
swapping A & C and B & D). When A->C and B->D the Jacobian becomes
singular. One way I thought of to avoid this is to use
s = A.exp(-B.t) + C.exp(-(B+d).t)
as the fitting function. It is identical, but the Jacobian cannot be
singular except when d->0, which can be excluded easily by bound
constraints. Simulation results show some slight improvement in
reliability using MATLAB's lsqcurvefit function. Just wondering if
anyone has any objection to this approach - or any better ideas?
I am no expert, but I believe the following:
What your second approach does
- restricting the variable d to positive values -
is to break the symmetry (A,B) <--> (C,D),
so that instead of (every pair of)
two equivalent global minima you get just one minimum (each).
The gain is obvious: While in your original formulation,
it could have been that the optimization algorithm
might oscillate between two equivalent minima - ending nowhere,
in your modified version these two equivalent minima are replaced
by one minimum, and the described kind of oscillation is prevented.
However, the other thing you write is wrong, I think:
Yes, the Jacobian of s = s(A,B,C,D) = A.exp(-B.t) + C.exp(-D.t)
is singular at B = D (your other condition, A = C, is not
needed for this to happen); BUT you cannot exclude
the effect of the singularity by setting D := B + d with d > 0.
The issue is that for some data
you must end up in singular points.
-
Just consider what happens when your data points have the form
(t_i, s_i = t_i * exp(b*t_i) ), i=1..N,
for some constant b.
In this case, the "best" approximation
s = A.exp(-B.t) + C.exp(-D.t)
will have
D -> B
and
|C-A| -> infinity.
The reason for this is:
What you have is the 4-parameter family
(A,B,C,D) |--> ( A.exp(-B.t_i) + C.exp(-D.t_i) )_{i=1..N}
in N-dimensional Euclidean space R^N. You want to pick (A,B,C,D)
so that the image point of (A,B,C,D) gets as near as possible
to a given point ( s_i )_{i=1..N}.
BUT, the image of your 4-parameter family is NOT CLOSED in R^N;
so it might happen that the best approxing point does not lie
in the image of the 4-parameter family, but only in its closure. -
There will be no parametrization in the form of your model for this
point.
So what you have to do is to add further models
which parametrize the degenerate cases.
I'd say there are three possibilities
at what points (A,B,C,D) the degeneration problem might occur:
Points where B and D get near to each other,
points where max(B,D) tends to +infinity,
and the corresponding factor A, resp. C becomes very large
to counter-affect the smallness of the exponential,
and points where min(B,D) tends to -infinity,
and the corresponding factor A, resp. C becomes very small.
What you can do is:
When the first case (B=D) occurs,
use the additional model
(A,a,B) |--> s = (A + a*t) * exp(-B*t).
It might be that approximation by this model is better
than by the bi-exponential model. However,
in this model it might also happen that B tends to + or -infinity.
When the second case, D=+infinity occurs,
try the additional model
(A,B) |--> s = A * exp(-B*t)
but instead of using the full data set
delete the data point (t_i,s_i) with smallest value t_i.
Similarly in the third case B = -infinity:
Try the additional mode
(C,D) |--> s = C * exp(-D*t)
on the given data points (t_i,s_i) with the point with largest t_i
deleted.
(Let us hope that I have not mixed up the last two cases... ;-( )
As I already said: I am not an expert; but I think this
should be the basic picture that you encounter when
doing least-square fits with exponential models.
Regards,
Thomas Mautsch |
|
|
| Back to top |
|
| Ray Koopman |
| Posted: Mon Mar 12, 2007 3:24 am |
|
|
|
Guest
|
spasmous wrote:
Quote: I've been thinking about curve-fitting a biexponential model to some
data:
s = A.exp(-B.t) + C.exp(-D.t)
One problem seems that the model has two identical solutions (by
swapping A & C and B & D). When A->C and B->D the Jacobian becomes
singular. One way I thought of to avoid this is to use
s = A.exp(-B.t) + C.exp(-(B+D).t)
as the fitting function.It is identical, but the Jacobian cannot be
singular except when D->0, which can be excluded easily by bound
constraints. Simulation results show some slight improvement in
reliability using MATLAB's lsqcurvefit function. Just wondering if
anyone has any objection to this approach - or any better ideas?
If you can plot z = f(x,y) in some fashion such as gray-levels,
and you know the rectangular region in which the optimal B and D lie,
then set up a (B,D) grid; for each point in that grid, plot R_{3,3}^2
from a QR decomposition of the n x 3 matrix in which row i is
[ exp(-B.t_i), exp(-D.t_i), s_i ]. Look for the overall minimum,
zooming in as needed. When you have identified the optimal B and D,
use ordinary least squares to find A and C. |
|
|
| Back to top |
|
| |
|
Page 1 of 1
All times are GMT - 5 Hours
The time now is Thu Jan 08, 2009 3:43 am
|
|