How to fit a biexponential decay function

I am having trouble fitting this biexponential decay function, any suggestions? Maybe lsqcurvefit is not the best for this purpose, I am not sure.

Where D1 through D4 are unknown fitting parameters, B is given x axis and B0 is the initial known B.

As some help, D1+D3 should approximately be = 1.

clear; clc; clf; close all;
%%
B = [23.11277703	27.4672621	55.28799957	77.30318258	128.4435712	144.59607	150.1007715	178.348134	178.3800029	186.8782579	262.1289782	289.5852833	387.9146522	404.3098886	414.7168762	438.7555353	460.9984083	463.9649687	622.263271	637.1680637	736.5544189	796.517528	823.94181	834.148104	875.3224784	883.5748864	1120.051524	1135.690878	1196.020072	1260.164236	1372.270872	1417.958217	1421.352213	1437.177887	1738.235663	1766.311611	1802.411799	1829.695659	2059.704061	2099.087613	2124.773971	2155.740229	2447.429036	2491.720482	2505.111798	2622.426034	2886.241378	2908.528677	2946.363137	3047.494139	3239.372347	3286.412652	3380.505981	3595.733583	3849.675405	3851.882823	3901.945387	4093.219946	4142.141544	4173.598221	4404.59216	4722.334446	4922.527799	4956.628395	4991.520719	5155.736628	5166.668506	5292.917652	5563.979018	6002.228624	6127.085857	6200.478095	6215.089135	6265.623506	6280.157597	6646.587256	6858.666555	7435.416116	7463.34958	7470.463221	7515.404453	7572.650633	7583.431923	8154.228757	8288.654772	8781.187652	8861.477195	8931.318967	9021.896921	9064.205214	9105.489879	9815.842157	9853.943669	10197.7968	10318.37582	10530.99402	10689.75288	10761.67104	10766.65196	11554.53325	11631.42746	11720.29066	11886.10034	12262.37474	12449.29363	12566.91817	12654.73847	13348.66924	13390.4235	13564.65074	13600.98465	14125.46112	14342.82746	14506.28851	14701.09922	15082.93253	15354.02702	15361.61444	15724.51375	16370.35437	16584.76298	16900.75329	16923.08054	17254.2292	17468.10605	18002.01474	18531.87437	18802.34157	18869.11326	19253.70066	19709.89835	20433.48763	20827.38744	20921.0307	21159.02429	22086.99132	23018.93242	23078.83285	23256.89361	23654.81114	24599.38498	25758.3491	25820.39285	26289.70212	27247.07931	28517.88518	28651.73769	30030.07432	31699.09817	34900.43055	38255.73483];
y = [1	1	1	1	1	1	0.931938171	1	0.956559088	1	0.870266512	0.847039689	0.819706047	0.828390985	0.820051375	0.817313288	0.859815135	0.854877764	0.690813826	0.654914521	0.683075897	0.638696623	0.622781558	0.630439319	0.724997201	0.655188661	0.474937192	0.507599749	0.53700173	0.487243326	0.452908865	0.453747137	0.587406267	0.486231391	0.332199872	0.403521204	0.347536005	0.370257005	0.330503095	0.450762837	0.35864512	0.327000256	0.28955166	0.232696502	0.281911772	0.230350445	0.245922383	0.325485137	0.246350406	0.239027624	0.191768251	0.212565359	0.160778372	0.14856721	0.222365205	0.178102728	0.196251022	0.171584923	0.129942067	0.155410041	0.106192265	0.093225369	0.142592864	0.119221524	0.136018225	0.082431061	0.107723572	0.1139216	0.061249575	0.056398596	0.086408465	0.07132462	0.08647058	0.069306234	0.050232351	0.067788085	0.036482203	0.032246077	0.049572052	0.041494183	0.029684983	0.049760318	0.038130662	0.036214727	0.018926179	0.022604606	0.017131299	0.027475932	0.017451191	0.026141832	0.018448348	0.017620037	0.009337477	0.012321971	0.009783651	0.014762149	0.012895101	0.009050812	0.008452548	0.004582352	0.008153299	0.006376435	0.005740571	0.008224445	0.006262219	0.003883995	0.004694069	0.003368699	0.002409039	0.003584668	0.003811612	0.005099296	0.003120894	0.001951461	0.002628413	0.001891251	0.002514979	0.001513604	0.001975856	0.001713393	0.001142865	0.00167866	0.001150661	0.001987164	0.001172911	0.001219744	0.001045124	0.000792151	0.000740858	0.001330426	0.000980825	0.000928473	0.000732892	0.000567673	0.000572429	0.000921887	0.000734829	0.000361819	0.000473032	0.000472086	0.000804721	0.000584736	0.00038979	0.000361897	0.000726435	0.000697616	0.000535132	0.000676624	0.000467918	0.00036719	0.000318515];
D_guess = [0.5 0.001 0.5 0.001];

