# Determing two Nonlinear model fits with 95% confidence bands from one data set

I have the following data which is postulated to have contributions from two different sources .....as can be seen from the confidence band, many points are out of the band. Question: How can I separate the data into two groups such that I can have 2 fitting curves with less data out of confidence band for each cluster?

data = {{17.63, 3794.90}, {19.88, 2410.30}, {21.67, 3282.10}, {23.92,
3153.80}, {24.37, 2846.20}, {24.82, 2487.20}, {25.27,
1692.30}, {26.16, 1564.10}, {27.06, 3205.10}, {27.06,
3692.30}, {27.51, 3384.60}, {29.76, 1307.70}, {30.20,
3615.40}, {31.10, 3282.10}, {33.80, 1179.50}, {40.53,
1051.30}, {47.27, 1871.80}, {47.71, 2589.70}, {49.51,
410.26}, {49.51, 3435.90}, {49.96, 1538.50}, {52.20,
410.26}, {52.65, 948.72}, {52.65, 1179.50}, {54.90,
717.95}, {55.35, 1384.60}, {55.35, 615.38}, {55.80,
846.15}, {56.69, 1897.40}, {57.14, 589.74}, {57.59,
2384.60}, {59.84, 2025.60}, {60.29, 538.46}, {60.29,
358.97}, {62.08, 538.46}, {62.98, 384.62}, {63.43,
1128.20}, {63.43, 615.38}, {63.88, 871.79}, {63.88,
717.95}, {64.78, 461.54}, {64.78, 1615.40}, {65.22,
3743.60}, {67.02, 794.87}, {68.82, 1538.50}, {70.16,
487.18}, {70.16, 897.44}, {72.86, 538.46}, {72.86,
641.03}, {73.76, 435.90}, {74.65, 1538.50}, {74.65,
1205.10}, {75.55, 769.23}, {76.00, 307.69}, {76.45,
410.26}, {76.90, 794.87}, {76.90, 641.03}, {76.90,
538.46}, {80.49, 820.51}, {81.84, 743.59}, {83.18,
512.82}, {86.33, 487.18}, {87.67, 384.62}, {88.12,
1102.60}, {89.02, 871.79}, {94.86, 461.54}, {95.31,
205.13}, {95.76, 333.33}, {96.65, 538.46}, {98.45,
256.41}, {98.45, 461.54}, {98.90, 128.21}, {99.80,
282.05}, {100.24, 666.67}, {100.69, 487.18}, {100.69,
384.62}, {100.69, 153.85}, {102.04, 564.10}, {102.04,
230.77}, {102.94, 974.36}, {102.94, 692.31}, {105.18,
512.82}, {105.18, 282.05}, {105.18, 974.36}, {106.08,
179.49}, {106.53, 1102.60}, {106.98, 615.38}, {107.43,
410.26}, {107.43, 333.33}, {111.92, 384.62}};

nlm = NonlinearModelFit[data, a Exp[-b (x - c)], {a, b, c}, x]

{bands95[x_], bands99[x_], bands999[x_]} =
Table[nlm["MeanPredictionBands",
ConfidenceLevel -> cl], {cl, {.85, .95, .999}}];

Show[ListPlot[data],
Plot[{nlm[x], bands95[x], bands95[x], bands99[x], bands999[x]}, {x,
1, 105}, Filling -> {2 -> {1}, 3 -> {2}, 4 -> {3}, 5 -> {4}}]]


