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.

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]

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]