figure
fun = @(D,B) D(1)*exp(-(B-B(1)).*D(2))+D(3)*(eBp(-(B-B(1)).*D(4)));
[fit_y, rsb_y] = lsqcurvefit(fun, D_guess, B, y);
semilogy(B, y, 'ko', B, fun(fit_y,B), 'r-')

Answer 1

[John D'Errico]

First, NEVER set up a solver to fit exponential models with the SAME initial rate parameters for both terms!!!!!!!! So this is flat out BAD:

D_guess = [0.5 0.001 0.5 0.001];

That will result in an immediate problem for the solver. You can have them the same for the amplitude terms out front, but the rate terms SHOULD NOT BE THE SAME.

B = [23.11277703	27.4672621	55.28799957	77.30318258	128.4435712	144.59607	150.1007715	178.348134	178.3800029	186.8782579	262.1289782	289.5852833	387.9146522	404.3098886	414.7168762	438.7555353	460.9984083	463.9649687	622.263271	637.1680637	736.5544189	796.517528	823.94181	834.148104	875.3224784	883.5748864	1120.051524	1135.690878	1196.020072	1260.164236	1372.270872	1417.958217	1421.352213	1437.177887	1738.235663	1766.311611	1802.411799	1829.695659	2059.704061	2099.087613	2124.773971	2155.740229	2447.429036	2491.720482	2505.111798	2622.426034	2886.241378	2908.528677	2946.363137	3047.494139	3239.372347	3286.412652	3380.505981	3595.733583	3849.675405	3851.882823	3901.945387	4093.219946	4142.141544	4173.598221	4404.59216	4722.334446	4922.527799	4956.628395	4991.520719	5155.736628	5166.668506	5292.917652	5563.979018	6002.228624	6127.085857	6200.478095	6215.089135	6265.623506	6280.157597	6646.587256	6858.666555	7435.416116	7463.34958	7470.463221	7515.404453	7572.650633	7583.431923	8154.228757	8288.654772	8781.187652	8861.477195	8931.318967	9021.896921	9064.205214	9105.489879	9815.842157	9853.943669	10197.7968	10318.37582	10530.99402	10689.75288	10761.67104	10766.65196	11554.53325	11631.42746	11720.29066	11886.10034	12262.37474	12449.29363	12566.91817	12654.73847	13348.66924	13390.4235	13564.65074	13600.98465	14125.46112	14342.82746	14506.28851	14701.09922	15082.93253	15354.02702	15361.61444	15724.51375	16370.35437	16584.76298	16900.75329	16923.08054	17254.2292	17468.10605	18002.01474	18531.87437	18802.34157	18869.11326	19253.70066	19709.89835	20433.48763	20827.38744	20921.0307	21159.02429	22086.99132	23018.93242	23078.83285	23256.89361	23654.81114	24599.38498	25758.3491	25820.39285	26289.70212	27247.07931	28517.88518	28651.73769	30030.07432	31699.09817	34900.43055	38255.73483];
y = [1	1	1	1	1	1	0.931938171	1	0.956559088	1	0.870266512	0.847039689	0.819706047	0.828390985	0.820051375	0.817313288	0.859815135	0.854877764	0.690813826	0.654914521	0.683075897	0.638696623	0.622781558	0.630439319	0.724997201	0.655188661	0.474937192	0.507599749	0.53700173	0.487243326	0.452908865	0.453747137	0.587406267	0.486231391	0.332199872	0.403521204	0.347536005	0.370257005	0.330503095	0.450762837	0.35864512	0.327000256	0.28955166	0.232696502	0.281911772	0.230350445	0.245922383	0.325485137	0.246350406	0.239027624	0.191768251	0.212565359	0.160778372	0.14856721	0.222365205	0.178102728	0.196251022	0.171584923	0.129942067	0.155410041	0.106192265	0.093225369	0.142592864	0.119221524	0.136018225	0.082431061	0.107723572	0.1139216	0.061249575	0.056398596	0.086408465	0.07132462	0.08647058	0.069306234	0.050232351	0.067788085	0.036482203	0.032246077	0.049572052	0.041494183	0.029684983	0.049760318	0.038130662	0.036214727	0.018926179	0.022604606	0.017131299	0.027475932	0.017451191	0.026141832	0.018448348	0.017620037	0.009337477	0.012321971	0.009783651	0.014762149	0.012895101	0.009050812	0.008452548	0.004582352	0.008153299	0.006376435	0.005740571	0.008224445	0.006262219	0.003883995	0.004694069	0.003368699	0.002409039	0.003584668	0.003811612	0.005099296	0.003120894	0.001951461	0.002628413	0.001891251	0.002514979	0.001513604	0.001975856	0.001713393	0.001142865	0.00167866	0.001150661	0.001987164	0.001172911	0.001219744	0.001045124	0.000792151	0.000740858	0.001330426	0.000980825	0.000928473	0.000732892	0.000567673	0.000572429	0.000921887	0.000734829	0.000361819	0.000473032	0.000472086	0.000804721	0.000584736	0.00038979	0.000361897	0.000726435	0.000697616	0.000535132	0.000676624	0.000467918	0.00036719	0.000318515];

Next, lets plot the data. How good is it? You have much data, but it seems pretty noisy.

plot(B,y,'.')

And that tells me little. Except what I should have known to do in the first place. (I did include the first plot, just to show my thinking.)

semilogy(B,y,'.')
grid on

Again, pretty darn noisy. It almost seems as if you have multiple sets of data overlaid on top of each other. But the first part is pretty linear in the logged plot. And that suggests we will be able to model at least one term well.

Now to try a fit.

min(B)

ans = 23.1128

mdl = fittype('D1*exp(-(B-23.11277703)*D2) + D3*exp(-(B-23.11277703)*D4)','indep','B')

mdl = 
     General model:
     mdl(D1,D2,D3,D4,B) = D1*exp(-(B-23.11277703)*D2) + D3*exp(-(B-23.11277703)*D4)

Dstart = [.5 .001 .5 .0001];
fittedmdl = fit(B',y',mdl,'start',Dstart)

fittedmdl = 
     General model:
     fittedmdl(B) = D1*exp(-(B-23.11277703)*D2) + D3*exp(-(B-23.11277703)*D4)
     Coefficients (with 95% confidence bounds):
       D1 =      0.4209  (0.1632, 0.6786)
       D2 =     0.00106  (0.0005892, 0.00153)
       D3 =      0.6216  (0.3566, 0.8867)
       D4 =   0.0003579  (0.0002785, 0.0004373)

Bplot = min(B):max(B);
ypred = fittedmdl(Bplot);

plot(B,y,'b.')
hold on
plot(Bplot,ypred,'r-')
hold off

And that seems to fit reasonably well. But now log the axes, and we see:

semilogy(B,y,'b.')
hold on
semilogy(Bplot,ypred,'r-')
hold off

Pretty crappy here. The fit is terrible in the log plot. It suggests again something I knew before, but I need to motivate what I am doing. Otherwise, I'll just get peppered with questions later on anyway.

We should be doing the fit in a log domain. Except, if we do that, then we have numerical problems, because we will be taking the logs of negative numbers at times. (Trust me here.)

So instead, I'll use weights that will mimic the same process.

W = 1./y';
fittedmdl2 = fit(B',y',mdl,'start',Dstart,'weight',W)

fittedmdl2 = 
     General model:
     fittedmdl2(B) = D1*exp(-(B-23.11277703)*D2) + D3*exp(-(B-23.11277703)*D4)
     Coefficients (with 95% confidence bounds):
       D1 =      0.3267  (0.2287, 0.4247)
       D2 =    0.001331  (0.0008104, 0.001852)
       D3 =      0.7221  (0.6181, 0.8261)
       D4 =   0.0003938  (0.0003696, 0.000418)

ypred = fittedmdl2(Bplot);

plot(B,y,'b.')
hold on
plot(Bplot,ypred,'r-')
hold off

semilogy(B,y,'b.')
hold on
semilogy(Bplot,ypred,'r-')
hold off

Still, the fit is crapola down low. And that suggests your model does not expllain the data. This is apparently not well fit by a sum of exponentials. There is other stuff happening in the tail, something different happening down low. Possibly that is just noise, and we should discard all the data above B==15000 or so. I'll try 20k.

k = B<=20000;

fittedmdl3 = fit(B(k)',y(k)',mdl,'start',Dstart,'weight',W(k))

fittedmdl3 = 
     General model:
     fittedmdl3(B) = D1*exp(-(B-23.11277703)*D2) + D3*exp(-(B-23.11277703)*D4)
     Coefficients (with 95% confidence bounds):
       D1 =      0.3233  (0.2197, 0.4268)
       D2 =    0.001345  (0.000779, 0.00191)
       D3 =      0.7258  (0.616, 0.8356)
       D4 =   0.0003948  (0.0003691, 0.0004205)

Bplotk = min(B):20000;
ypredk = fittedmdl2(Bplotk);

semilogy(B(k),y(k),'b.')
hold on
semilogy(Bplotk,ypredk,'r-')
hold off

Finally, that seems quasi-reasonable. The amplitude parameters D(1) and D(3) sum to quite close to 1.

Answer 2

[Walter Roberson]

To be honest, you should probably give up on trying to use lsqcurvefit to fit that model, unless you happen to have reason to know quite good approximations of the model parameters.

The sum of exponentials is notoriously difficult to fit using least squared approaches. There is a very large tendency for one of the exponentials to become very wide, effectively a constant line or fairly slight slope, and for the other exponential to be left to try to do all of the fitting work.

This is especially true if the coefficients are not constrained: with the forms of the two terms being the same, exchanging the three model parameters between the two terms gives you the same function value, so the fitting cannot possibly decide whether (for example) it is the first term or the second term that should be multiplied by 0.83. When you are doing a fitting you either need the terms to be biased or you need constraints to prevent the proposed coefficients from swapping places.

You might perhaps be able to use the curve fitting app asking for sum of two exponentials. Or use some other approach.

A year or two ago I did see a paper that proposed a way to fit sum of exponentials, but I could not tell how practical it was.

In the short term, you could consider ga or other evolutionary algorithm approaches, in that they run many different model parameter sets simultaneously, increasing the chances that at least one set will happen to land in the sweet spot.

[John D'Errico]

I tried a few things, not all of them shown in my answer. The data does not seem to be fully fit well by that model, though with some effort, I did manage to find a fit to the first 2/3 of the curve. Down in the bottom end it seems like noise (or some extraneous factor) is corrupting the data beyond the point that the data in the tail is of any value for the model.