Chapter 6
Biasing for a Fair Return

Code from the Appendix

Kronecker Tensor Product of Matrices

SetAttributes[KroneckerTensor, OneIdentity]; 
    KroneckerTensor[u_?MatrixQ, v_?MatrixQ] := 
      Module[{w = Outer[Times, u, v]}, 
        Partition[
          Flatten[Transpose[
              w, {1, 3, 2, 4}]], Dimensions[
                w][[2]] Dimensions[
                w][[4]]]];
KroneckerTensor[u_, v_, w__] := 
        Fold[KroneckerTensor, u, {v, w}];
CircleTimes = KroneckerTensor;

Supporting Interval Functions

mid[X_] := (Min[X] + Max[X])/2;
extremes[X_] := {Min[X], Max[X]};
diam[X_] := Max[X] - Min[X];

Subdivision of Intervals and Rectangles

subdivide1D[X_] := 
        Interval /@ {{Min[X], mid[X]}, {mid[X], Max[X]}};

Count of Digit Agreement and Form for Pretty Interval Output

DigitsAgreeCount[a_, 
          b_] := (prec = 
            Ceiling@Min[Precision /@ {a, b}]; {{ad, ae}, {bd, 
                be}} = RealDigits[#, 10, prec] & /@ {a, b}; 
          If[ae ≠ be ∨ a b ≤ 0, Return[0]]; 
          If[ad == bd, Return@Length[ad]]; {{com}} = 
            Position[MapThread[Equal, {ad, bd}], False, 1, 1] - 1; 
          com);
DigitsAgreeCount[Interval[{a_, b_}]] := 
        DigitsAgreeCount[a, b];

Interval Arithmetic-Geometric Mean

AGMStep[{a_, b_}] := {(Sqrt[a b]), (1/2) (a + b)}; 
    EnclosingAGM[{a_Real, b_Real}] := 
      Interval@FixedPoint[
          extremes@AGMStep[Interval /@ #] &, {a, 
            b}];
Unprotect[ArithmeticGeometricMean];
ArithmeticGeometricMean[A_Interval, B_Interval] := 
        Block[{Experimental`$EqualTolerance = 0, 
            Experimental`$SameQTolerance = 0}, 
          Interval@
            extremes[EnclosingAGM /@ {Min /@ #, Max /@ #} &@{A, 
                  B}]];
Protect[ArithmeticGeometricMean];

Interval Bisection Method

IntervalBisection[f_, {a_, b_}, tol_] := (X = Interval[{a, b}]; 
        pos = {1}; 
        While[Max[diam /@ X] > tol, 
          pos = Flatten@
              Position[
                f /@ (X = 
                      Join @@ subdivide1D /@ 
                          X), _?(IntervalMemberQ[#, 0] &)]; 
          X = X[[pos]]]; 
        IntervalUnion @@ X)

§6.2 Using Linear Algebra

Matrices[
          n_] := Matrices[n] = (m = 2 n + 1; (Subscript[p, E]) = (1/4) + ϵ0; 
            (Subscript[p, W]) = (1/4) - ϵ0; (Subscript[p, N]) = (1/4); (Subscript[p, S]) = (1/4); 
            Id = SparseArray[{i_, i_} -> 1, {m, m}]; 
            PEW = SparseArray[{{i_, j_} /; j == i + 1 -> 
                    (Subscript[p, E]), {i_, j_} /; j == i - 1 -> (Subscript[p, W])}, {m, 
                  m}]; PNS = 
              SparseArray[{{i_, j_} /; j == i + 1 -> 
                    (Subscript[p, N]), {i_, j_} /; j == i - 1 -> (Subscript[p, S])}, {m, 
                  m}]; {PEW⊗Id + Id⊗PNS, 
              Id⊗Id});
\(ReturnProbability[ϵ_Real, 
        n_Integer] := \(({A, Id} = Matrices[n]; 
        Block[{ϵ0 = ϵ}, m = 2\ n + 1; ctr = n\ m + n + 1; 
          r = A[[All, ctr]]; 
          A[[All, ctr]] = 0; 
          q = LinearSolve[Id - A, r]; 
          q[[ctr]]])\); \({n == #, 
            N == \((2  # + 1)\)\^2 - 
                1, \(ϵ\_*\) == \((ϵ /. 
                  FindRoot[
                    ReturnProbability[ϵ, #] == 
                      1\/2, {ϵ, 0.06, 0.07}, 
                    WorkingPrecision -> MachinePrecision, 
                    PrecisionGoal -> 12\ ])\)} &\) /@ {10, 20, 40} // 
      TableForm\)

§6.4 Reducing the Complexity

§6.4.1 Using Combinatorics

\(ExpectedVisits[ϵ_Real, K_] := \((p\_E = 0.25 + ϵ; 
        p\_W = 0.25 - ϵ; p\_N = \(p\_S = 0.25\); w = \(p\_E\) p\_W; 
        z = \(p\_N\) p\_S; \((a = 1)\) + 
          Sum[\((a *= \(\(2  k \((2  k - 1)\)\)\/k\^2\) z)\) \((\((b = 1)\) + 
                  Sum[b *= \(\((\(k - j + 1\)\/j)\)\^2\) w\/z, {j, 1, 
                      k}])\), {k, 1, 
              K - 1}])\); \({K == #, \(ϵ\_*\) == \((ϵ /. FindRoot[ExpectedVisits[ϵ, #] == 2, {ϵ, 0.06, 
                      0.07}, WorkingPrecision -> MachinePrecision, 
                    PrecisionGoal -> 12\ ])\)} &\) /@ {125, 250, 500, 
          1000} // TableForm\)

§6.4.2 Using Symbolic Computation

ExpectedVisits[ϵ_Real] := ((Subscript[p, E]) = (1/4) + ϵ; 
          (Subscript[p, W]) = (1/4) - ϵ; (Subscript[p, N]) = (Subscript[p, S]) = (1/4); 
          w = (Subscript[p, E]) (Subscript[p, W]); 
          z = (Subscript[p, N]) (Subscript[p, S]); 
s = (Subscript[p, old]) = 1; 
          (Subscript[p, new]) = 2 (w + z); k = 2; 
          While[s =!= (s += (Subscript[p, new])), 
            (Subscript[p, new]) = ((-4 (2 k - 1) (2 k - 3) (((w - z))^2) (Subscript[p, old]) + 2 (((2 k - 1))^2) (w + z) ((Subscript[p, old]) = (Subscript[p, new])))/((k++)^2))]; s);
\(\({Precision == #, \(ϵ\_*\) == \ \((ϵ /. 
                  FindRoot[
                    ExpectedVisits[ϵ] == 2, {ϵ, 0.06, 
                      0.07}, WorkingPrecision -> #, 
                    AccuracyGoal -> #])\)} &\) /@ {13, MachinePrecision, 
          19, 22, 25} // TableForm\)

§6.5 The Joy of Special Functions

LegendreKernel[x_, 
          k_] := (1/π) (((x + (Sqrt[(x^2) - 1]) Cos[θ]))^k);
Integrate[
      Sum[Binomial[2  k, k] 
          LegendreKernel[((((Subscript[p, EW])^2) + ((Subscript[p, NS])^2))/(((Subscript[p, EW])^2) - ((Subscript[p, NS])^2))), 
            k] (((((Subscript[p, EW])^2) - ((Subscript[p, NS])^2)))^k), {k, 0, ∞}], {θ, 
        0, π}, 
      Assumptions -> {0 ≤  (Subscript[p, EW]) ≤  (Subscript[p, NS]), 
          2 ((Subscript[p, EW]) + (Subscript[p, NS])) < 1}]

\(ExpectedVisits[ϵ_] := 
      1/ArithmeticGeometricMean[\@\(1 - 4 \((p\_EW + p\_NS)\)\^2\), \@\(1 - 4 \((p\_EW - p\_NS)\)\^2\)] /. {p\_EW -> \ \@\(\((1\/4 + ϵ)\) \((1\/4 - ϵ)\)\), 
          p\_NS -> \ 
            1\/4}; \({Precision == #, \(ϵ\_*\) == \ \((ϵ /. FindRoot[
                    ExpectedVisits[ϵ] == 2, {ϵ, 0.06, 
                      0.07}, WorkingPrecision -> #, 
                    AccuracyGoal -> #])\)} &\) /@ {13, MachinePrecision, 
          19, 22, 25} // TableForm\)

§6.5.1 Using Interval Arithmetic

prec = 
      10010; ϵTry = ϵ /. 
        FindRoot[ExpectedVisits[ϵ] - 2, {ϵ, 0.06, 0.07}, 
          WorkingPrecision ->  prec, 
          AccuracyGoal -> prec]; ϵInterval = 
      IntervalBisection[
        ExpectedVisits[#] - 2 &, ϵTry + {-1, 1} 
            (10^(-prec + 5)), (10^(-prec))];
DigitsAgreeCount[ϵInterval]

§6.7 Harder Problems

Hypercubic Lattices.

Scaled Modified Bessel Function

The Theory

presets...

Off[InverseFunction::ifun, Solve::ifun]; << 
      Graphics`Graphics`
thick = 0.01;
{Red, Green, 
          Blue} = {{Thickness[thick], RGBColor[1, 0, 0]}, {Thickness[thick], 
            RGBColor[0, 1, 0]}, {Thickness[thick], 
            RGBColor[0, 0, 1]}};
SetOptions[Plot, 
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 12}, ImageSize -> 300, 
        PlotStyle -> {Blue, Red, Green}];
SetOptions[LogPlot, 
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 12}, ImageSize -> 300, 
        PlotStyle -> {Blue, Red, Green}];

The scaled modified Bessel function (e^(-x )) (Subscript[I, 0])(
          x) of the first kind of order zero is calculated with the help of an  asymptotic  expansion with error bounds that can be found at [Olv74, p. 269, Ex. 13.2]. Specifically, we obtain for x > 0 that

\(TraditionalForm\`e\^\(-x\)\ \(\(I\_0\)(
            x)\) = \(1\/\@\(2  π\ x\)\) \((∑\+\(k = 0\)\%\(n - 1\)\((\(\((2\ k - 1)\)!!\))\)\^2\/\(\(k!\)\ \((8\ x)\)\^k\) + R\_n)\)\)

where, with χ(
          n) = (π^(1/2)) Γ((1/2) n + 
                1)/Γ((1/2) n + (1/2)),

\(TraditionalForm\`\(\(|\)\(R\_n\)\(|\)\(\(≤\)\(\ \)\(2 \( χ(n)\) e\^\(π/\((8  x)\)\)\ \(\((2  n - 1)\)!!\)\^2\/\(\(n!\)\ \((8  x)\)\^n\)\)\)\) = 
        2\ e\^\(π/\((8  x)\)\)\ \(Γ(n + 1\/2)\)\^2\/\(\(Γ(\(1\/2\) \((n + 1)\))\)\^2\ \((4  x)\)\^n\)\).

derive the last equality...

FullSimplify[
        2  (π^(1/2)) Γ((1/2) n + 
                1)/Γ((1/2) n + 
                (1/2)) (E^(π/(8  x))) ((((2  n - 1)!!)^2)/(n! (((8  x))^n))), {x > 0, n > 0, n ∈ Integers}] // 
      TraditionalForm

Stirling's formula yields the asymptotics

\(TraditionalForm\`\(Γ(n + 1\/2)\)\^2\/\(\(Γ(\(1\/2\) \((n + 1)\))\)\^2\ \((4  x)\)\^n\) ∼ \((\(2  n\)\/e)\)\^n\) ,

check this asymptotics...

Limit[(((Γ(n + (1/2)))^2)/(((Γ((1/2) (n + 
                        1)))^2) ))/(((((2  n)/E)))^n), n -> ∞]

in fact by monotonicity

\(TraditionalForm\`\(Γ(n + 1\/2)\)\^2\/\(\(Γ(\(1\/2\) \((n + 1)\))\)\^2\ \((4  x)\)\^n\) ≤ \ \(\(e\ π\)\/8\) \((\(2  n\)\/e)\)\^n\),	n ≥ 1.

check this claim...

(((Γ(n + (1/2)))^2)/(((Γ((1/2) (n + 
                        1)))^2) ))/(((((2  n)/E)))^n) /. n -> 1

Plot[((E π)/8) (((((2  n)/E)))^n)/(((Γ(n + (1/2)))^2)/(((Γ((1/2) (n + 1)))^2) )),{n, 1, 100},AxesLabel->{"n","ratio"}]

Summarizing, we have

\(TraditionalForm\`\(\(|\)\(R\_n\)\(|\)\(\(≤\)\(\(\(e\ π\)\/4\) e\^\(π/\((8  x)\)\)\ \((n\/\(2  e\ x\))\)\^n\)\(≤\)\(\ \)\(e\ \((n\/\(2  e\ x\))\)\^n\)\)\)\),	x ≥ 2.

check the last inequality...

((E π)/4) (E^(π/(8  x))) ≤  E /. x -> 2

Thus, we can guarantee the error \(TraditionalForm\`\(\(|\)\(R\_n\)\(|\)\(\(≤\)\(\ \)\(e\ ·10\^\(-d\)\)\)\)\), which essentially guarantees d correct decimal digits, if for a given x the summation index n is chosen in accordance to the solution \(TraditionalForm\`\(\(η\_*\)\(\ \)\)\)of the equation

\(TraditionalForm\`\(\(\((η\/\(2  e\ x\))\)\^η\)\(==\)\(10\^\(-d\)\)\(\ \)\)\),	that is 	\(TraditionalForm\`η\ \ \(log(\(2  e\ x\)\/η)\) == d\ log\ 10\( = \&def\)\ 2  x\_0\) .

Generally, if any, there will be two solution that a grouped around a minimum error. This behavior is typical for asymptotic expansions that do not converge as n -> ∞.

have a look at a typical situation...

LogPlot[((((η/(2  E x))))^η) /. 
        x -> 60. , {η, 3. , 300. }, 
      AxesLabel -> {"η", "error"}]

It turns out that the minimum is located at (Subscript[η, min]) = 2  x.

check this claim...

Solve[
      D[((((η/(2  E x))))^η), η] == 
        0, η]

We obtain

\(TraditionalForm\`\(η\_*\) = \(2\ x\_0\)\/\(-\(W(\(-x\_0\)/\((e\ x)\))\)\)\),

calculate the solution \(TraditionalForm\`\(η\_*\)\)

Solve[η  log(((2  E x)/η)) == 2  (Subscript[x, 0]), η]

where W denotes Lambert's W-function. There are exactly two branches of -W(-z) that are real for z ∈ (0, (e^(-1))]. We choose the smaller solution \(TraditionalForm\`\(η\_*\)\), that is, the branch for which -W(-z) ≥  1,  z\(TraditionalForm\`\(\(∈\)\((0, e\^\(-1\)]\)\)\). (The other has -W(-z) ≤ 
        1, z∈(0, (e^(-1))] .) In particular, -W(-(e^(-1))) = 1 and -W(-z)  ->  +∞ as z ↘ 0.This branch is given in Mathematica by -ProductLog[-1,-z]. Since -W(-z) is complex for z > (e^(-1)), we must restrict ourselves to \(TraditionalForm\`x\_0/\((e\ x)\) ≤ e\^\(-1\)\), that is x ≥ (Subscript[x, 0]).

have a look at the branches...

Plot[{-ProductLog[-1, -z], -ProductLog[0, -z], 
          1},{z, 0, (E^(-1))},AxesLabel->{"z","-W(-z)"}]

With that choice of the branch we have (Subscript[η, * ]) ≤  
          2  (Subscript[x, 0]) ≤  2  x = (Subscript[η, min]). Thus, the error is descending here and we choose the summation index n as the ceiling of \(TraditionalForm\`\(η\_*\)\),

\(TraditionalForm\`n = ⌈\(2\ x\_0\)\/\(-\(W(\(-x\_0\)/\((e\ x)\))\)\)⌉\),		x ≥ (Subscript[x, 0]) = (d/2) log 10.

The Code

BesselITilde[x_?Positive] := 
      Module[{x0 = Precision[x] Log[10]/2, n}, 
        If[x < x0, Exp[-x] BesselI[0, x], 
          n = Ceiling[
              2  x0/-ProductLog[-1, -x0/(E x)]]; (1/(Sqrt[2  π x])) 
            Sum[(( ((((2  k - 1)!!))^2))/(k! (((8  x))^k))), {k, 0, n - 1}]]]

a test...

BesselITilde[63.`50]

(E^(-x)) (Subscript[I, 0])(x) /. 
      x ->  63.`52

Solution of Problem 6 (once more)

ExpectedVisits[ϵ_Real] := ((Subscript[p, E]) = (1/4) + ϵ; 
          (Subscript[p, W]) = (1/4) - ϵ; (Subscript[p, N]) = (Subscript[p, S]) = (1/4); 
          (Subscript[p, EW]) = (Sqrt[(Subscript[p, E]) (Subscript[p, W])]); 
          (Subscript[p, NS]) = (Sqrt[(Subscript[p, N]) (Subscript[p, S])]); σ = 2 ((Subscript[p, EW]) + (Subscript[p, NS])); 
          NIntegrate[
            Exp[-(1 - σ) t] BesselITilde[2  (Subscript[p, EW]) t] 
              BesselITilde[2  (Subscript[p, NS]) t], {t, 0, ∞}, 
            PrecisionGoal ->  15, 
            Method ->  DoubleExponential]);
FindRoot[
      ExpectedVisits[ϵ] == 2, {ϵ, 0.06, 0.07}, 
      AccuracyGoal -> 15]

Table 6.1

ReturnProbability[d_Integer] := 
      1 - 1/NIntegrate[BesselITilde[t/d]^d, {t, 0, ∞}, 
            WorkingPrecision -> 20, PrecisionGoal ->  17, 
            Method ->  DoubleExponential]
{#, NumberForm[ReturnProbability[#], {15, 15}]} & /@ 
        Join[Range[3, 20], {30, 40, 100}] // TableForm