1
$\begingroup$

We have a measured data set with a relationship with the following formula:

$$Y = A\cdot \exp(B\cdot T) + C\cdot \exp(D\cdot T)$$

How can we approximate the values ​​of A, B, C and D? The method used should be fast and implementable with slow processors. The main purpose is programming in C language.

I approximated with the help of MATLAB software. The result was excellent, but I don't understand the working algorithm of MATLAB software.

clc
clear
%%
A = [
   43,   116,   389,   554,   674,   803,   861,   731,  1021,  1011,  1176,  1243,  1306,  1393,  1479, ...
 1459,  1635,  1694,  1745,  1795,  1803,  1850,  1780,  1869,  1806,  2003,  2140,  2001,  2056,  2154, ...
 2211,  2298,  2245,  2256,  2373,  2365,  2333,  2457,  2368,  2417,  2355,  2393,  2492,  2460,  2575, ...
 2492,  2492,  2570,  2533,  2637,  2601,  2576,  2596,  2677,  2536,  2640,  2632,  2660,  2767,  2784, ...
 2829,  2795,  2781,  2825,  2694,  2905,  2782,  2725,  2864,  2780,  2587,  2849,  2909,  2918,  2913, ...
 2878,  2885,  2837,  2909,  2836,  2899,  2972,  2835,  2897,  2999,  2851,  3019,  2970,  2941,  2932, ...
 2969,  3011,  2941,  2853,  2840,  2881,  3012,  2929,  2909,  3022,  2975,  2974,  2977,  2970,  3107, ...
 3043,  2996,  2985,  2859,  2940,  2982,  2941,  3051,  3066,  3049,  3107,  2997,  2976,  3089,  3023, ...
 3088,  3027,  3008,  2973,  2926,  2999,  2975,  3112,  3112,  3102,  3117,  2970,  3002,  3162,  3020, ...
 3075,  3020,  2989,  3043,  3017,  3066,  3216,  3089,  3250,  3205,  3080,  3137,  3151,  3078,  3095, ...
 3110,  3044,  3112,  3099,  3151,  3202,  3157,  3048,  3145,  3200,  3117,  3157,  3110,  3076,  3080, ...
 3011,  3233,  3194,  3165,  3214,  3225,  3053,  3169,  3090,  3204,  3189,  3119,  3165,  3084,  3017, ...
 3169,  3264,  3245,  3258,  3193,  3227,  3181,  3217,  3126,  3265,  3186,  3165,  3232,  3193,  3026, ...
 3306,  3325,  3240,  3289,  3240, ];

T = [
11316, 11658, 11986, 12328, 12670, 12998, 13340, 13668, 14010, 14338, 14680, 15022, 15350, 15692, 16020, ...
16362, 16690, 17032, 17374, 17702, 18044, 18372, 18714, 19042, 19384, 19726, 20054, 20396, 20724, 21066, ...
21394, 21736, 22078, 22406, 22748, 23076, 23418, 23746, 24088, 24430, 24758, 25100, 25428, 25770, 26098, ...
26440, 26782, 27110, 27452, 27780, 28122, 28450, 28792, 29134, 29462, 29804, 30132, 30474, 30802, 31144, ...
31486, 31814, 32156, 32484, 32826, 33154, 33496, 33838, 34166, 34508, 34836, 35178, 35506, 35848, 36190, ...
36518, 36860, 37188, 37530, 37858, 38200, 38542, 38870, 39212, 39540, 39882, 40210, 40552, 40894, 41222, ...
41564, 41892, 42234, 42562, 42904, 43246, 43574, 43916, 44244, 44586, 44914, 45256, 45598, 45926, 46268, ...
46596, 46938, 47266, 47608, 47950, 48278, 48620, 48948, 49290, 49618, 49960, 50302, 50630, 50972, 51300, ...
51642, 51970, 52312, 52654, 52982, 53324, 53652, 53994, 54322, 54664, 55006, 55334, 55676, 56004, 56346, ...
56674, 57016, 57358, 57686, 58028, 58356, 58698, 59026, 59368, 59710, 60038, 60380, 60708, 61050, 61378, ...
61720, 62062, 62390, 62732, 63060, 63402, 63730, 64072, 64414, 64742, 65084, 65412, 65744, 66086, 66414, ...
66756, 67084, 67426, 67754, 68096, 68438, 68766, 69108, 69436, 69778, 70106, 70448, 70790, 71118, 71460, ...
71788, 72130, 72458, 72800, 73142, 73470, 73812, 74140, 74482, 74810, 75152, 75494, 75822, 76164, 76492, ...
76834, 77162, 77504, 77846, 78174, ];
%%
BBZ = 3165;
T = T'/1000;
YN = BBZ-A';
%%
plot(T, YN)
[population2,gof] = fit(T, YN, 'exp2');
plot(population2,T,YN);

