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

Question

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}}]]

Answer 1

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.

Answer 2

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

Answer 3

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.