This method uses single prediction confidence intervals to determine and select outliers. The confidence levels are set to 1, 2 and 3 standard deviations. The points outside 3 SD can be found in o3 and the points outside 1 SD can be found in o1 and o3 combined. The points in o1 and o3 are plotted in the chart in green and red respectively.
data = Transpose@{experimentalCR, theoreticalCR};
{mx, my} = 1.1*Max /@ {experimentalCR, theoreticalCR};
{sd1, sd2, sd3} =
2*(CDF[NormalDistribution[0, 1], #] - 0.5) & /@ {1, 2, 3};
lm = LinearModelFit[data, {1, x }, x];
getseries[sd_] := Module[{cb},
cb = lm["SinglePredictionBands", ConfidenceLevel -> sd];
{lower, upper} = Transpose[cb /. x -> # & /@ experimentalCR];
p1 = Position[
MapThread[#1 <= #2 <= #3 &, {lower, theoreticalCR, upper}], True];
Extract[data, p1]]
s3 = getseries[sd3];
o3 = Complement[data, s3];
s1 = getseries[sd1];
o1 = Complement[s3, s1];
{bands68[x_], bands95[x_], bands99[x_]} = Table[lm["SinglePredictionBands",
ConfidenceLevel -> cl], {cl, {sd1, sd2, sd3}}];
Show[ListPlot[data],
Plot[{lm[x], bands68[x], bands95[x], bands99[x]},
{x, 0, mx}, Filling -> {2 -> {1}, 3 -> {2}, 4 -> {3}}],
ListPlot[o1, PlotStyle -> Directive[Green, PointSize[Large]]],
ListPlot[o3, PlotStyle -> Directive[Red, PointSize[Large]]],
AxesOrigin -> {0, 0}, PlotRange -> {{0, mx}, {0, my}},
ImageSize -> 480, Frame -> True]
Edit
It seems appropriate to add an alternative method, (although the data in this case does not suggest the suitability of a multivariate fit):-
data = Transpose@{experimentalCR, theoreticalCR};
prange = Sort[#][[{1, -1}]] & /@ {experimentalCR, theoreticalCR};
{{xmin, xmax}, {ymin, ymax}} = {#1, #2*1.35} & @@@ prange;
(* For values within two standard deviations,(approx 95.45% of values) *)
sd = 2;
cl = 2*(CDF[NormalDistribution[0, 1], sd] - 0.5);
Needs["MultivariateStatistics`"];
e = EllipsoidQuantile[data, cl];
ctr = e[[1]];
{r1, r2} = e[[2]];
inc = ArcTan[e[[3, 1, 2]]/e[[3, 1, 1]]]*180/Pi;
Print["Ellipse center = " <> ToString@ctr];
Print["Ellipse radii (r1, r2) = " <> ToString@{r1, r2}]; Print[
StringJoin["Ellipse inclination = ", ToString@inc, " degrees"]];
(* Find the foci of the ellipse *)
f = Sqrt[r1^2 - r2^2];
dx = f*Cos[inc Degree];
dy = f*Sin[inc Degree];
f1 = ctr - {dx, dy};
f2 = ctr + {dx, dy};
edge = ctr + r1*e[[3, 1]];
rlim = EuclideanDistance[edge, f1] + EuclideanDistance[edge, f2];
(* nod to belisarius here *)
inside[{x_, y_}, {f1_, f2_}] :=
Sum[EuclideanDistance[{x, y}, i], {i, {f1, f2}}];
sd = Select[data, inside[#, {f1, f2}] < rlim &];
Show[RegionPlot[inside[{x, y}, {f1, f2}] < rlim,
{x, xmin, xmax}, {y, ymin, ymax}],
ListPlot[data], Graphics[{Green, Point@sd}],
Graphics@e,
Graphics[{Black, Thick, Dashing[0.05],
Rotate[Circle[ctr, {r1, r2}], inc Degree]}],
Graphics[{Red, Line[{ctr + r1*e[[3, 1]], ctr, ctr + r2*e[[3, 2]]}]}],
Graphics[{Yellow, PointSize[Large], Point[{f1, f2}]}],
PlotRange -> {{xmin, xmax}, {ymin, ymax}},
AspectRatio -> (ymax - ymin)/(xmax - xmin), ImageSize -> 300]
