Geary Rao PMC Formula for the Binomial

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

Needs["Statistics`DiscreteDistributions`"]

PMC: Direct Computation

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

PMC(\(p\&_\), \(p\&^\)| p) = \(Pr { | p\&_ - p | \( < \( | \ p\&^ - p | \)\)\ }\).

we can compute it directly with the Mathematica function,

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

However, due to a minor bug in Mathematica (Version 3.01) this is not computed correctly when p = 0 and n = 0.  The general formula,

PDF[BinomialDistribution[n,  p],  X]

really does not apply,

PDF[BinomialDistribution[n,  0],  0]

This is not a big problem with plotting or evaluating at other values.

Table[PMC[p,  15], {p, 0.05, 0.95, 0.05}]

Geary-Rao Formula

Using Greary-Rao Theorem,

\(TraditionalForm
      \`PMC(p\&_, \ p\&^\  | p)\  = 
        ∑\+\(x = 0\)\%\(x\_1 - 1\)\ \(f\_x\)(n, p) + 
          ∑\+\(x = x\_2 + 1\)\%n\ \(f\_x\)(n, p)\)

where

\(TraditionalForm
      \`x\_1 = ⌈\(2\ p\ \(n(n + 2)\) - n\)\/\(2 
              \((n + 1)\)\)⌉\)

and

(Subscript[x, 2]) = ⌊n/2⌋.

PMCGR[p_,  n_] := 
      Module[{x1 = 
            Ceiling[(2 p (n + 2) n  - n)/(2 (1 + n))], 
          x2 = Floor[n/2]}, 
        If[p <= 1/2, 
          CDF[BinomialDistribution[n, p], x1 - 1] + 1 - 
            CDF[BinomialDistribution[n, p], x2], PMCGR[1 - p, n]]];

TableForm[
      Table[{p, PMC[p, 15], PMCGR[p,  15]}, {p, 0.05, 0.95, 0.05}], 
      TableHeadings -> {None, {"p", "Direct", "Geary-Rao"}}]

As a check, we also use simulation to compute the PMC.

\(PMCsim[p_, \ n_] := 
      Module[{NSIM = 10000, X, errbar, errhat, sd, phat}, 
	
        X = RandomArray[BinomialDistribution[n, p], NSIM]; 
		
        errbar = Abs[p - \((X + 1)\)/\((n + 2)\)]; 
		
        errhat = Abs[p - X/n]; 
	
        phat = Count[Thread[Less[errbar, \ errhat]], \ True]/NSIM // N; 
        sd = Sqrt[p \((1 - p)\)/NSIM]; 
				{phat - 3  sd, phat + 3\ sd}]
        \)

TableForm[
      Table[Flatten[{p, PMCGR[p, 15], PMCsim[p, 15]}], {p, 0.05, 0.95, 
          0.05}], TableHeadings -> {
          None, {"p", "PMC-GearyRao", "est.PMC+3sd", 
            "est.PMC-3sd"}}]