Matlab Fit

Another used method is described in the clip below. It works well for noise-free data, but gives wrong results with the smallest noise. https://www.youtube.com/watch?v=BL7k2s3tsqU

A solution to this question is given in link "Exponential regression with two terms and constraints". I put its implementation below and compared the result with the MATLAB method.

clc
clear
cla
%%

A = [
   43,   116,   389,   554,   674,   803,   861,   731,  1021,  1011,  1176,  1243,  1306,  1393,  1479, ...
 1459,  1635,  1694,  1745,  1795,  1803,  1850,  1780,  1869,  1806,  2003,  2140,  2001,  2056,  2154, ...
 2211,  2298,  2245,  2256,  2373,  2365,  2333,  2457,  2368,  2417,  2355,  2393,  2492,  2460,  2575, ...
 2492,  2492,  2570,  2533,  2637,  2601,  2576,  2596,  2677,  2536,  2640,  2632,  2660,  2767,  2784, ...
 2829,  2795,  2781,  2825,  2694,  2905,  2782,  2725,  2864,  2780,  2587,  2849,  2909,  2918,  2913, ...
 2878,  2885,  2837,  2909,  2836,  2899,  2972,  2835,  2897,  2999,  2851,  3019,  2970,  2941,  2932, ...
 2969,  3011,  2941,  2853,  2840,  2881,  3012,  2929,  2909,  3022,  2975,  2974,  2977,  2970,  3107, ...
 3043,  2996,  2985,  2859,  2940,  2982,  2941,  3051,  3066,  3049,  3107,  2997,  2976,  3089,  3023, ...
 3088,  3027,  3008,  2973,  2926,  2999,  2975,  3112,  3112,  3102,  3117,  2970,  3002,  3162,  3020, ...
 3075,  3020,  2989,  3043,  3017,  3066,  3216,  3089,  3250,  3205,  3080,  3137,  3151,  3078,  3095, ...
 3110,  3044,  3112,  3099,  3151,  3202,  3157,  3048,  3145,  3200,  3117,  3157,  3110,  3076,  3080, ...
 3011,  3233,  3194,  3165,  3214,  3225,  3053,  3169,  3090,  3204,  3189,  3119,  3165,  3084,  3017, ...
 3169,  3264,  3245,  3258,  3193,  3227,  3181,  3217,  3126,  3265,  3186,  3165,  3232,  3193,  3026, ...
 3306,  3325,  3240,  3289,  3240, ];

T = [
11316, 11658, 11986, 12328, 12670, 12998, 13340, 13668, 14010, 14338, 14680, 15022, 15350, 15692, 16020, ...
16362, 16690, 17032, 17374, 17702, 18044, 18372, 18714, 19042, 19384, 19726, 20054, 20396, 20724, 21066, ...
21394, 21736, 22078, 22406, 22748, 23076, 23418, 23746, 24088, 24430, 24758, 25100, 25428, 25770, 26098, ...
26440, 26782, 27110, 27452, 27780, 28122, 28450, 28792, 29134, 29462, 29804, 30132, 30474, 30802, 31144, ...
31486, 31814, 32156, 32484, 32826, 33154, 33496, 33838, 34166, 34508, 34836, 35178, 35506, 35848, 36190, ...
36518, 36860, 37188, 37530, 37858, 38200, 38542, 38870, 39212, 39540, 39882, 40210, 40552, 40894, 41222, ...
41564, 41892, 42234, 42562, 42904, 43246, 43574, 43916, 44244, 44586, 44914, 45256, 45598, 45926, 46268, ...
46596, 46938, 47266, 47608, 47950, 48278, 48620, 48948, 49290, 49618, 49960, 50302, 50630, 50972, 51300, ...
51642, 51970, 52312, 52654, 52982, 53324, 53652, 53994, 54322, 54664, 55006, 55334, 55676, 56004, 56346, ...
56674, 57016, 57358, 57686, 58028, 58356, 58698, 59026, 59368, 59710, 60038, 60380, 60708, 61050, 61378, ...
61720, 62062, 62390, 62732, 63060, 63402, 63730, 64072, 64414, 64742, 65084, 65412, 65744, 66086, 66414, ...
66756, 67084, 67426, 67754, 68096, 68438, 68766, 69108, 69436, 69778, 70106, 70448, 70790, 71118, 71460, ...
71788, 72130, 72458, 72800, 73142, 73470, 73812, 74140, 74482, 74810, 75152, 75494, 75822, 76164, 76492, ...
76834, 77162, 77504, 77846, 78174, ];
%%
%plot(T, YN)
%[population2,gof] = fit(T, YN, 'exp2');
%plot(population2,T,YN);
%%
x = T'/1000;
BBZ = 3165;
y = BBZ-A';


