Biometrika Trust

Sampling Theory of the Negative Binomial and Logarithmic Series Distributions

Author(s): F. J. Anscombe
Source: Biometrika, Vol. 37, No. 3/4 (Dec., 1950), pp. 358-382
Published by: Biometrika Trust
Stable URL: http://www.jstor.org/stable/2332388 .
Accessed: 19/06/2014 19:56

Your use of the JSTOR archive indicates your acceptance of the Terms & Conditions of Use, available at .
http://www.jstor.org/page/info/about/policies/terms.jsp

.
JSTOR is a not-for-profit service that helps scholars, researchers, and students discover, use, and build upon a wide range of
content in a trusted digital archive. We use information technology and tools to increase productivity and facilitate new forms
of scholarship. For more information about JSTOR, please contact support@jstor.org.

Biometrika Trust is collaborating with JSTOR to digitize, preserve and extend access to Biometrika.

http://www.jstor.org

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 [ 358 ]
SAMPLING THEORY OF THE NEGATIVE BINOMIAL AND
LOGARITHMIC SERIES DISTRIBUTIONS
BY F. J. ANSCOMBE, Statistical
Laboratory.
University
ofCambridge
1. INTRODUCTION

The negativebinomialdistribution
dependson two parameters,
whichformanypurposes
maybe convenientlytakenas themeanm and theexponentk.The chanceofobservingany
non-negative r
integer is In'k
r
+(m ) ( 1)

Sometimesit is moreconvenient
to replacem byp or X definedby
m-
m)
_= p1=2m
-l+p m+k' (1.2)