• Can you edit your code properly to make it easier for folks to help you out? Jan 7, 2014 at 1:41
• OK, just inserted a missing line which I left out earlier Jan 7, 2014 at 1:53
• You've posted several other questions on this site, so I'm sure you know how to present code properly. I edited your question and dropped some of the precision of data for purposes of clarity. Jan 7, 2014 at 2:07
• When you say "two different sources" do you mean you wish to model the data with a Exp[-b (x - c)]+ d Exp[-e (x - f)] ? Jan 7, 2014 at 2:39
• So each data point fits EITHER a Exp[-b(x-c) or d Exp[-e(x-f]? If this is the case you may need to come up with some criteria to separate the data, as Mathematica won't know a priori how to bin the data. Jan 7, 2014 at 3:01

The data you represent is cloud-like. Indeed, a brief look at the ListPlot[data]shows that at the approximately same x values the y coordinates may differ 2-3 fold. Such cloudy appearance are typical e.g. for biological data or can be met in polymer physics. In such cases people typically say that the difference by factor 2-3 plays no role, and they are only interested in the exponents, rather than in coefficients.

What I can propose here is to represent the data in the semi-log scale:

    Clear[dataWithin];
Manipulate[
dataWithin =
Select[data, #[[2]] <= a1*Exp[-b1*#[[1]]] && #[[2]] >=
a2*Exp[-b2*#[[1]]] &];
Show[{
ListLogPlot[data, PlotRange -> {100, 10000}],
LogPlot[{a1*Exp[-b1*x], a2*Exp[-b2*x]}, {x, 20, 150},
PlotStyle -> {Red, Darker[Green]}, Filling -> {1 -> {2}}]},
Epilog -> Inset[Column[{
Row[{Style["Points inside the band:   ", 12],
Style[Round[Length[dataWithin]/Length[data]*100 // N, 0.1],
12], Style["%", 12]}],
Row[{Style["Points outside the band: ", 12],
Style[Round[100 - Length[dataWithin]/Length[data]*100 // N,
0.1], 12], Style["%", 12]}]
}], Scaled[{0.7, 0.9}]]
],
{{a1, 9900}, 7000, 10000}, {{b1, 0.022}, 0, 1},
{{a2, 2000}, 1000, 8000}, {{b2, 0.022}, 0, 1}
]


You should see the following on the screen: Most of the points of your data lie within the stripe fixed by the amplitudes of two exponents with the same factor b1=b2=0.022. Play with it. The fraction of points inside and outside the band is shown in the top right corner of the Manipulate panel. Hope it helps.

How about using "SinglePredictionBands" instead of "MeanPredictionBands"?

nlm = NonlinearModelFit[data, a Exp[-b (x - c)], {a, b, c}, x]

{bands95[x_], bands99[x_], bands999[x_]} =
Table[nlm["SinglePredictionBands",
ConfidenceLevel -> cl], {cl, {.85, .95, .999}}];

Show[ListPlot[data],
Plot[{nlm[x], bands95[x], bands95[x], bands99[x], bands999[x]}, {x,
1, 105}, Filling -> {2 -> {1}, 3 -> {2}, 4 -> {3}, 5 -> {4}}]]


See also demonstration: Mean and Single Prediction Bands for a Nonlinear Model

According to the chat mentioned in the comments to the question, the fit required is to the model a*Exp[-b(x-c)] + d*Exp[-e(x-f)], a superposition of two exponential decays. It was expected that one component amplitude would account for 20% of the signal, and the other 80%. The constants Exp[b c] and Exp[e f] may be incorporated into the leading constants a and d, respectively, forming A and B. Then

NonlinearModelFit[data, A*E^(-x*c1) + B*E^(-x*c2), {A,B,c1,c2},
x, MaxIterations->1000]


gives, roughly, 5300*Exp[-0.026 x] + 0.003*Exp[+0.1 x]. The following fit to a single exponential

NonlinearModelFit[data, A*E^(-x*c1), {A,c1}, x, MaxIterations->1000]


gives an almost identical 5100*Exp[-0.025 x]. In other words, there may be good physical reasons to believe these 1980s data are composed of two exponentially decaying sources, but the data do not support such a model. Given the noise and the serious non-orthogonality of exponential decays, fitting to just a single component is the only reasonable course, according to William of Ockham.

Edit

1. The two cluster model may be an artefact of the reduced sampling from about x=30 to x=45. If the sampling density were uniform from x=17 to x=110, would two clusters appear?

2. If two components are required, then what are the amplitude and relaxation time of the small, fast-relaxing component? Consider something analogous to a bootstrap, or a cross-validation, to generate many equivalent data sets. Analyse each with a bi-exponential model. When a thousand such results are combined, the mean+/-stdev values for the large, slow-relaxating component are amplitude: 5120 +/- 77, relaxation time: 40 +/- 0.6. The mean+/-stdev values for the small, fast-relaxating component are amplitude: -0.9 +/- 6, relaxation time: 74 +/- 78. In other words, specifying a second component does not make it measurable. Statistically, the second component is noise.

• If you do a bivariate cluster analysis using the following, there is a trend which shows one source acting at lower values of x, and the other is more influential at higher (> 40) values, so it is not so simple as fitting with just one single component....... Show[SmoothDensityHistogram[data, ColorFunction -> "TemperatureMap", Mesh -> Automatic], ListPlot[FindClusters[data, 2], PlotRange -> {{1, 100}, {0, 4000}}, PlotStyle -> {Black, Red}]] Jan 7, 2014 at 5:12