%figure(1)
hold off
%plot(x, y, 'o')
hold on
%%
SumXn     = 0;
SumYn     = 0;
SumX2n    = 0;
SumX3n    = 0;
SumX4n    = 0;
SumXnYn   = 0;
SumSn     = 0;
SumS2n    = 0;
SumSSn    = 0;
SumSnXn   = 0;
SumSnYn   = 0;
SumSS2n   = 0;
SumX2nYn  = 0;
SumSSnSn  = 0;
SumSSnYn  = 0;
SumSnX2n  = 0;
SumSSnXn  = 0;
SumSSnX2n = 0;
Sn        = 0;
SnLast    = 0;
SSn       = 0;
Items     = 0;
for I=1:length(x)
    Yn = y(I);
    Xn = x(I);
    X2n = Xn*Xn;
    X3n = Xn*X2n;
    X4n = Xn*X3n;
    if(I == 1)
        Sn      = 0;
        SnLast  = 0;
        SSn     = 0;
    else
        SnLast = Sn;
        Sn  = Sn  + 0.5*(y(I)+y(I-1))*(x(I)-x(I-1));
        SSn = SSn + 0.5*(Sn + SnLast)*(x(I)-x(I-1));
    end
    S2n  = Sn  * Sn;
    SS2n = SSn * SSn;

    SumSS2n   = SumSS2n   + SS2n;
    SumSSnSn  = SumSSnSn  + SSn * Sn;
    SumSSnX2n = SumSSnX2n + SSn * X2n;
    SumSSnXn  = SumSSnXn  + SSn * Xn;
    SumSSn    = SumSSn    + SSn;

    SumS2n    = SumS2n    + S2n;
    SumSnX2n  = SumSnX2n  + Sn * X2n;
    SumSnXn   = SumSnXn   + Sn * Xn;
    SumSn     = SumSn     + Sn;

    SumX4n    = SumX4n    + X4n;
    SumX3n    = SumX3n    + X3n;
    SumX2n    = SumX2n    + X2n;

    SumXn     = SumXn     + Xn;
    
    Items     = Items     + 1;

    SumSSnYn  = SumSSnYn  + SSn * Yn;
    SumSnYn   = SumSnYn   + Sn  * Yn;
    SumX2nYn  = SumX2nYn  + X2n * Yn;
    SumXnYn   = SumXnYn   + Xn  * Yn;
    SumYn     = SumYn     + Yn;
end

AA = ([ ...
    SumSS2n   SumSSnSn   SumSSnX2n   SumSSnXn   SumSSn; ...
    SumSSnSn  SumS2n     SumSnX2n    SumSnXn    SumSn ; ...
    SumSSnX2n SumSnX2n   SumX4n      SumX3n     SumX2n; ...
    SumSSnXn  SumSnXn    SumX3n      SumX2n     SumXn ; ...
    SumSSn    SumSn      SumX2n      SumXn      Items ; ...
]);
BB = [ ...
    SumSSnYn; ...
    SumSnYn ; ...
    SumX2nYn; ...
    SumXnYn ; ...
    SumYn   ; ...
];

RES = AA \ BB;

A = RES(1);
B = RES(2);

p = (B + sqrt(B^2 + 4*A)) / 2;
q = (B - sqrt(B^2 + 4*A)) / 2;

