# 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? – bobthechemist Jan 7 '14 at 1:41
OK, just inserted a missing line which I left out earlier – thils Jan 7 '14 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. – bobthechemist Jan 7 '14 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)] ? – bobthechemist Jan 7 '14 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. – bobthechemist Jan 7 '14 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}]] – thils Jan 7 '14 at 5:12