Thus we may write r= (l-X)kF(k+r)Xr(

We assumek,m,p > 0, 0 < X < 1. The factorial-cumulant-generating

function
is

InE{(1+ t)r} =-kIn

E K[ilt1Ii! (1-pt), (1.4)
1=1
and the ith factorialcumulantis Ail = (i-i1)! kpi. (1.5)
The generating
function
ofordinarycumulants*
is

InE(er) _E Kit/i! k1n{1 p(et - 1)}, (1.6)

i=1

and thefirstfourare K= kp =m,

K2= kp(l +p) = k,
M+m12/ (1.7)
K3 =kp(l +p) (l +2p),
K4 = kp(l +p9) (1 +6p+6p2).j
From(1.4) or (1.6) we see thatthesumofN independent observations fromthedistribution
has stilla distribution
ofnegativebinomialform,withmeanNm and exponentNk.
The logarithmic ofR. A. Fisheris obtainedby a limiting
seriesdistribution processfronm
the negativebinomialdistribution a
by considering sampleof N readings,lettingN tend
to infinityand kto zero,and neglectingthe zero readings.It is a multivariate
distribution,
consisting ofa setofindependent Poissondistributionswithmeanvalues
cLX, xaX2, IaX3 (1.8)
A 'sample' comprisesone readingfromeach Poisson distribution.
* For a discussionof ordinaryand factorialcumulantsof a related distributionsee Wisliart(1947).
Aitken (1939) and Haldane (1949) have pointed out that discrete distributionsare often more con-
venientlydescribedby factorialthan ordinarycumulants,and this proves to be so forthe distributions
consideredin ? 2. The followingrelationsmay be noted:
K1 = Kill,

K2 = K121 + Ktlls

K3 = K131+ 3K121+ KellI

K4 = K141+ 6K[31+ 7K21 + Kill.

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 359
The mainpurposeof thispaperis to carrysomewhat further Fisher'sinvestigations
(Fisher,1941;Fisher,Corbet& Williams, 1943)intothesampling propertiesofthesedis-
tributions.Thefollowing is a brief
summary ofcontents.
In ?2 thenegativebinomial formof distribution is compared withsevenothertwo-
parameter forms ofdistribution thathavebeenproposed byvarious It isshown
writers. that
theycanbe arranged inorderofincreasing skewness andtaillength, andthattheyvaryin
thenumber ofmodespossiblein thefrequency function. ThuswhileNeyman'sTypeA
contagious distribution mayhavean unlimited number ofmodes,a distribution givenby
Polyamayhaveeither oneortwomodes, andthenegative binomial anda discreteform of
thelognormal distribution havealwaysonemode.The estimation ofthedistribution of
localmeanvaluesinheterogeneous Poissonsampling is considered.
In ?3 theestimation oftheparameters ofa negative binomial distributionfrom a single
largesampleis considered.* Alternativestothemaximum-likelihood method aredescribed
andtheirefficienciesindicated. Threesuchmethods arefound tobeofpractical importance:
estimation bythefirst twosamplemoments, estimation bythefirst samplemoment andthe
observed proportion ofzeroreadings, andestimation withtheaidofa transformation ofthe
observations whichmakesthevariance independent ofthemean.
In ?4 twolarge-sample testsaredescribed fordiscriminating between alternative forms
ofparentdistribution. Eachtestis fully efficient
incertain circumstances.
In ?5 theestimation ofa common exponent froma seriesofsamplesis considered, when
theparentpopulations possibly differ
intheirmeans.Theresults of??3 and5 havealready
beensummarized bymeelsewhere (1949),andtheirusediscussed.
Finally,??6 and 7 dealwiththelogarithmic seriesdistribution.Theestimation ofarby
maximum likelihood, and somealternative formulae forits sampling variance, are dis-
cussed.Two testsofdeparture fromthelogarithmic seriesformofdistribution are con-
sidered,oneofthembeingduetoFisherandtheothernew.
Notation. Thefollowing notation willbe usedforthenegative binomialdistribution:
Mti, p, X, PI as defined above.
N = totalnumber ofobservations insample.
nr= number ofobservations equaltor (forr> 0).
00
r= E nrr!N= meanofsample.
r=0
00

s2 E nr(r- i)2/(N )-variance estimate.

r=o

Forotherdistributions
in?2,kandp aredefinedso thatthemean= lp, variance kp(1+ 1)).
Thenotation forthelogarithmic
seriesdistribution
willbe:
a., p, X, n,.as definedabove.
00

S= nr.
r=1
00

I-E nrr.
r~1
* The investigation may be compared with that of Shenton
(1949) for Neyman's two-parameter
Tvpe A distribution.

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
360 Negative binomial and logarithmicseries distributions
An estimate of a parameterwill be denoted by the same symbol with circumflexadded.
Differentestimates of the same parameterare not distinguishedin the notation, but only
by context.

2. COMPARISON OF NEGATIVE BINOMIAL WITH OTHER DISO1TR-tiBLT[ONS

A numberof ways are knownin whichthe negative binomialdistributionca&1arise:

(1) Inverse binomial sampling. If a proportion0 of individuals in a population possess
a certaincharacter,the numberof observationsin excess of k that must be taken to obtain
just k individuals with the characterhas a negative binomial distributionwith exponent k
(Yule, 1910; Haldane, 1945).*
(2) HeterogeneousPoisson sampling. If the mean A of a Poisson distributionvaries ran-
domly fromoccasion to occasion, a 'compound Poisson distribution'results (Feller, 1943).
We obtain a negative binomialwithexponentk ifA has a Type III distribution,proportional
to a x2 distributionwith 2k degreesof freedom(Greenwood & Yule, 1920).
(3) Randomlydistributedcolonies. If colonies or groups of individuals are distributed
randomlyover an area (or in time) so that the number of colonies observed in samples of
fixedarea (or duration)has a Poisson distribution,we obtain a negativebinomialdistribution
forthe total countifthe numbersofindividualsin the coloniesare distributedindependently
in a logarithmicdistribution(Lfiders,1934; Quenouille, 1949).
(4) Immigration-birth-death process. A certain simple model of population growth,in
which there are constant rates of birth and death per individual and a constant rate of
immigration,leads to a negative binomialdistributionforthe population size (McKendrick,
1914; Kendall, 1949). The model has been applied to the growthof some livingpopulations,
e.g. populations of bacteria, and to the spread of an infectiousdisease in a community.
The firstof these,inverse binomial sampling,is the simplest mathematically,and is the
onlyone wherethe mathematicalmodel is likelyto hold exactly in practice. While 0 may be
unknownand requireestimation,k is known,and the estimationproblemsdiscussed in the
present paper are irrelevant. Inverse binomial sampling will thereforenot be considered
further.
In general heterogeneousPoisson sampling A may be supposed to have a distribution
dU(A) with mean m. If A has a cumulant-generating function
CIO

s()In E eilA= KEAtilP,! (2 1)

> -1

functionof the distributionof the observed

then q5(t)is the factorial-cumulant-generating
count r, and O(et- 1) the generatingfunctionof ordinarycumulants. Hence if A has a Type
III distributionwith cumulant-generatingfunction

0(t) = - k11 (I - mt/k). (2i2)

r has the negative binomial distributionrequired.
With the model of randomlydistributedcolonies,let the numberof colonies observedper
and let the numberof individuals p per
sample have a Poisson distributionwith mean inm,
* It appears fromI. Todhunter'sHistorythat the earliestgeneralstatementof tHenegativebinomial
distributionin this connexionwas by Montmortin 1714. As to the other methodsof derivingthe dis-
tribution,a moredetailed reviewhas been givenby Irwin (1941).

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 361
colonyhave a distributionwithfrequencyfunctionUP. If the latterdistributionhas factorial-
cumulant-generatingfunction

~f(t) = In E( 1 + t)P= i=1

Liti/ii!, (2.3)

then ml(eAt- 1) is the factorial-cumulant-generating functionof the distributionof the

total count r, and the firstfourfactorialcumulantsof r are
MAL1, ml(L2 +L2), + 3L2L1+ L3),
nln(L3 ml(L4 + 4L3L1+ 3L + 6L2L2 + L4).
WVe obtain the negative binomialwithparametersp and k ifml = kIn (1 + p) and ifp has the
logarithmicdistribution 1 P
P pln(l+p) (l+p) P 24
whichhas mean m2 = p/ln(1 + p) and factorial-cumulant-generating
function
I
00t = In t1- (1 +-pt) (:25)

Models (3) and (4) forthe negative binomial are closely associated, in that we may use
(4) to justify(3). But it will be convenientbelow to considermodels ofrandomlydistributed
colonies withoutspecifyingan evolutionarystochasticprocessthat could give rise to it, and
so the two models have been separated.
While we may expect that distributionsclosely resemblingthe negative binomial will
oftenin factbe observedin population countsand in the samplingofheterogeneousmaterial,
it will not be surprisingif sometimesthe specificassumptionsmade in the above models are
so wide of the mark that a substantially differentform of distributionappears. Before
embarkingon a detailed study of the samplingpropertiesof the negative binomialit will be
as well to considerbrieflywhat otherdistributionshave been proposedthat mightperhapsfit,
such observationsbetter.Attentionwill be confinedto distributionshaving onlytwo adjust-
able parameters; there seem to be seven of these outstandingin addition to the negative
binomial. A convenientmethod of comparisonis to express each distributionin terms of
parametersk andp such that the mean and variance are kp and kp(1 + p), and thenevaluate
the thirdand fourthfactorialcumulants. Results are shown in Table 1. The distributions
can also be compared by computingspecimen frequencyfunctions;this is done in Tables
2 and 3.
The two-parametercontagious distributionof Type A of Neyman (1939) arises fromtile
model of randomly distributedcolonies in which the number of colonies per sample has
a Poisson distributionwithmean ml,ifwe assume that the numberofindividualsper colony
also has a Poisson distribution,say with mean Mn2 mnland m.2 being positive constants.
A derivation along these lines has been given by Cernuschi & Castagnetto (1946), who,
however,appear not to have recognizedwhat they derived. The distributioncan also arise
fromheterogeneousPoisson sampling,if A has a discretedistributionand is equal to in2x,
wherex has a Poisson distributionwithmean mnl.This is moreor less the model Neyman used
in derivingthe distribution.The frequencyfunctionis, forr > 0,
_ C
Pr= 2 e-ml (2.6)
J(erM2)Y'

and its factorial-cumulant-generating

functionis

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
362 Negativebinomialand logarithmic
seriesdistributions
To expressthe mean and variancein the formkp and kp(1? p) we must set ml = k,
M2 = P.
Neyman'stwo-parametercontagiousdistributions
ofTypesB and C werederivedfrom
a morecomplicated
model.The factorial-cumulant-generating
function
oftheTypeB is
(ern2t~1
l (2.8)

and we findIm, = k, 3m2 -pp. The factorial-cumulant-generating

functionof the Type C
distribution
is em2t-1-r t (29)
mlj (m t)2

and 2m1 = k,im2=p. Thesetwogenerating functions,

(2.8) and (2.9),can be derivedfrom
thatofTypeA, (2.7),by a suitableintegration.Thusto get (2.8) we replacem2in (2.7) by
x and ml by (ml/m2) dx,and integrateforx between0 and M2; whileto get (2 9) we do the
same except that ml is replaced by [2in1(i2- x)/m2]dx.The observedvariable r can therefore
be regarded,in each case, as the limitof the sum of a largenumberof randomvariables
followingindependent TypeA distributions withvaluesofthesecondparameter distributed
in a range (0, M2).
Recently,Thomas(1949)has proposeda distributionverysimilarto Neyman'sTypeA.
Withthe modelof randomlydistributedcolonies,wherethe numberof coloniesin the
samplehas a Poissondistributionwithmeanml,the numberofindividualsper colonyis
assumedto be oneplus an observation withmeanM2-1. m2iS
froma Poissondistribution
nowa constant> 1. The factorial-cumulant-generating
functionofthedistribution
is
ml(I1+ t) e(M2-l)t- I}. (2m
10)
If M2 is largethe distribution
is closeto the NeymanType A withthe same valuesofthe
parameters ml, m2. For M2-1 small,we mayset
m2= 1+ jp +jip2 +O(p4), mlm2 =kp.
Polya (1930) gives a distribution
that arisesfromthe modelof randomlydistributed
colonieswhenthe numberof individualsper colonyhas a geometricdistribution with
frequencyfunction (1-T) P-i. (2d1)

p takes positiveintegervalues, and ris a constant,0< -r< 1. The mean numberofindividuals

per colonyis2 = (1 - r)-l. The frequencyfunctionof the total observed count r per sample
is given by 1.r m1''-'r)(
PO= e-ml, Pr = ernt rr (i- ( l ( ) (2.12)

and the factorial-cumulant-generating

functionis
m1lt (2.13)
1- -r -,rt
We findml/(2r)= k,2r/(1 - r) = p. Polya states that the distributionwas givenby A. Aeppli
in a thesisin 1924. It will accordinglybe referredto here as the Polya-Aeppli distribution.
Preston (1948) has considered a distributionderived fromthe model of heterogeneous
Poisson sampling,whereit is supposed that A has a lognormaldistribution,i.e. that In A
has a normal distribution,say with mean 6 and variance 0-2. The distributionmay be

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 363
convenientlyreferredto as the discrete lognormal distribution.* It suffersfromthe dis-
advantage that its frequencyfunctioninvolves an untabulated integral;t forr > 0,

m-~ ~ )I expF U
(2.14)
mrreldu,
=r
r!V/(27T
Inr) j oo 2 In -r

wherem = exp (6 + jO-2),r = exp (o.2). However, the firstfewfactorialcumulantsare easily

found,forthey are the ordinarycumulants of the distributionof A, and these are obtained
(Finney, 1941) fromthe momentsof A about the origin,
ju = exp (if+ ji2o2) = mirii(i-1.) (2.15)

The mean and variance of A are thereforem and m2(r- 1); and we have

(r- 1)-i = k, m(T- 1) = p.

Anotherratherintractabledistributionderivedfromthe model of heterogeneousPoisson

sampling has been given by Fisher (1931), who supposed that A was distributedlike the
square root of a Type III variable. The frequencyfunctioncan be expressed in terms of
Hh-functions, whichhave been tabulated to some extent.The cumulantsinvolve F-functions.
For the entryin Table 1 the limitp -+0 with kp constanthas been considered.The distribu-
tion will be referredto as Fisher's Hh-distribution.

Table 1

Distribution Kl/[3(kp3) K141/(kp4)

Thomas 4+ N1p+O(p2) 2?O(p)

FisherHh 1 +k-'+ O(p2) O + O(p)
Neyman A 1 1
Neyman B 9 27
Neymani C 6S 8
P61ya-Aeppli 4 3
Negative binomial 2 6
Discrete lognormal 3 + k-1 +
16+ 15k-1+ 6k2k-3

The thirdand fourthfactorialcumulantsof the above distributionsare given in Table 1.

It willbe seen that,apart fromFisher'sHh distribution,theyforma sequence ofdistributions
ofincreasingskewnessand tail length(leptokurtosis)in the ordershown.The positionofthe
Hh distributionrelative to the others is ambiguous and variable (ambiguous in that it
depends on whetherwe rank by K[3]or by K[4, variable because it depends on the values ofthe
parameters),but we may say at least that it should come somewheretowards the frontof
the list.
* Prestondoes not give any exact samplingtheory.Otherwriters(e.g. Williams,1937; Gaddum,1945)
who have alluded to lognormaldistributionsin connexionwithfrequencycountshave contentedthem-
selves withrecommendingthat the data should be transformed by a logarithmictransformation of the
formy = log (r + c), so as to appear approximatelynormal. It should also be noted that Prestonis con-
cernedwith the situationwherezero counts are not recordedand thereforethe total sample size N is
unknown.This will be discussedin ??6 and 7.
t A usable approximationto the frequencyfunctionof the discretelognormaldistributionhas been
developed by Dr P. M. Grundy.

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
364 Negativebinomialand logarithmic
seriesdistributions
The difference is clearlysubstantialifp is large.To
in shape betweenthe distributions
demonstratethisfurthersomeexpectedfrequencies are givenin Table 2 forthreedistribu-
tions havingp = 10, mean = 20, variance = 220, namely,
(a) NeymanTypeA, withml = 2, m2= 10;
withml = 31, X = a;
(b) P6lya-Aeppli,
(c) negative binomial,with k = 2, m = 20.
havingmean= 20,variance= 218,namely,
Alsoshownis a distribution
(d) Thomas,withml = 2, M2 = 10-
To save space,thefrequencies have beengrouped.The Neymandistribution (a) has modes
orpeaksat r = 0, 10,20,whileat r = 30 thereis a modein firstdifferences
whichis notlarge
enoughto producea modein the frequencies themselves.The Thomasdistribution (d) is
practicallyindistinguishable
from(a), thedifferencebeingthatthemodesof(d) are slightly
morepronounced thanthoseof (a). The Polya-Aepplidistribution (b) has twomodesonly,
at r = 0 and 11.The negativebinomial(c) has onemode,at r = 9 and 10 (equal frequencies).
If a discretelognormal wereaddedto Table 2 (withm = 20,T = l-5,6 = 2 7930,
distribution
and oC= 0.6368)it wouldresemble(c) in havingonlyone mode,but wouldbe rathermore
skew;thefrequencies forthefirstfewvaluesofr wouldbe lower.
Table 2

Percentagefrequency Percentagefrequency
r , , , i r
(a) (b) (c) (d) (a) (b) (c) (d)

0 13-53 3.57 0.83 13-53 17-18 5 04 5.52 5.77 5.09

1- 2 007 4*18 3.55 0)03 19-20 5.37 5.18 5.28 5.52
3- 4 0.72 4-96 5.31 0 54 21-22 5.23 4.81 4.79 5.33
5- 6 2.74 5.52 6.34 2.56 23-24 4.74 4-42 4.31 4.74
7- 8 5,54 5.89 6.86 5.66 25-26 4*19 4-02 3.86 4.15
9-10 7.01 6.07 7.01 7.29 27-28 3.78 3-63 3-43 3.77
11-12 6-41 6.10 6.90 6-47 29-30 3.48 3.25 3 03 3651
13-14 5.18 6.00 6.62 5.00 31-32 3.19 2.89 2.67 3.23
15-16 4.73 5'80 6.22 4.59 33-oo 19.05 18-20 1722 18-98

A lessextremecomparison
is showninTable 3,whichcomparesdistributions
havingp = 3,
mean = 6, variance = 24, namely,
(a) NeymanTypeA, withml =2, m2= 3;
(b) Polya-Aeppli,
withml = 2-4,X = 0-6;
(c) negative binomial,with k = 2, m = 6.
(a) and (b) have twomodes,(c) has one.
In general,theNeymanTypeA and Thomasdistributions can haveanynumberofmodes
fromoneupwards,andifthereareseveralmodestheywilloccurat valuesofr approximately
equal to multiplesofin2. The Polya-Aepplidistribution,
on the otherhand,has eitherone
or two modes-two if 2 < ml< (1 T)-, one otherwise.
- The negativebinomialdistribution
has alwaysone mode. Presumably,by analogy,Fisher'sHh distribution has one or two
modes,and thediscretelognormal one mode.

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 365
For the merepurposeofgraduatingdata thereis littleto choose betweenthe distributions
in shape ifp ( = m/ik)
is not large, and then considerationsof ease of handlingweighheavily
in favour of the negative binomial. Experimental discriminationbetween the different
formsof distributionis practicable, however,if p is large. An interestingattempt at such
discriminationhas been made by Beall (1940), who fittedNeyman Types A, B, C, Polya-
Aeppli, and negative binomial distributionsto eleven seriesof counts ofinsect larvae. Some
of the seriesseemed to indicate a bimodal population, and he concluded that they were well
fittedby the Neyman formsof distribution,but not by the othertwo. Mr D. A. Evans has
pointed out to me, however,that Beall fittedthe latter distributionsincorrectly,having
mixed up the two parameters,and he was consequentlyunfairto them.

Table 3

Percentage frequency Percentage frequency

(a) (b) (c) (a) (b) (c)

0 14.95 9*07 6*25 9 5 23 4 92 4 69

1 4.47 8*71 9*37 10 4*42 4*12 3.87
2 7.36 9 41 1055 11 3.67 3 40 3*17
3 8 77 9 49 1055 12 3 00 2*78 2*57
4 8-83 9*12 9 89 13 2 42 2 25 2*08
5 8.29 8.46 8*90 14 1*92 1 80 1*67
6 7*60 7.63 7.79 15 1.51 1 43 1 34
7 6.86 6*71 6 67 16 1*17 1 13 106
8 6 06 5'80 5 63 17-co 3 48 378 3 95

In analysingpopulation countswe may have two quite distinctobjects. On the one hand,
the counts may have been made on plots in an experiment,and we desire some means of
interpretingthemso that the effectsoftreatmentscan be judged. What is usually done is to
apply a transformationto make the method of analysis of variance appropriate. We study
the distributionof the originalcountsin orderto finda suitable transformation.It does not
mattergreatlywhetherthe formof distributionfitted,if any, is very accurate or not.
On the otherhand, we may be interestedin relatingthe observed counts to some theory
of population growthor spread. In that case we shall endeavour to use only formsof dis-
tributionthat are 'biologicallysignificant'.Neyman's distributionswere intendedto repre-
sent populations ofinsect larvae observed shortlyafteremergencefromeggs,the eggs being
supposed laid in clustersof a fixedsize. The models seem ratherspecial, and not likelyto be
widely applicable to other sorts of population counts. The characteristicfeature of the
distributions(at any rate, that of Type A) is the possibilityof a series of three or more
equally spaced modes. Unless such a seriesof modes is demonstratedconclusivelyby obser-
vation, one may reasonably feel reluctantto use such a model. Thomas's distributionwas
intended to represent plant quadrat counts, but no evidence has been adduced to make
plausible the special formof up assumed, and again one may reasonably feel reluctant to
use it. The derivationofthe Polya-Aepplidistributionfromthemodelofrandomlydistributed
colonies is much more promising. Kendall (1949) has shown that the progenyof a single
individual aftera fixedlapse of time will followa geometricdistributionwith modifiedzero

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
366 seriesdistributions
Negativebinomialand logarithmic
term,in certainfairlygeneralconditions ofno competition. Henceifprogenitors (e.g.plant
seeds)are releasedrandomlyoveran area at onetimeand theirprogeny(freely increasing
by vegetativereproduction) are observedat a latertime,we shall expectthe numberof
individualsper quadratto followthe Polya-Aepplidistribution. The parameterml willbe
themeannumberofprogenitors perquadratof whichsomeprogenysurvive.If,insteadof
beingreleasedall at onetime,theprogenitors arereleasedwithuniform distributionin time
froma particulartimeup to thepresent,and ifthebirth-and death-rates per individualare
constant, wegetthenegativebinomial,as alreadyremarked at thebeginningofthissection.
We maytherefore expectthatcloseapproximations to boththesedistributionswillin fact
be observedin thestudyofgrowing populations.Ofcourse,somepopulationcountswillnot
resembleany ofthedistributions we have considered, on accountofovercrowding or,with
mobilefauna,aggregating forreproduction, defence,or othersocialpurposes.It is unlikely
thatanytwo-parameter distributionwilldescribesuchcountsadequately.
In viewofthe difficulty ofdiscriminating experimentally betweenformsofdistribution
arisingfromdifferent mathematical models,thestudyofcountsmadeall at onetimeis not
likelyto givereliableinformation on lawsofpopulationgrowth.For thispurpose,repeated
observations onthesamepopulationareneeded,ifpossiblewithidentification ofindividuals.
Sometimes, in samplinginvestigations, it is reasonable to supposethat the observations
arisefromheterogeneous Poissonsampling,but theremay be no definite groundsforpre-
dictingthedistribution ofA. If theobservations numerous,
are sufficiently thedistribution
ofA can be estimated(Newbold,1927). Let Al,k2,..., be the k-statistics calculatedfroma
sampleofN observations onr (see,forexample,Kendall,1943).Thenwecantakeas unbiased
estimatesofthefirstfourcumulants K1,K2,K3,K4,ofAthefollowing:
A

Ki = kl,
AI
K2 = k2-1kl,
A (2.16)
K3- k3- 3k2+ 2kl,
AI
K4 = k4-6k3+1lk2--6kl.
The right-handsides are in factunbiasedestimatesof the factorialcumulantsof the dis-
ofr, analogousto thewell-known
tribution forordinarycumulants.It is quite
k-statistics
to calculate
straightforward thevariances (and other cumulants)oftheseestimates.If,for
we find(forN large)
example,K1 = m,K2 = cr2, and all theKi fori > 3 arezeroornegligible,

var(K1) = N

var(K2) - 2(m+of2)2+ 20-2 (2.17)

var(K) 6 (m+ 0-2)3 + 18oT2(M + 302)

theprocesscurveoftheoutputof
One possibleapplicationofthismethodis to estimating
a productionline frominspectionrecordswhichgive the numberof detectivesfoundin
samplesof a fixedsize takenfromeach lot produced.A similarbiologicalproblemin sur-
veyinga districtforpresenceof potato-rooteelwormhas been describedby me elsewhere
(1950). Newbolddevelopedthemethodin a studyofaccident-proneness.

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 367
3. FITTING A NEGATIVE BINOMIAL IYSTRIBUTION TO A LARGE SAMPLE

Haldane (1941) and Fisher (1941) have considered the maximum-likelihoodmethod of

fittinga negative binomial distributionto a large sample. If the distributionis expressedin
termsof the parametersm and k, the maximum-likelihoodestimatesof the parametersturn
out to be independent(asymptotically). For the estimate Anof m we have simply
Mn=r. (3 1)
The estimate k of k is the root of the followingequation in x:
I
Nn (I +/x)= = Xn +I + x+) (3.2)
j=1 \XX+lX +j -l
It is easy to showthat theleft-handside is greaterthan theright-handside ifx is largeenough,
provided (N -1) S2> NF (or, ignoringthe difference betweenN and N - 1, if 82 > F); and also
that the left-handside is less than the right-handside if x is small enough (but positive),
provided that no< N, so that there are some non-zero observations. Since both sides are
continuousfunctionsforx > 0, the equation musthave at least one finitepositiveroot. I have
not proved that thereis only one root in this case, and that if S2< F thereis none, but after
an unsuccessfulsearch for a gegenbeispielI suppose both these statementsto be true. If
82< F (and if,indeed, thereis no finiteroot), the excess of the right-handside over the left-
-z oo,and we may say k = oo,implyingthat a Poisson distribution
hand side tendsto zero as x
is being fitted.If no = N, we may conventionallytake k = 1, say.
For the variance of A we have
var() = (n+) (33)

From the matrix of expectations of second derivatives of the log-likelihoodfunction,we

obtain the large-sampleformulae:
cov(S., k) O.0 (3 4)

varkV 1) I( + 4X
+ 2)+ (k+ 2) (k-+ 3)
Na kk+ )( IX 3(kc
_I_
-
___

vak(kk-1X (3.5)
The summationinvolved in derivingthe second of these is due to Fisher; the seriesin curly
brackets may be written

1 += j + I (lo+ 29)(k+ 3) ... (k )

It may be noted in passing that large-samplevariances and covariances foundin this way
relate to the asymptoticnormaldistributionofthe estimates,and are not necessarilyasymp-
toticallyequal to the actual variances and covariances of the estimatesforfiniteN. In the
presentcase, in samplingfromany negative binomial distribution,thereis always a positive
chance, albeit perhaps a minuteone, of finding82 < i, when we should set k = oo. kidoes not
thereforehave a proper distribution,nor variance. It does, however,have an asymptotic
distributionwiththe variance given, as N-a-oo.
While equation (3.1) gives Anverysimply,equation (3 2) fork is tedious to solve, and it is
worthwhileto consideralternativemethodsofestimatingk. We beginwitha generalmoment
method. Let fr be any specifiedconvex or concave function(not linear) of the non-negative

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
368 and logarithmicseries distributions
Negative bin~iomial
00
integerr.Thenwe mayconsiderE nrgf
as a statisticfork. Let E (fr), expressedas a function
r= A
ofmandk,bedenotedbyF(m,k).Thenweshalltakeas ourestimatektherootofthefollowing
equationin x: 1
yrnrfr=F(r, x). (3.7)

Theright-hand sideofthisequationis a monotonefunction ofx,withno constantstretches,

x
for positive,providedr > 0. For if x is increasedto x + 6x,the change MIrin Pr(m, x) is
negativefor0 < r< a and forr> b+ 1, and non-negativefora + 1 < r< b, being strictly
positiveforone ofthesevaluesofr at least,wherea and b are integerssatisfying 0 < a < b.
0 o00 00

Moreover,E &PJ.E r8r = 0. It followseasilythatiffr is convexfrombelowE fr6I, < 0,

r=() r=O r=O
whileiffris concavefrombelowtheinequalityis reversed. Thusineithercase theright-hand
side of (3.7) is monotone,as stated; and therefore the equationhas at mostone positive
solutionk. On the otherhand,in repeatedsamplesfromthe same negativebinomialdis-
tribution, theprobability is smallthat (3.7) has no positivesolution,whenN is large. For
theleft-hand sidehas a highprobability ofbeingcloseto F(m,k),and F(i',x) is a continuous
function ofI-and x witha rangeofvalues,as x variesin (0,so),differing littlefromtherange
ofvaluesofF(m,x) ifr is nearto m.
We can findthelarge-sample varianceofthe estimatek givenby (3.7) by treatingk- k
and r- m as infinitesimals. Denotingthe latterby 8k and 6mrespectively, and nr/N- Pr
by &LP,we have co
An = E r8P,
r=O
00

Ak~k = X (frr-Amr) 6J,

-=O

a ~~~~a
where A= --F(m,k),
am Ak- -F(mk-).
ik
Now var(4PT)= P( 1-Pj)/N, cov(Pi,,Pj)-- PjPj/N (i,j > 0, i *j),

whilefrom(1 1) we find (m+--) Ain= E (rfr)-mE (fr).

Hencewe obtainthedesiredresults:
cov(r,k) - 0, (3.8)

var
(k) r~ [E(fr)]2-(m+m2/k)Am (3 9)

The ratioof the variances(3.5) and (3.9) is the large-sample efficiencyof method(3 a7)
ofestimating k. As alreadyremarked, mis easilyestimatedwithfullefficiency by equation
(3s1). Sincetheseestimatesare in largesamplesindependent, it is appropriateto consider
separately.In general,whenconsidering
theirefficiencies inefficientestimatesofthepara-
metersofa distribution, a reasonablesingle measure oflarge-sample is provided
efficiency
bythesquarerootoftheratioofdeterminants ofsamplingvariancematricesforthemaxi-
inum-likelihood estimatesand forthealternative estimatesunderconsideration, ifthereare
two parametersto be estimated,or the cube rootif thereare threeparameters, etc. This

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 369
measureis theratioofthenumbers ofobservations
required
bythemaximum-likelihood
methodand thealternativemethodto achievethesameerrorvariancedeterminant.
Tn
thepresentcaseitwouldbeequaltothesquarerootoftheefficiency
ofestimation
ofk.
Letusnowconsider someexamples.
Method 1. ftr=,
= so thatkis estimated
fromthesamplevariance
$2, i.e.
^ 2
k =2 -a (3.10)
On evaluatingE (fr),E (fl), Ak,Am,we findthelarge-sample
variance
vAr 2k(k+ 1).
var(k) NX2 X (3.11)
andtheefficiency
ofestimation ofkis equaltothereciprocal oftheexpression
(3.6).
Method2. fo= 1,f,= 0 forr > 1, so thatk is estimatedfromtheobservedproportion of
zeros,i.e. nO/N = (1+r/k)-. (3.12)
Thelarge-samplevarianceis

var(k)-N[ -Il (X)-_X]2. (313)

Method3. fr I/(r+1). k is given by
1 ~r where X=
A

(1.~X)~(~XA (3.14)
Nror+ (k-1) Xk
Thelarge-sample
variance
ofkis givenby(3 9),where
E (fr)=( (k)-1()X ), E (fr) =1X)k E (r(+ +)!
1)j(r ) Xrl

Am
k k(k - ) X2[{
(k-1)2X
+ (k-1I 1Xk(-) (3 15)
I
A -l2X{ k(k-1)In (I -X) -(k -1)2X) (I -X)k -(I1 X)] .

This is fork * 1. Correspondingexpressionswhen k = 1 are easily founddirectly.

Method4. f = Cr, wherec is a positive constantnot equal to 1. We find

j = [1 + (1-c)
2flnrCr /k]-k; (3.16)
r=O
and forthe large-samplevariance
E(fr) = [L+p(jC)]-k, E(fr2) = [j+p(lc2)]-k,
Am = (1-c) [1 +p(1-c)]-k-1, (3I17)
Ak = [1 +p(c)-k( ln [ +p(l c)] +lP (1(C)) (

The above seem to exhaust the tractable applications of the moment method. Several
other formsforthe functionfrsuggestthemselvesas possibly worthinvestigating,such as

(i)fo=O,fr=+I+f+...-forr1,(ii)fr=ln(r+1)forr O,and(iii)fr=f/rforr O;
2 3 r
but unfortunatelythese do not lead to a simple expressionforE (fr). Their practical use is
thereforeruled out unless special tables are constructedforestimatingk fromEnrfr/N.
Biometria37 24

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
370 Negative binomial and logarithmicseries distributions
A methodof a different kind is the following.
Method5. We guess a value of k, apply a transformation(dependingon k) to the observa-
tions to make the variance independentof the mean, and then obtain an improvedestimate
of k by equating the observed variance to the expected; the process is repeated until the
answer becomes stable. Suitable transformations were investigatedby me (1948) and sum-
marized by me (1949). From a considerationof equation (4.37) of the formerpaper, which
gives an asymptoticexpansion forthe variance of the transformedvariable when k = 1, it
appears that the transformationmethod is unlikelyto be serviceable fork< 1. For higher
values of k,the methodis possible fornot-too-lowvalues ofmi. In principlethe methodcould
be used forany values of m and k,if the expected variance of the transformedobservations
were knownas a functionof m and k in a convenientform,e.g. by an adequate double-entry
table. But such a table is not available, and so the methodis restrictedto those values of mt,
forwhichthe limitingvariance as m--oo is near enough attained.
The large-samplevariance of the estimateof k derivedby the transformation methodwill
now be found. Let y denote the transformedvariable, a functionof the observed count r,
whenthe truevalue ofkis insertedin the formulafory; and let 9 be the transformedvariable
when an estimate k of kis used. Then in samples froma fixednegative binomialdistribution
with parametersm and k,we have, fromequation (4.23) of my 1948 paper (settingA = Ok),

var(y) =V/(lc)+(kI+_ +o(-2i), (3.18)

as mo-* co, if k > 2.* Let S2 be the variance estimate calculated fromthe observationsafter
transformation to y-variables,and s that calculated afterthe observationshave been trans-
formedto 9-variables. We choose k, by successive approximation,so that
S = r'(k). (3.19)
Let as2/akdenote the derivative of 82 (forthe given sample of observations)with respectto
the parameterk enteringinto the y-transformation, when k is set equal to its truevalue. For
any arbitraryk, not necessarilyclose to k, we have in probabilityas N -? oo
8 = var(9)+ 0(1/N),

and therefore * = [-2var(9)] +o( ).

Hence, if k is determinedby (3 19) and we set k - k = 8k and - var (y) = &O-2,

we have in
probabilityforlarge N and m
2
0'(k) + #(k) 6k = i/r'(k)
-S2 = 82 + 82-(k=
alc ?'(k) + '
2rn,'
~~~~~(kc-1)

and therefore 8f2= (Vr"(k)+(Ik) m}8k.

The large-samplevariance of the left-handside is foundfromequations (4.23) and (4.30) of

my 1948 paper, and we obtain the result
..(k)++ 2[?*(k)]2(32)
var(ick(k) /"() 2/']2'k (3.20)
N[Pr"(k) + (k1-1)2 m]
* *-(t),*'(t), ... denotethesuccessive
derivatives
ofInr(t).

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 371
This should be sufficiently accurate forpractical purposes if m is above 50 and k above 5,
assumingthat the appropriateinversehyperbolicsine transformation is used. It is not clear
fromexistingcalculations how far (3.20) may be relied on outside this region.
We thus have an assortment of alternatives to the maximum-likelihoodmethod of
estimatingk. Let E1, E2, E3, E4, E5 denote the large-sampleefficiencies of the fivemethods.
By tabulating these efficiencieswe can see under what conditions,if any, each method is
likely to prove useful. In the figureare shown 50, 75, 90 and 98 % contoursof E1 and E2,
and 90 and 98 % contoursofE5. The 75 % contourofE5 has been foundonlyform verylarge,
whenit is close to the line k = 1; it is not shownin the figure.Since Method 5 can hardlybe
used when k < 1, no attempthas been made to determinethe 50 % contour.

40 _4__

'l0

2 04 09-8204 1 2 4 102040 100 200400

Mean m
Fig. 1. Large-sample
efficiencies
ofestimation
ofIc.
Method1 Method2---- Method5-.- -
th
90%
In consideringthe lower right-handregionof the figure,where m is large and icis small,
it is helpfulto note* that if we set a =--lln ( 1-X) and let lc-., X- 1, witha constant,
the expression (3.6) is asymptoticallyequal to 2(1 - ea)/k, whence the limitingefficiencies
of Methods 2, 3 and 4 are all equal to a2 ea/(ea -1)2. Otherlimitingformsof the efficiencies,
formovementaway fromthe centreof the figurein various directions,are easily foundand
need not be quoted.
No contoursare shown in the figureforMethods 3 and 4, since, as it turnsout, these are
nowheremore efficientthan the more efficientof Methods 1 and 2, and the latter are more
convenientto use. Method 4 becomes equivalent to Method 2 as co-+ , and to Method 1
as co-1. If we imagine c increasingcontinuouslyfrom0 to 1, a constant-valuecontour of
te9
K4 such as /contour, departs fromthe correspondingcontourofK2 by a translation
of the uppermostpart of the curve to the right(' east ') and a pullingof the middle part of
the curve downwards to the left ('south-west'). (The firstof these movements is easily
expressed.The contourhas a verticalasymptotefork large,of the formmn= constant.This
value of mnis equal to the value forthe correspondingcontourof K2 multipliedby ( 1- c)-2.)
As c increases,the contoureventuallybreaks up intotwo disjointparts,and two new parallel
* Proofsofthis and otherstatementsin the remainderof thissectionare omitted,to save space.
24-2

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
372 Negativebinomialand logarithmic
seriesdistributions
asymptotesappear,of the form n/k= constant.As c-r-1, the upperpart of the curve
approachestheE1 contour,whilethelowerpartrecedesand vanishesin thelimit.The E3
contoursare similarin characterto E4 contours.

4. TESTS FOR DEPARTURE FROM THE NEGATIVE BINOMIAL FORM OF DISTRIBUTION

We have consideredhowa negativebinomialdistribution can be fittedto a largesample.

We turnnow to testinggoodnessof fit,in particularto detectinga departurefromthe
negativebinomialformtowardsone oftheotherformsdiscussedin ? 2. Testsare required
thatarereasonablyconvenient to use.
Particularinterestattachesto discriminating
betweenthe negativebinomialand the
Polya-Aeppliformsof distribution.Let us supposeto beginwiththat m/k(=p) is small.
Thenthelog-likelihoodfunctionoftheobservations
on thehypothesisofa negativebinomial
distribution
is
00 00 00 00
(r -

Xnrln(r!)+ EnrI(")"'P
) 1

nrlnPr= Enrrlnm-Nm-
=
Lo-E
r=O r=O r=2 r=O m 12

- n[r(r r) _3r+ 2m + (p3), (41)

as p -* withmfixed.The maximum-likelihood equationsforestimating k andp areasymp-

toticallythoseofMethod1 ofthelast section.Considernowthehypothesis thattheobser-
vationsaredrawnfroma Polya-Aepplidistribution. In termsofparameters k andp defined
as in ? 2, the log-likelihood
functionofthe observations,
L1 say, is the same as the above
expression forLo exceptthattheterminp2 is

r on[r (r1 -')r+m] (4.2)

Againthemaximum-likelihood methodofestimating k andp is asymptotically equivalent
to fitting
by thefirsttwosamplemoments.The likelihood ratiocriterionfordiscriminating
betweenthe two distributions is foundby maximizing Lo and L1 separatelywithrespect
to k and p, and thensubtracting them. Its leadingterminvolvesthe firstthreesample
moments.We are thusled to proposethe following testfordeparturefromthe negative
binomialformtowardsthePolya-Aeppliform,to be usedwhenp is small:
Test1. Estimatetheparametersofthenegativebinomialdistribution fromthefirsttwo
samplemoments(Method1),and thencomparethethirdsamplemomentwiththeestimated
thirdmomentofthedistribution.
Thereis no pointin actuallyevaluatingthelikelihoodratiocriterion justdescribed,since
wedo notknowa priorithattheparentdistribution is necessarily
ofone or otherofthetwo
formsconsidered. We may,however, applya testsimilarto Test 1 to see whether ornotthe
observationsare in agreementwiththe Polya-Aepplior any otherhypothetical formof
distribution.
We can similarlyinvestigatethelikelihoodratiocriterion in anothersimplelimiting case,
namely,forPo-? 1 withm constant.We find,forbothdistributions, that the maximum-
likelihoodestimationof the parametersinvolvesasymptotically onlythe two statisticsr
and nO/N(Method2 above). The leadingtermofthelikelihoodratiocriterion involves,in
00

additionto these,a further

statistic,E nrInr. This is not a convenient
statisticon which
r=2

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 373
tobasea test,sinceitsexpected valuecannotbeexpressed simply.Theonlysimple testthat
suggests itselfis,infact,onebasedonthesamplevariance, thus:
Test2. Estimatetheparameters ofthenegative binomial distributionfrom thesample
meanand theobserved proportion of zeros(Method2), and thencomparethe sample
variance withtheestimated varianceofthedistribution.
Test2 arisesmoredirectly in another problem ofdiscrimination.We supposethatthe
parent distributionoftheobservations is a compound (orheterogeneous)negative binomial,
i.e. thateachobservation is drawnfrom a negative binomial distributionhavingthesame
exponent kbuta meanthatvariesrandomly from observation toobservation ina distribution
withmeanrnand variancee, say. The resulting distributiondepartsfromthenegative
binomial form towards thatofthediscrete lognormal, havinghighskewness andkurtosis
coefficients.We nowtestthehypothesis thate = 0. If thelog-likelihood functionofthe
observations is expanded in ascending powersofe,thecoefficient ofe is a functionofthe
first twosamplemoments (andoftheunknown parameters mandk).Theoptimum large-
sampletestofthehypothesis thusconsists offitting kandmbymaximum likelihood,assuming
thate = 0,andthencomparing theobserved andestimated variances.In theregion where
Method2 offitting kis efficient,wewoulduseTest2 above.In theregion whereMethod1
is efficient,
i.e. fork large,a smallheterogeneity in themeanofthenegative binomial has
theeffect ofapparently reducing kbutnototherwise changing theshapeofthedistribution,
sothattheheterogeneity isnoteasilydetectable. Test1couldbeusedtodetecta pronounced
degree ofheterogeneity inthemean.
It remains toconsider howTests1and2 arecarried out.Weshallsupposethesamplesize
N tobelarge.Thecriterion ofTest1is thedifference between thethirdsamplemoment and
itsestimated value,i.e. 1 () ( 282
n,Or(r-r)-( i l) (4

Usingtheknown formulae(quotedbyKendall,1943)forthevariances andcovariances of

k-statistics,
or,alternatively,
workinginterms ofthesamplefactorialmoments (ofwhichthe
variances
andcovariances areeasiertofinddirectly),weobtainthelarge-sample result
Nvar(T) 2k(k+ l)p3(l +p)2 [2(3 + 5p)+ 3k(1+p)]. (4.4)
ofTest2 isthedifference
Thecriterion between thesamplevariance anditsestimatedvalue,
i.e. U =: 82_-_2/' (4.5)

wherekA
is theestimate
ofk byMethod2. Wefind,
forN large,
IVcov(k,r) O.0 N cov(k,8S2) -k(k +r 1) p2/{ - In(I 1-X)- XI, (4 6)
whichand(3413)weobtain
from
(I _ X)-k
Nvar(U) 2k(k+ 1) p2(1 +p)2[1- I_ 1X 4 (4.7)

As an example,we mayconsider a frequency quotedby WVilliams

distribution (1944,
Table6) ofthenumberofhead-liceofallstagesfoundinthehairofHindumaleprisoners
on
admissiontoCannanore
jail,SouthIndia,overa period to 1939
fromr1937 (seeBuxton, 1940).*
* Prof.Buxton has verykindlyallowed me to check the frequencydistributionagainst the original
records.Thereseems to be one errorin Williams'stable, the numberofheads withoutlice beinggivenas
622 instead of 612. Williamshas also made a numericalslip in fittinga negative binomialdistribution
to the observations,so that the fitappears worsethan it should.

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
374 seriesdistributions
Negativebinomialand logarithmic
We have N = 1073, no= 612, r = 6-93569, s2 = 583-8. By Method 2, k = 0-144198, and
hence U = 243.3, with estimated standard error38-9. We conclude that the observations
are not adequately fittedby a negative binomial distribution,being too skew. The resultis
not surprising,in view of the heterogeneityof the prisonersin caste and in other respects
likelyto affectpersonal hygiene.

5. FITTING A COMMON EXPONENT TO A SERIES OF SAMPLES

We suppose now that we have a series of samples, one fromeach of io negative binomial
distributions.Characteristicsof the ith distributionand the sample fromit will be denoted
by the usual symbol with suffixi added. We may be interestedin investigatinghow ki is
related to ma. If the kiare not too small we may estimatethem by Method 1 and plot them
against ri. We shall obtain a less skew distributionof errorsif the reciprocalof the estimate
of ki,ratherthan the estimate itself,is considered,i.e. (SA- rA. In either case, there is
a bias of orderN-1 whichmay be worthremovingif v is large. The followingestimateof kj-1,
I _82- __
1 (1 + (5.1)

is approximately unbiased, having expected value kI1 + O(N-2). This is easily proved
(dropping now the suffixi, for convenience) by writingr = kp+ &m,82 = kp(l +p) + &y2,
and expanding (5. 1) in ascending powers of &mand &o.2. To the firstorderwhen N is large,
the variance of thisestimateof ki-1is equal to the right-handside of (3.11) divided by 14,and
its correlationwithr is zero, in samples fromthe same distribution.
Sometimesit is reasonable to suppose that all the kiare equal. We then desirean efficient
estimate of the commonvalue. A method that suggestsitselfis to choose a value of k such
that the sum (or a weightedsum) of82 -_ r/k - vanishes. The expected value ofthisexpres-
sion, when k is replaced by the true ic,is O(N 71), and it can be reduced to O(N -2) if we use
instead ^
t (k) = s8?-(s + i/k) (1- [Nk]-A).
(5A2)

If weightswi are used, our equation fork becomes

= 0.
Z wAtR(k) (5.3)

Treatingk - k as infinitesimaland workingto the firstorderforNilarge, we findeasily

var (k) - A" t )2) (5.4)

N var(ti) , k(k+ 1) pgi( + p)2' (5-5)

stlake (5.6)
var (k) is a minimumwhen the weightswi are taken to be proportionalto
I
aki i.e. to (i

-I
We thereforechoose wi = Ns (5.7)

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 375
The numeratorNi - 1 is more convenientthan N*,since (Ni - 1) s2 is foundin the course of
calculatingsA.Thus the methodis to choose k to satisfyITi(k) = 0, where

()(Ni - 1)A - (Ni - I-k-1) (rij+ r91k) 58

Tj(k) -(?)2(5.8)
(ri + k)
It is easy to verifythat T1(k)has expected value O(Ni-1). The methodis equivalent to taking
an appropriatelyweightedaverage of the estimates (5-1); and, to the firstorderforNi large,
the variance of k is the reciprocalof the sum of reciprocalsof the variances (3.11) foreach
sample. In samples froma singlepopulation,the correlationbetween T(k) and r is O(N-1).
We can deal similarlywith estimates of k based on Method 2. Correspondingto (5.2) we
may consider /\ k
+
A,2 %i( +j i) (N-(k 1) 9
k 2(r + k) I ~~~~~~(5.9)
the expected value of whichis O(Nj-'). The optimumweightfactorwi (forNi large) is

W= ___ (5.10)
- 1 -
ti =(-1 - X&) [(Jl- X)-AC kXi] (* 0
or ratheran estimateofthis. Since thereis littleto be gained by exactitudein weightfactors,
we may preferto take instead
wi= n(1+ s/k), (5.11)
whichhas roughlythe same effectand is easier to calculate. If (5 10) is used, we findthat,
to the firstorderforNi large, the variance of k is the reciprocalof the sum of reciprocalsof
the variances (3 13) foreach sample. To test whetherk changes progressivelywith m, we
may plot Uj = wjui against ij and look fora correlation.In samples froma singlepopulation,
the correlationbetween U and r is O(N-1).*

6. FITTING THE LOGARITHMIC SERIES DISTRIBUTION

Fisher's logarithmicseries distributionwas proposed as a model forthe relative abundance

of differentspecies found in trap catches or other methods of sampling; in particular,to
describethe relationbetweenthe numbersofmothsofdifferent species caughtin a light-trap
over a period of time. Suppose there are N species that mightbe observed, and that their
abundances (numbersexpected to be caught in a unit period of time) are distributedas if
they were a sample of size N froma Type III distribution,proportionalto a %2 with 2k
degreesof freedom.We suppose furtherthat the individuals of each species move indepen-
dently,so that the numberof individuals of any species caught has a Poisson distribution
with mean value equal to the abundance multiplied by the lengthof time of observation.
The numbersofindividuals caught per species thenforma randomsample ofsize N fromthe
negative binomial distribution(1 1). The resultsof the observationcan be expressed by the
numbersni of species representedby i individuals forall i > 0. We defineS and I as in ? 1,
so that S is the numberof species representedby at least one individual, and I is the total

* In my 1949 paper thereis an oversighton thispoint. Having noticeda correlationbetween U, and

ri in some data, I added (p. 172): 'The effectis too markedto be attributedto the negativecorrelation
betweennoand r that occursin repeatedsamplingofthe same population.' This is literallytrue,but mis-
leading,since the relevantcorrelationis that between U and r, not noand r, and that is O(N-1).

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
376 Negativebinomialand logarithmic
seriesdistributions
number of individuals of all species observed. The probabilitydistributionof nj fori > 1,
given N, k and X, is
N! (I -X)kNkn, k(k2{ 1) ((k
(2 + 1) (k + 2)&
xi.
I

(N -S)! fin1i!
i=1

If we set Nk = a, and considerthe limit k-0 , N --a , with a constant,we findeasily that
the above breaks up into a productofPoisson frequencyfunctionsforni (i > 1), as indicated
at (1-8). If the timeofexposureor attractivepowerofthe trap weremultipliedby a factorc,
withoutthe abundances of the species being affected,it followsfromthe above assumptions
that p would be changed to cp, and thereforeX to cp/(cp + 1), while a would be unaltered.
a is thus a propertyofthe biologicalassociationthat is beingexamined,and has been termed
by C. B. Williams the 'index of diversity'of the association.
The log-likelihoodfunctionof the observationsis

L=caln(1-X)+Slna+llnX- 2{rlnr+lnnr!}. (6 2)
r=1
Thus S and I are jointly sufficientfor estimatinga and X, and the maximum-likelihood
equations are
I= &X/(1-X) = Lp } (6.3)
S=I-&ln (1-X) = (I up).

On invertingthe matrixof expectations of second derivativesof L, to findthe variances of

these estimatesin the usual way, we get

1(6-4)
1
[-In (1-X)-X] var(X) ~--X( 1-X)2 In(1-X),l
A, ~~ ~ ~
[-In(I1- X) - X] cov(X, i'), - X(I1- X),(64
[In (1 -X) -X] var(&^) -ua.
These formulae would certainly be correct asymptoticallyif the right-handsides were
divided by v and we were consideringpooled estimates from v completelyindependent
samples, with v tendingto infinity.But we are actually concernedwith one sample. Let us
see in what sense, if any, (6.4) can still be regarded as correct. If a and p are both large,
differentiationof any derivative of L by p changes its orderof magnitudeby a factor l/p,
and differentiation by a changes its order by a factor I/a. The second derivatives are
effectivelyconstantifin probability
8a+ a = ol),
a p
where&a = a-a, &p= p-p. Assuming(6.4) to be true,

da = o(Jja) p ( )

in probability. Hence formulae(6 4) are correctasymptoticallyif -a* c whilep is constant

or increases (or, moregenerally,has a positive lower bound).
But this result is not entirelysatisfactory,since a is a constant of the association being
observed,and althoughin many of the examples cited by Fisher etal. (1943) and Williams

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 377
(1944) a is fairlylarge, there is no logical necessityforit to be so.* The only adjustable
featureof the sample is the time of exposure (or the attractive power) of the trap, by
increasingwhich p may be increased. It is thereforeof interest to develop asymptotic
formulaevalid as p - so with ac constant.To the firstorder when p is large, the variance
of & given at (6.4) is (6
var(a) lnp (6.5)
while to the next orderof approximation

var(c)1
( I. (6.6)

Let us see whetherin fact one or both of these is correct,as p -* 00.

From (6 3), & is determinedby the equation

S = {ln(I+&)-ln}. (6 7)
The distributionof S is Poisson, with mean c ln(1 + p); the distributionof I is negative
binomial with mean ocpand exponent ac. Thus while the formerapproaches normalityas
p -x-00, and has coefficientof variation tending to zero, the latter does not. The distribution
of in (I + ac) also does not approach normalityas p -+ oo, but the asymptoticdistributionis
known (Anscombe, 1948), and we have
E{ln(I+co)} = lnp+i(ac)+O(p-fl), var{In(I +c)} = *'(z)+O(p-fl), (6 8)
where, is any quantitysuch that O< fi< 1, / < ac.We writenow &S, &I, &ln (I + ac) forthe
differences
between S, I, In (I + ac),and theirrespectivemean values. Then (6 7) can easily
be shown to give, in probabilityas p -* oo,

lnp-8{ + {lna-/f(a) -&ln (I +x)}_ (6-9)

Consideringthefirsttermsofnumeratorand denominator,we see at once thatthedistribution
of & is asymptoticallynormal and that (6.5) is correct.To see whether(6.6) is also correct,
we need to evaluate E{4S 41n(I + c)} and E{(aS)2 a1n (I + c)}. The joint distributionof S
and I has probability-generatingfunction
E (tSuI) = exp {ac(t-1 ) In (1 + p)-cat In (1 + p-pu)}. (6.10)
Writing x = {S- aln(I+p)}{aln(I+p)}-i and y = Ip-1)
we find,on expanding the characteristicfunctionof x and y forp large and applying the
Fourierinversionformula,the followingasymptoticcontinuousdistributionof x and y:

,a LI i) {8x(x23)+ax(lny-?x(a))}+O( )]dxdy. (6.11)

* The fact that cais not adjustable does not in itselfbar the use of an asymptoticformulavalid as
o--e o, since in any case such asymptoticformulaeare used as approximations.Even if the parameter
concernedis adjustable, onlyone value is usuallyavailable forconsideration,and not an infinitesequence
of values. An example of confusionon this point is the criticismby Kendall (1948) of a limitsituation
consideredby Jones (1948) in the theoryof systematicsampling. Jones gives a formulaforthe error
variancewhichis asymptoticallycorrectas thepopulationextenttendsto infinity, withconstantspacing
betweensamplepoints. Kendall, remarkingthatthepopulationextentis not in generaladjustable, gives
a formulaasymptoticallycorrectas the spacingbetweensamplepointstendsto zero,withthepopulation
extentconstant.As an approximationto the actual situation,Jones'sformulais the better(see Kendall,
1948,equations (20) and (25)).

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
378 Negative binomial and logarithmicseries distributions
Hence easily E {1S 6 In (I + cx)}= cY'(c) + O{(ln p)-A}) (6 12)
E {(&S)2 In
l (I +a)} = 0(1).

Finally, from(6.9), E(^) = a I+lna- (a)+o(1)} (6.13)

a
var
( )a ~~~~~~~~~~~~~~(6
Inp + a7f '(a) -2(1 + Ina - ib (a)) + o( l) *(- 4
On insertingthe asymptoticexpansions of 32(a) and #f'(a)in powers of a-', we obtain the
right-handside of (6.6) with remainderterm in the denominatorwhich is o(l) for both p
and a large. Thus (6 6) is not correct,to the ordersuggested,ifp is large,unlessa is also large.
The formulaefor var(5) that we have just discussed differfromone another only in
accuracy; they are all approximationsto the same true value, which has not been found.
An essentiallydifferent formulahas been givenby Fisher (Fisheretal. 1943). When expressed
in a formsimilarto (6 6), it is aln2
var(^) (6.15)
(lnp-1)2'
This formulais appropriateto a special type of comparison,namely, between estimates of
a for the same biological association derived fromsimilar nearby traps, where it may be
assumed that the individual species have exactly the same abundances (or at least the same
relative abundances), and the differencebetween the catches at any two traps arises solely
fromPoisson variation in the numberscaught of each species. In such a samplingprocess,
the estimate of a given by (6.7) is substantiallybiased, since only one set of relativeabund-
ances is involved. If we considerthe variation in the bias for all possible sets of relative
abundances followingthe limitingType III distributionassumed, the overall variance of &
is increasedfromFisher's value to that already considered.The largervarianceis appropriate
to comparingthe values of a fromobservationson different sorts of biological association,
involvingperhaps entirelydifferentfamiliesof species, and also to comparingvalues of a
fromobservationson the same sort of biological association observed at different seasons of
the year or in different years,when,even if the familiesof species are the same, the relative
abundances of the species are differentand may be supposed in aggregate to constitute
independent samples fromthe limitingType III distribution.Fisher's formula,in fact,is
not likely to be often useful,since if we desire to test whetherthe abundances of the in-
dividual species are the same (or in proportion)at a numberof traps it will be correctin the
firstplace to make direct comparisonsof counts of individual species, by a x2contingency-
table test, or by analysis of variance aftermaking a square-roottransformation.If it has
been established that the relative abundances of the species are not the same at different
traps, it is likely that they will differsufficiently
to appear in aggregateto be independent
samples fromthe hypotheticalparentType III distribution.They may indeed be so different
as to suggestquite different parent distributions,and a test of this point would be based on
the total variance of a.
For example, Williams gives figuresfor captures of Noctuidae during a period of three
months in 1933 at two traps (Fisher et al. 1943). One trap, on a roof-top,gave S = 58,
I = 1856, a = 11b37, = 163; the other,in a fielda quarter of a mile away, gave S = 40,
I = 929, a = 8.51, 5 = 109. Fisher's formula(6.15) forthe standard errorof eitherestimate
of a gives approximately0-67,indicatinga significantdifference betweenthem. Whetheror
not the relativeabundances of the species differedat the two traps would be moreefficiently

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 379
tested by comparingcounts ofindividual species. Formula (6.6) gives forthe standard error
of either estimate of a approximately 1-59,against which the observed differenceis not
significant.A significantdifferencein the richnessof the associations observed at the two
traps mightbe demonstrated,perhaps, by showingthat the fortyspecies caught in the field
trap had the same relative abundances (withinthe limits of Poisson variation in numbers
caught) as in the rooftrap, while the remainingspecies caught in the rooftrap were signi-
ficantlymore abundant, relativelyto the others,than in the fieldtrap (where the catches
werezero). Apart fromsome such argumentbased on comparingcatches ofindividualspecies,
we cannotconcludefromthefiguresforS and I alone that the biologicalassociationsobserved
at the two traps differediil diversityindex x.

SERIES FORMOF DISTRIBUTION

7. TESTS FORDEPARTUREFROMTHE LOGARITHMIC
Numerousalternativesto the logarithmicseriesdistributionsuggestthemselves. Fisher has
consideredthe negative binomial form(6.1) with k > 0 and X, k and NVunknown.We can
obtain otherthree-parameterdistributionsby replacingthe Type III distributionof species
abundances by any other distributionof a non-negativerandom variable. The situation is
the same as forthe heterogeneousPoisson samplingconsideredin ? 2, except that nois not
observed and N is an unknownparameterrequiringestimation.
Fisher gives a test fordeparturefromthe logarithmicformof distributiontowards that
at (6. 1) withk > 0. The appropriatestatistic,in addition to S and I, is

r=2 r( 2 3 r-1) 71
ofwhichtheexpected value is I + may be takenas a testcriterion,
and its samplingvariance can be investigated
(1 p)]2. Thus the methodsalready indicated. If we are
byJ-S2/(2
asymptoticresult when a -b oo, we may consider the matrix of
content with a first-order
acln &)
expectations of second derivatives of the logarithm of the likelihood function(6.1), the
being with respectto X, a, and k. Settingk = 0, we obtain
differentiations

La JS2 ) X-7{In (I1-X)}3 {In (I1-X) -4X} + 0(X) { -In (I1- X) - XI] (7
2"a" I~-n(1- X) - X,()

where (X r I2( +I+ 4 +(/-1)2)+ .r (7.3)

When p is large Os(X) A lnp -B, (7.4)

where,in termsof the Riemann c-function,
A = (2) = -6449, B = 2C(3) = 2-4041; (7.5)
a factor1+ O(p-1),as in (6 6),
and hence,ignoring

var _I\2 ca[-2(1np)3

- (lnp -4) + (A lnp -B) (lnp - 1)] (76)
I~~~np-1I
Fisher denotes the right-handside of (7.2) by i and has tabulated it.
Applicationsof the distributionmade by VWilliams suggestthat anothersort of departure
fromthe logarithmicseries formmay be worthinvestigating. In a complex association it
mighthappen that while certain componentsof the association exhibit logarithmicseries

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
380 Negativebinomialand logarithmic
seriesdistributions
distributionsthe whole does not, since the component distributionshave differentpara-
meters. Let cj, Xi, pi relate to the ith componentassociation (i = 1,2, ..., 4), and consider
a logarithmicseriesdistributionwithparametersao, X0,po, chosento give the same expected
numbersof species and individuals as in the whole association. Then
V

ao In (1 + po) =E: In(1 +pi))

*= 1 (7 7)
IJ

Logo EXii

Hence n(+p0) = Wi (1+Pi) (7.8)

Po i=1 Pi

where = cXiPi/> jpj~

so that In (1 +po)/po is a weighted mean of In (1 +pi)/pi. If E and E* denote respectively
expectations for the actual distributionand for the fittedlogarithmicseries distribution,
we have E )r-1
}i (n,.)= ltPi r p) 7

E*(nfrV){= cxiPi} r( +p )r |

We thereforeconsider Y =pr-l(1 +p)-r as a functionof Z =pal In (1 +p). When r = 1, we

findd2YldZ2 > 0 forall p, so that Y is a convex functionof Z, and thereforeE*(nj) < E (nj),
withequality onlyifthe pi are all equal. For r > 2, the signof d2YldZ2 depends on p and can
be studied in detail foreach r. It is not difficultto show that E* (ar) < E (nr)fora range of
small values of r and also forlarge values of r, while the inequality is reversedin a middle
range (assumingthe pi not all equal). If all the pi are large and not veryunequal, the values
of r wherethe inequality changes are roughly

2 lPo and 2po.

The appropriatestatisticfordetectinga small degreeof inequality in the pi is easily seen

fromthe likelihoodfunctionto be Z nrr(r- 1), ofwhichthe expectedvalue is acp2 whenthere
r=2

is no heterogeneity,i.e. whenpi = p forall i and a = a. Hence we may take

i= 1
00
x nrr(r-1)
w r=2
(7*10)
12

as test criterion.It is easy to show that, asymptoticallyforlarge a,

E(W)1, var(W)A[, - -ln(l2X)-X] (7.11)

or, whenp is large also, var (IW) a[2- l i -]. (7.12)

The test is a limitingformof Test 2 of ? 4.

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
 F. J. ANSCOMBE 381
Several distributionsof classificationof species into genera given by Williams (1944) seem
to showthis kind of heterogeneity.Thus forCoccidae of the worldclassifiedby MacGillivray
(Williams's Table 11), S = 352, I = 1763, a' = 132*2,p = 13 34, W = 2 76, s.e. (W) = 0.11.
The observednumberofmonotypicgeneran1is 181,whichis much above its expected value,
123*0,on the basis of the logarithmicseries distribution. From r = 2 to r = 25, roughly,
the nrare on the whole less than expected, whileforhigherr theyappear again to be greater
than expected.
W is closelyrelated to the characteristicK of Yule (1944). Thus

K = 1O,000W = 10,000 .n(2- (7.13)

It is easy to show that, asymptoticallyforlarge a,

E (K) = 10,00O0'-, var (K) = 2 x 108[a3X2]-1, (7.14)
if the observations are drawn froma logarithmicseries distribution. Simpson (1949) has
shown that, asymptoticallyforlargep, but with axnot assumed large,
E (K) = 10,000(x? 1)-'. (7.15)
As a statisticforestimatinga, K is oflow efficiency.Its properfunctionis to testdlistribution
shape.

I owe my interestin this subject to stimulatingconversationswith Dr C1.B. Williams.

Severalpersonshave offered in particularMr1). G. Kendall,
helpfulsuggestionsand criticisms,
Mr J. G. Skellam and my colleagues at Cambridge,Dr J. Wishart,Dr H. E. Daniels and
Mr D. V. Lindley. The work was begun at Rothamsted Experimental Station, and some
help with the computationson whichthe figureis based was given by Mr B. M. Church.

REFERENCES

AITKEN, A. C. (1939). StatisticalMathematics.Edinburgh: Oliver and Boyd.

ANSCOMBE, F. J. (1948). Biometrika,35, 246.
ANSCOMBE, F. J. (1949). Biometrics, 5, 165.
ANSCOMBE, F. J. (1950). Ann. Apple. Biol. 37, 286.
BEALL, G. (1940). Ecology,21, 460.
BUXTON, P. A. (1940). Parasitology,32, 296.
CERNUSOHI, F. & CASTAGNETTO, L. (1946). Ann. Math. Statist.17, 53.
FELLER, W. (1943). Ann. Math. Statist.14, 389.
FINNEY, D. J. (1941). J. R. Statist.Soc. Suppl. 7, 155.
FISHER, R. A. (1931). BritishAssociationMathematicalTables, 1, xxvi.
FISHER, R. A. (1941). Ann. Eugen., Lond., 11, 182.
FISHER, R. A., CORBET, A. S. & WILLIAMS, C. B. (1943). J. Anim. Ecol. 12, 42.
GADDUM, J. H. (1945). Nature,Lond., 156, 463.
GREENWOOD, M. & YULE, G. U. (1920). J. R. Statist.Soc. 83, 255.
HALDANE, J. B. S. (1941). Ann. Eugen., Lond., 11, 179.
HALDANE, J. B. S. (1945). Biometrika, 33, 222.
HALDANE, J. B. S. (1949). J. R. Statist.Soc. B, 11, 1.
IRWIN, J. 0. (1941). J. R. Statist.Soc. Suppl. 7, 101.
JoNEs,A. E. (1948). Biometrika,35, 283.
KENDALL, D. G. (1949). J. B. Statist.Soc. B, 11, 230.
KENDALL, M. G. (1943). The AdvancedTheoryofStatistics,1. London: Griffin.
KENDALTLt, M. G. (1948). Biometrika, 35, 291.

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions
382 Negativebinomialand logarithmic
seriesdistributions
LUDERS,R. (1934). Biometrilca, 26, 108.
McKENDRICK, A. G. (1914). Proc. Lond. Math. Soc. 13, 401.
NEWBOLD,E. M. (1927). J. R. Statist.Soc. 90, 487.
NEYMAN, J. (1939). Ann. Math. Statist.10, 35.
P6LYA, G. (1930). Ann. Inst. Poincar4,1, 117.
PRESTON, F. W. (1948). Ecology, 29, 254.
QUENOUILLE, M. H. (1949). Biometrics, 5, 162.
SHENTON, L. R. (1949). Biometrika,
36, 450.
SIMPSON,E. H. (1949). Nature,Lond., 163, 688.
THOMAS, M. (1949). Biometrika,36, 18.
WILLIAMS, C. B. (1937). Ann. Apple. Biol. 24, 404.
WiLuLms, C. B. (1944). J. Ecol. 32, 1.
WISHART, J. (1947). J. Inst. Actu. Stud. Soc. 6, 140.
YuIE, G. U. (1910). J. R. Statist.Soc. 73, 26.
YuLE, G. U. (1944). The StatisticalStudyofLiteraryVocabulary.CambridgeUniversityPress.

This content downloaded from 193.104.110.40 on Thu, 19 Jun 2014 19:56:38 PM

All use subject to JSTOR Terms and Conditions