Bernoulli Trials

A.I. McLeod
Department of Statistical and Actuarial Sciences,
University of Western Ontario,
e-mail: aim@uwo.ca

B. Quenneville,
Statistics Canada
e-mail: quenne@statcan.ca

We will now examine the performance of these two estimators in the estimation of the parameter p in a sequence of n Bernoulli trials where X is the observed number of successes and p is the probability of success.  The probability function is

(Subscript[f, x])(n, p)= (({{n}, {x}})) (((p^x)(1 - p))^(n - x)) .

So if X successes are observed in n trials, the likelihood function may be written L(p) = (((p^X)(1 - p))^(n - X)) and the mle may be derived by calculus, \(TraditionalForm\`p\&^ = X/n\).  Using Mathematica it is easily shown that the mele of p is \(TraditionalForm\`p\&_ = \((X + 1)\)/\((n + 2)\)\) and that \(TraditionalForm\`\(R(p\&_, p\&^ | p) > 1\ \)\) provided

\(TraditionalForm
      \`p\ ϵ\ 
        \((\(2\ n - \@\(2\ n\^2 + 3\ n + 1\) + 1\)\/\(2\ \((2\ n + 1)\)\), 
          \(2\ n + \@\(2\ n\^2 + 3\ n + 1\) + 1\)\/\(2\ \((2\ n + 1)\)\))
          \)\).

As shown in the figure below, the mele as well as the Bayes estimator is more efficient over most of the range but the relative efficiency tends to 1 as n -> ∞.

In terms of the PMC, the mele and Bayes estimator are better near the center. The figure below also suggests,

\(lim\+\(n -> ∞\)\ PMC(p\&_, \ p\&^\  | p)\) = {({{1, , p = 0.5, }, {(1/2), , p ≠ 0.5, 0, 1, }, {0, , p = 0, 1. }})

Using the Geary-Rao Theorem (Keating, Mason and Sen, 1993, p.103) the result can be established.

In conclusion, if no prior information is available there are no strong grounds for preferring one estimator over another.

BINOMIAL DERIVATIONS

Needs["Statistics`DiscreteDistributions`"]

Derivation of MELE Estimator

Now obtain the mean likelihood estimator, \(TraditionalForm\`p\&_\).

pbar = 
      Divide[Integrate[p*PDF[BinomialDistribution[n,  p], X], 
             {p,  0,  1},  GenerateConditions -> False], 
          Integrate[PDF[BinomialDistribution[n,  p], X],  {p,  0,  1},  
            GenerateConditions -> False]] // FullSimplify

Bayes Estimate with Jeffrey's Prior

The Jeffrey's noninformative prior is (Box and Tiao, p.35)

π[p] = 1/(Sqrt[p (1 - p)])

and hence the posterior distribution for p,

\(TraditionalForm
    \`π(p | X) = \(\(p\^\(X - 1\/2\)\)(1 - p)\)\^\(n - X - 1\/2\)\)

The Bayes estimator of p is then

Jprior[p_] := p^(X - 1/2) (1 - p)^(n - X - 1/2);

numer = 
      Integrate[
          p*(p^X)*(((1 - p))^(n - X))*(p^(X - 1/2) )*
            ((1 - p)^(n - X - 1/2)),  {p,  0,  1},  
          GenerateConditions -> False] // FullSimplify

denom = 
      Integrate[
          (p^X)*(((1 - p))^(n - X))*(p^(X - 1/2) )*
            ((1 - p)^(n - X - 1/2)),  {p,  0,  1},  
          GenerateConditions -> False] // FullSimplify

pBayes = numer/denom  // FullSimplify

Mean-Square Derivations

msebar[p_,  n_]  :=  
      Sum[PDF[BinomialDistribution[n,  p],  k] 
          (p - (k + 1)/(n + 2))^2,  {k, 0, n}]

msehat[p_,  n_]  :=  
      Sum[PDF[BinomialDistribution[n,  p],  k] (p - k/n)^2,  {k, 0, n}]

sols = Solve[FullSimplify[msebar[p, n]/msehat[p, n]] == 1,  p]

sols /. n -> 2 // N

Limit[(p /. sols),  n -> ∞] // N

<<  Graphics`Legend`

g1 = 
      Plot[{msehat[p,  10]/msebar[p,  10], 
          msehat[p,  30]/msebar[p,  30]}, {p, 0, 1},   
        PlotLabel -> "MELE vs MLE", RotateLabel -> False, 
        FrameLabel -> {"p", "R"}, 
        PlotStyle -> {{Thickness[0.001]}, {Thickness[0.001], 
              Dashing[{0.01, 0.02}]}}, Frame -> True, 
        PlotLegend -> {"n=10", "n=30"}, 
        LegendPosition -> {0, -0.37}, LegendShadow -> None, 
        GridLines -> {
            None, {{1, {RGBColor[0.8, 0.8, 0.8], Thickness[0.01]}}}}, 
        PlotRange -> {Automatic, {0, 1.5}}];

Bayes Estimate

mseBayes[p_,  n_]  :=  
      Sum[PDF[BinomialDistribution[n,  p],  k] 
          (p - ((1 + 4 k)/(2 + 4 n)))^2,  {k, 0, n}]

g2 = 
      Plot[{msehat[p,  10]/mseBayes[p,  10], 
          msehat[p,  30]/mseBayes[p,  30]}, {p, 0, 1},   
        RotateLabel -> False, PlotLabel -> "Bayes vs MLE", 
        FrameLabel -> {"p", "R"}, 
        PlotStyle -> {{Thickness[0.001]}, {Thickness[0.001], 
              Dashing[{0.01, 0.02}]}}, Frame -> True, 
        PlotLegend -> {"n=10", "n=30"}, 
        LegendPosition -> {0, -0.37}, LegendShadow -> None, 
        GridLines -> {
            None, {{1, {RGBColor[0.8, 0.8, 0.8], Thickness[0.01]}}}}, 
        PlotRange -> {Automatic, {0, 1.5}}];

Show[GraphicsArray[{g1, g2}]];

Show[GraphicsArray[{{g1}, {g2}}]];

Interval where Bayes estimator better than mle

sols = Solve[FullSimplify[mseBayes[p, n]/msehat[p, n]] == 1,  p]

TeXForm[p /. sols[[1]]] >> 
      "c:\\p\\math\\1999\\mele\\bayes1.tex"

TeXForm[p /. sols[[2]]] >> 
      "c:\\p\\math\\1999\\mele\\bayes2.tex"

sols /. n -> 2 // N

Limit[(p /. sols),  n -> ∞] // N

PMC Derivations: MLE and MELE

Using the definition of \(TraditionalForm\`PMC(p\&_, p\&^, p)\), we can compute it directly.

Needs["Statistics`DiscreteDistributions`"]

PMCScore[p1_,  p2_] := 
      If[Abs[p1 - p] < Abs[p2 - p], 1, 0] + 
        If[Abs[p1 - p] == Abs[p2 - p], 1/2, 0]

PMC[p_,  n_] := 
      Sum[PDF[BinomialDistribution[n,  p],  X] 
          PMCScore[(X + 1)/(n + 2), X/n], {X, 0, n}]

Needs["Graphics`Legend`"]

ga = 
      Plot[{PMC[p, 10], PMC[p, 30]}, {p, 0, 1},   RotateLabel -> False, 
        FrameLabel -> {"p", "PMC"}, PlotLabel -> "MELE vs MLE", 
        PlotStyle -> {{Thickness[0.001]}, {Thickness[0.001], 
              Dashing[{0.01, 0.02}]}}, Frame -> True, 
        PlotLegend -> {"10", "30"}, 
        LegendPosition -> {-0.10, -0.35}, LegendShadow -> None, 
        LegendOrientation -> Horizontal, LegendSize -> {0.5, 0.5}, 
        GridLines -> {
            None, {{0.5, {RGBColor[0.8, 0.8, 0.8], Thickness[0.01]}}}}, 
        PlotRange -> {Automatic, {0, 1}}];

PMC Derivations: Bayes Estimator

PMC2[p_,  n_] := 
      Sum[PDF[BinomialDistribution[n,  p],  X] 
          PMCScore[(4 X + 1)/(4 n + 2), X/n], {X, 0, n}]

gb = 
      Plot[{PMC2[p, 10], PMC2[p, 30]}, {p, 0, 1},   RotateLabel -> False, 
        PlotLabel -> "Bayes vs MLE", FrameLabel -> {"p", "PMC"}, 
        PlotStyle -> {{Thickness[0.001]}, {Thickness[0.001], 
              Dashing[{0.01, 0.02}]}}, Frame -> True, 
        PlotLegend -> {"10", "30"}, 
        LegendPosition -> {-0.10, -0.35}, LegendShadow -> None, 
        LegendOrientation -> Horizontal, LegendSize -> {0.5, 0.5}, 
        GridLines -> {
            None, {{0.5, {RGBColor[0.8, 0.8, 0.8], Thickness[0.01]}}}}, 
        PlotRange -> {Automatic, {0, 1}}];

Show[GraphicsArray[{ga, gb}]];

Show[GraphicsArray[{{ga}, {gb}}]];