Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

Mathematica graphics

share|improve this question
    
Can you edit your code properly to make it easier for folks to help you out? –  bobthechemist Jan 7 at 1:41
    
OK, just inserted a missing line which I left out earlier –  thils Jan 7 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 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 at 2:39
2  
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 at 3:01

3 Answers 3

up vote 5 down vote accepted

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: enter image description here 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.

share|improve this answer

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

enter image description here

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

share|improve this answer

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.

share|improve this answer
    
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 at 5:12

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.