This method uses [single prediction confidence intervals][1] 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]

![enter image description here][2]

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

![enter image description here][3]


  [1]: http://reference.wolfram.com/mathematica/tutorial/StatisticalModelAnalysis.html
  [2]: https://i.sstatic.net/1LiSu.png
  [3]: https://i.sstatic.net/XSrYt.png