SumExPXn = 0;
SumExQXn = 0;
SumEx2PXn = 0;
SumEx2QXn = 0;
SumExPaQXn = 0;
SumYnExpPXn = 0;
SumYnExpQXn = 0;

for I=1:length(x)
    Xn = x(I);
    Yn = y(I);
    SumExPXn = SumExPXn + exp(p*Xn);
    SumExQXn = SumExQXn + exp(q*Xn);
    SumEx2PXn = SumEx2PXn + exp(2*p*Xn);
    SumEx2QXn = SumEx2QXn + exp(2*q*Xn);
    SumExPaQXn = SumExPaQXn + exp((p+q)*Xn);

    SumYnExpPXn = SumYnExpPXn + exp(p*Xn) * Yn;
    SumYnExpQXn = SumYnExpQXn + exp(q*Xn) * Yn;
end

AA = ([ ...
    Items    SumExPXn   SumExQXn  ; ...
    SumExPXn SumEx2PXn  SumExPaQXn; ...
    SumExQXn SumExPaQXn SumEx2QXn ; ...
]);
BB = [ ...
    SumYn      ; ...
    SumYnExpPXn; ...
    SumYnExpQXn; ...
];

RES = AA \ BB;

sA = RES(1);
sB = RES(2);
sC = RES(3);

YY = sA + sB*exp(p*x) + sC*exp(q*x);
%%
plot(x, YY, '-b');
disp('END');

ft = fittype('A*exp(B*x) + C*exp(D*x) + E');
[pop2,gof,output] = fit(x, y, ft, 'StartPoint', [1e3, -0.010, 1e4, -0.100, -500]);

plot(pop2, x, y);

Result:The output result of the code using the MATLAB method

$\endgroup$
4
  • $\begingroup$ $D$ does not appear in your function? $\endgroup$ Commented Oct 25, 2024 at 18:40
  • $\begingroup$ I need the value of all 4 constants in the relation. (A, B, C, D) $\endgroup$ Commented Oct 26, 2024 at 15:34
  • $\begingroup$ @Rose Your algorithm doesn't find the optimal fit ($sse=8.837e5$). MATLAB's "Trust-region" algorithm fits better ($sse=8.824e5$). You also should avoid for-loops in MATLAB. $\endgroup$ Commented Nov 1, 2024 at 6:07
  • $\begingroup$ @gpmath Yes, You are right. I am still looking for more optimal algorithms. Algorithms that can be implemented in other programming languages. I do not understand the algorithm of MATLAB method. $\endgroup$ Commented Nov 1, 2024 at 7:25

4 Answers 4

Reset to default
1
$\begingroup$

In order to compare with your results :

The not iterative method and without initial guess is shown below.

From https://fr.scribd.com/doc/14674814/Regressions-et-equations-integralesand as used in Exponential regression with two terms and constraints

Calculus with MathCad software :

enter image description here

Details of the calculus :

enter image description here

Of course the results depends slightly on the criteria of fitting.

If you prefer a better fit at low values of y instead than at high values it is better to use a criteria of fitting of LMSRE kind than LMSE kind :

enter image description here

If even more accuracy is wanted depending of a particular criteria of fitting one have to use a nonlinear regression software (Thus iterative starting from guessed values of parameters). Then the values computed above are very good initial values.

$\endgroup$
1
  • $\begingroup$ I agree with you that your approach delivers good initial values for a Trust-region or Levenberg-Marquardt minimizer. $\endgroup$ Commented Nov 2, 2024 at 19:23
0
$\begingroup$

Be aware that you cannot determine all constants accurately, especially if there is no guarantee that $y$ is bias-free (it must tend to exactly zero for infinite $t$).

Better work on the logarithm of the signal, which shows a smooth transition between two straight lines of slopes $b$ and $d$. (Check the plot below.) The logarithm can be computed quickly by means of a lookup table.

If the range of $t$ is large enough, you will see straight portions and you can use standard line fitting on these portions only. But if you only see the "knee", there is little that you can do to separate the two terms.

enter image description here

$\endgroup$
0
$\begingroup$

The measured data set has a negative bias, hence we should extend the model formula with bias parameter $E$ to

$$Y = A\cdot e^{B\cdot T}+C\cdot e^{D\cdot T}+E$$

Modification of MATLAB script:

ft = fittype('A*exp(B*x) + C*exp(D*x) + E');
[pop2,gof,output] = fit(T, YN, ft, 'StartPoint', [1e3, -0.010, 1e4, -0.100, -500]);

Result pop2:

 General model:
 pop2(x) = A*exp(B*x) + C*exp(D*x) + E
 Coefficients (with 95% confidence bounds):
   A =        1203  (533.8, 1871)
   B =    -0.01259  (-0.03699, 0.01181)
   C =   1.337e+04  (1.157e+04, 1.516e+04)
   D =     -0.1481  (-0.167, -0.1292)
   E =      -503.2  (-1629, 622.8)

Result gof (goodness of fit):

       sse: 8.824355863574741e+05 (sum-of-squares error)
   rsquare: 0.989753269533616
       dfe: 195
adjrsquare: 0.989543080190715
      rmse: 67.270429605206914

The original approach without parameter E has a poorer goodness of fit:

       sse: 9.857743494119085e+05
   rsquare: 0.988553312881687
       dfe: 196
adjrsquare: 0.988378108487019
      rmse: 70.918692644374715

Measured data plus fit

$\endgroup$
5
  • $\begingroup$ Your answer is excellent. thank you Now, the problem is that I have to write this code in C language. And I need the algorithm of this method. You can help with this too. $\endgroup$ Commented Oct 30, 2024 at 10:09
  • $\begingroup$ @Rose Do you really need C or better C++ or even better C# ? And which library do you want to use for the Levenberg-Marquardt fitting? $\endgroup$ Commented Oct 30, 2024 at 15:56
  • $\begingroup$ @Rose I have a C++ implementation of your curve fitting problem using the ALGLIB library. $\endgroup$ Commented Oct 31, 2024 at 6:24
  • $\begingroup$ I'm curious to see how to implement it with C++ or even C#. Can you share the result here? $\endgroup$ Commented Oct 31, 2024 at 16:24
  • $\begingroup$ @Rose It's a MS Visual Studio 2022 Project: MyCppTest1 $\endgroup$ Commented Oct 31, 2024 at 16:34
0
$\begingroup$

Given a vector of measurements $\boldsymbol{y} \in \mathbb{R}^{n}$:

$$ {y}_{i} = a + b \exp \left( p {t}_{i} \right) + c \exp \left( q {t}_{i} \right) + {w}_{i} $$

Where ${w}_{i}$ is a white noise independent of ${y}_{i}$.

This model fits the Non Linear Least Squares model:

$$ \arg \min_{\boldsymbol{p}} \frac{1}{2} {\left\| {f} \left( \boldsymbol{t} \mid \boldsymbol{p} \right) - \boldsymbol{y} \right\|}_{2}^{2} $$

Or in its scalar form:

$$ \arg \min_{\boldsymbol{p}} \frac{1}{2} \sum_{i = 1}^{n} {\left( {f}_{i} \left( {t}_{i} \mid \boldsymbol{p} \right) - {y}_{i} \right)}^{2} $$

Where ${f}_{i} \left( {t}_{i} \mid \boldsymbol{p} \right) = {p}_{1} + {p}_{2} \exp \left( {p}_{4} {t}_{i} \right) + {p}_{3} \exp \left( {p}_{5} {t}_{i} \right)$ and $p \in \mathbb{R}^{5}$ is the parameters vectors: $\boldsymbol{p} = {\left[ a, b, c, p, q \right]}^{T}$.

It can be solved with many solvers. One of the simplest ones would be the Levenberg Marquardt Algorithm.

Since the problem is not convex, the solution is sensitive to the initial point.
A simple method to get an adequate initialization is given by the answer of @JJacquelin and his paper / book - Jean Jacquelin - Regression and Equations of Integrals.

enter image description here

As can be seen, the initial solution is very similar to the tuned one (Compare the RMSE).
In the graph, the tuned solution, based on the Non Linear LS solver overrides the initial solution.

I implemented both methods in a simple Julia code (The Data CSV is also available).
Given you have access to a Linear Algebra library in your C Code, it will be easy to translate it into C code.

Similar ideas are available in SciKit Guess and Non Linear Least Squares Minimization and Curve Fitting for Python (lmfit).

The code is available on my StackExchange Code GitHub Repository (Look at the Mathematics\Q4989670 folder).

$\endgroup$

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.