Vidakovic B. Statistics for Bioengineering Sciences: With Matlab and WinBugs Support

Подождите немного. Документ загружается.

146 5 Random Variables

Example 5.8. Surviving Pairs.



Daniel Bernoulli (1700–1782), a nephew

of Jacob Bernoulli, posed and solved the following problem. If among N mar-

ried pairs there are m random deaths, what is the expected number of intact

marriages?

Suppose that there are N pairs of balls denoted by 1,1, 2,2, . . . , N,N. If

m balls are selected at random and removed, what is the expected number of

intact pairs? Consider the pair i. Deﬁne a Bernoulli random variable Y

equal

to 1 if pair i remains intact after the removal of m balls, and 0 otherwise. Then

the number of unaffected pairs N

would be the sum of all Y

, for i =1,..., N.

The probability that pair i is not affected by the removal of m balls is

2N−2

(2N−2)(2N−3)...(2N−2−m+2)(2N−2−m+1)

2N(2N−1)...(2N−m+2)(2N−m+1)

(2N −m)(2N −m −1)

2N(2N −1)

and it is equal to

. If N

is the number of unaffected pairs, then

= Y

+···+Y

= EY

+EY

+···+EY

= NEY

(2N −m)(2N −m −1)

2(2N −1)

For example, if among N

= 1000 couples there are 100 random deaths,

then the expected number of unaffected couples is 902.4762. If among N

1000 couples there are 1936 deaths, then a single couple is expected to remain

intact.

Even though N

is the sum of N Bernoulli random variables Y

, each with

the same probability p

(2N− m)(2N−m−1)

2N(2N−1)

, it does not have a binomial distribu-

tion due to the dependence among Y

5.3.3 Hypergeometric Distribution

Suppose a box contains m balls, k of which are white and m −k of which are

black. Suppose we randomly select and remove n balls from the box without

replacement, so that when sampling is ﬁnished, there are only m

−n balls left

in the box. If X is the number of white balls in n selected, then the probability

that X

= x is

(x) =

¢¡

m−k

−x

, x ∈{0,1,...,min{n, k}}.

Random variable X is called hypergeometric and denoted by X

∼H G (m, k, n),

where m, k, and n are integer parameters.

5.3 Some Standard Discrete Distributions 147

This PMF can be deduced by counting rules. There are

different ways

of selecting the n balls from a box with a total of m balls. From these (each

equally likely), there are

ways of selecting x white balls from the k white

balls in the box and, similarly,

m−k

−x

ways of choosing the black balls. The prob-

ability is the ratio of these two numbers. The PDF and CDF of

H G (40,15,10)

are shown in Fig. 5.4.

It can be shown that the mean and variance for the hypergeometric distri-

bution are, respectively,

µ =

and

¶µ

m −k

m −n

m −1

The MATLAB commands for hypergeometric CDF, PDF, quantile, and a ran-

dom number are

hygecdf, hygepdf, hygeinv, and hygernd. WinBUGS does not

have a built-in command for a hypergeometric distribution.

Example 5.9. CASES. In a group of 40 people, 15 are “CASES” and 25 are

“CONTROLS.” A sample of 10 subjects is selected [(A) with replacement and

%Solution

%(A) - with replacement (binomial case);

%Let X be the number of CASES. The event

%X is at least 2 is the complement of X <= 1.

disp(’(A) Bin(10, 15/40): P(X >= 2)’); 1 - binocdf(1, 10, 15/40)

% ans = 0.9363

% or

1 - binopdf(0, 10, 15/40) - binopdf(1, 10, 15/40)

% ans = 0.9363

%B - without replacement (hypergeometric case) hygecdf(x, m, k, n)

% where m size of population,

% k - number of cases among m, and n sample size.

disp(’(B) HyGe(40,15,10): P(X >=2)’); 1 - hygecdf(1, 40, 15, 10)

% ans = 0.9600, or

1 - hygepdf(0, 40, 15, 10)- hygepdf(1, 40, 15, 10)

% ans = 0.9600

Example 5.10. Capture-Recapture Models. Suppose that an unknown num-

ber m of animals inhabit a particular region. To assess the population size,

ecologists often apply the following capture-recapture scheme. They catch k

animals, tag them and release them back into the region. After some time,

when the tagged animals are expected to be mixed well with the untagged, a

second catch of size n is made. Suppose that x animals in the second sample

are found to be tagged.

If catching any animal is assumed equally likely, the number x of tagged

animals in the second sample is hypergeometric

H G (m, k, n). Ecologists use

(B) without replacement]. Find the probability

P(at least 2 subjects are CASES).

148 5 Random Variables

0 1 2 3 4 5 6 7 8 9 10

0.1

0.2

0 2 4 6 8 10

0.2

0.4

0.6

0.8

Fig. 5.4 The PDF and CDF for hypergeometric distribution with m =40, k =15, and n =10.

the observed ratio x/n as an approximation to k/m, from which m is estimated

k ×n

A statistically better estimator of m (known as the Schnabel formula) is given

(k +1) ×(n +1)

(x +1)

−1.

In epidemiology and public health capture-recapture methods use multi-

ple, routinely collected, computerized data sources to estimate various popu-

lation indices.

For example, Gjini et al (2004) investigated the number of matching

records of pneumococcal meningitis among adults in England by comparing

data from Hospital Episode Statistics (HES) and the Public Health Labora-

tory Services reconciled laboratory records (RLR). The time period covered

was April 1996 to December 1999. The authors found 646 records in RLR and

737 in HES, and matching based on demographic information was possible in

296 cases.

By the capture-recapture method the estimated incidence is

= 646 ·

737/296 = 1608.5 ≈ 1609. If Schnabel’s formula is used, then

m ≈1607. Thus,

the total incidence of of pneumococcal meningitis in England between April

1996 to December 1999 is estimated to be 1607.

5.3 Some Standard Discrete Distributions 149

For large m, the hypergeometric distribution is close to binomial. More

precisely, when m

→∞ and k/m → p, the hypergeometric distribution with pa-

rameters (m, k,n) approaches a binomial with parameters (n, p) for any value

of x between 0 and n.

format long

disp(’(A)=(B) for large population’);

1 - binocdf(1, 10, 150000/400000) %ans = 0.936335370875895

1 - hygecdf(1, 400000, 150000, 10) %ans = 0.936337703719839

We will use the hypergeometric distribution later in the book, (i) in the

Fisher exact test (p. 546), and in Logrank test (p. 713).

5.3.4 Poisson Distribution

The PMF for the Poisson distribution is

(x) =

−λ

, x =0,1,2,... ,

which is denoted by X

∼ P oi(λ). This distribution is named after Simeon

Denis Poisson (1781–1840), French mathematician, geometer and physicist

(Fig. 5.3b).

From the moment-generating function m

(t) =exp{λ(e

−1)} we have EX =

λ and Var X =λ; the mean and the variance coincide.

The sum of a ﬁnite independent set of Poisson variables is also Pois-

son. Speciﬁcally, if X

∼ P oi(λ

), then Y = X

+···+ X

is distributed as

P oi(λ

+···+λ

). If X

∼ P oi(λ

) and X

∼ P oi(λ

) are independent, then

the distribution of X

given that X

= n is binomial B in

+λ

(Exer-

cise 5.5).

Furthermore, the Poisson distribution is a limiting form for a binomial

model, i.e.,

lim

n,np→∞,λ

(1 − p)

n−x

−λ

. (5.6)

The MATLAB commands for Poisson CDF, PDF, quantile, and random num-

ber are

poisscdf, poisspdf, poissinv, and poissrnd. In WinBUGS the Poisson

distribution is denoted as

dpois(lambda).

Example 5.11. Poisson Model for EBs. After 7 days of aggregation, the mi-

croscopy images of 2000 embryonic bodies (EBs) are used to assess their sur-

face area size. The probability that the area of a randomly selected EB exceeds

the critical size S

is 0.001.

(a) Find the probability that the areas of exactly three EBs exceed the

critical size.

150 5 Random Variables

(b) Find the probability that the number of EBs exceeding the critical size

is between three and eight, inclusively.

We use a Poisson approximation to binomial probabilities since n is large,

p is small, and product n p is moderate.

%Solution

disp(’Poisson(2): P(X=3)’); poisspdf(3, 2)

%ans= 0.1804

disp(’Poisson(2): P(3 <= X <= 8)’); poisscdf(8, 2)-poisscdf(2, 2)

%ans= 0.3231

Figure 5.5 shows the PMF and CDF of the P oi(2) distribution.

0 1 2 3 4 5 6 7 8 9 10

0.1

0.2

0.3

0 2 4 6 8 10

0.2

0.4

0.6

0.8

Fig. 5.5 (Top) Poisson probability mass function. (Bottom) Cumulative distribution function

for

λ =2.

In the binomial sampling scheme, when n →∞ and p →0, so that np →λ,

binomial probabilities converge to Poisson probabilities.

The following MATLAB simulation demonstrates the convergence. In the

binomial distribution n is increasing as 2,000, 200,000, 20,000,000 and p is

decreasing as 0.001, 0.00001, 0.0000001, so that np remains constant and

equal to 2. These binomial probabilities for X

=3 are compared to the Poisson

distribution with the parameter

λ =2.

disp(’P(X=3) for Bin(2000, 0.001), Bin(200000, 0.00001), ’);

disp(’ Bin(20000000, 0.0000001), and Poi(2) ’);

format long

5.3 Some Standard Discrete Distributions 151

binopdf(3, 2000, 0.001) % 0.180537328031786

binopdf(3, 200000, 0.00001) % 0.180447946554779

binopdf(3, 20000000, 0.0000001) % 0.180447058859339

poisspdf(3, 2) % 0.180447044315484

format short

Example 5.12. Cold. Suppose that the number of times during a year that an

individual catches a cold can be modeled by a Poisson random variable with an

expectation of 4. Further, suppose that a new drug based on vitamin C reduces

this expectation to 3 (but the distribution still remains Poisson) for 90% of the

population but has no effect on the remaining 10% of the population. We will

calculate

(a) the probability that an individual taking the drug has two colds in a

year if they are part of the population that beneﬁts from the drug;

(b) the probability that a randomly chosen individual has two colds in a

year if they take the drug; and

part of the population that beneﬁts from the drug, given that they had two

colds in the year during which they took the drug.

poisspdf(2,3) %(Cold (a))

%ans = 0.2240

poisspdf(2,3)

0.90 + poisspdf(2,4)

0.10 %(Cold (b))

%ans = 0.2163

poisspdf(2,3)

0.90/(poisspdf(2,3)

0.90 + ...

poisspdf(2,4)

0.10) %(Cold (c))

%ans = 0.9323

5.3.5 Geometric Distribution

Suppose that independent trials are repeated and that in each trial the prob-

ability of a success is equal to 0

< p < 1. We are interested in the number of

failures X before the ﬁrst success. The number of failures is a random variable

with a geometric

G e(p) distribution. Its PMF is given by

(x) = p(1− p)

, x =0,1,2,... .

The expected value is

EX =(1−p)/p and the variance is Var X =(1−p)/p

. The

moments can be found either directly or by the moment-generating function,

which is

(t) =

1 −(1 − p)e

152 5 Random Variables

The geometric random variable possesses a “memoryless” property. That is,

if we condition on the event X

≥ m, for some nonnegative integer m, then for

≥ m, P(X ≥ n|X ≥ m) = P(X ≥ n − m) (Exercise 5.17). The MATLAB com-

mands for geometric CDF, PDF, quantile, and random number are

geocdf,

geopdf, geoinv, and geornd. There are no special names for the geometric dis-

tribution in WinBUGS; the negative binomial can be used as

dnegbin(p,1).

If instead of the number of failures before the ﬁrst success (X ) one is inter-

ested in the total number of experiments until the ﬁrst success (Y ), then the

relationship is simple: Y

= X +1. In this formulation of the geometric distri-

bution, Y

∼G eom(p), EY =EX +1 = 1/p, and Var Y =Var X =(1 − p)/p

≤160 mg/dL, respectively. According to NHANES III, the prevalence of

CASES among white male Americans aged 20 or older (target population) is

=20%. Subjects are sampled (when the population is large, it is unimportant

if the sampling is done with or without replacement) until the ﬁrst CASE is

found. The number of CONTROLS sampled before ﬁnding the ﬁrst CASE is a

geometric random variable with parameter p

=0.2 (Fig. 5.6).

(a) Find the probability that seven CONTROLS will be sampled before we

come across the ﬁrst CASE.

(b) Find the probability that the number of CONTROLS before the ﬁrst

CASE will fall between four and eight, inclusively.

disp(’X ~ Geometric(0.2):P(X=7)’);

geopdf(7, 0.2)

%ans=0.0419

disp(’X ~ Geometric(0.2):P(4 <= X <= 8)’);

geocdf(8, 0.2) - geocdf(3,0.2)

%ans=0.2754

5.3.6 Negative Binomial Distribution

The negative binomial distribution was formulated by Montmort (1714). Sup-

pose we are dealing with independent trials again. This time we count the

number of failures observed until a ﬁxed number of successes (r

≥1) occur. Let

p be the probability of success in a single trial.

If we observe r consecutive successes at the start of the experiment, then

the count of failures is X

=0 and P(X =0) = p

If X

= x, then we have observed x failures and r successes in x + r tri-

als. There are

x+r

different ways of arranging x failures in those x + r tri-

als, but we can only be concerned with those arrangements in which the last

Example 5.13. CASES I. Let a subject constitute either a CASE or a CON-

TROL depending on whether the subject’s level of LDL cholesterol is

160 mg/dL

5.3 Some Standard Discrete Distributions 153

0 5 10 15 20

0.05

0.1

0.15

0.2

0 5 10 15 20

0.2

0.4

0.6

0.8

Fig. 5.6 Geometric (Top) PMF and (Bottom) CDF for p =0.2.

trial ended in a success. So there are really only

x+r−1

equally likely arrange-

ments. For any particular arrangement the probability is p

(1−p)

. Therefore,

the PMF is

(x) =

r +x −1

(1 − p)

, x =0,1,2,... .

Sometimes this PMF is stated with

r+x−1

−1

in place of the equivalent

r+x−1

This distribution is denoted as X

∼N B(r, p).

From its moment-generating function

(t) =

1 −(1 − p)e

the expectation of a negative binomial random variable is

EX = r(1− p)/p and

its variance is

Var X = r(1 − p)/p

Since the negative binomial X

∼ N B(r, p) is a convolution (a sum) of r

independent geometric random variables, X

+···+Y

, Y

∼G (p), the

mean and variance of X can be easily derived from the mean and variance

of its geometric components Y

, as in (5.1) and (5.2). Note also that m

(t) =

(

(t)

)

, where m

(t) =

1−(1−p)e

is the moment-generating function of the

component Y

in the sum. This is a consequence of the fact that a moment-

generating function for a sum of independent random variables is the product

of the moment-generating functions of the components, see (5.4).

The distribution remains valid if r is not an integer, although an inter-

pretation involving r successes is lost. For an arbitrary nonnegative r, the

154 5 Random Variables

distribution is called a generalized negative binomial or a Polya distribution.

The constant

r+x−1

(r+x−1)!

x!(r−1)!

is replaced by

Γ(r+x)

x!Γ(r)

, keeping in mind that

Γ(n) = (n −1)! when n is an integer. The generalized negative binomial dis-

tribution is used in ecology for inference about the abundance of species in

nature.

The MATLAB commands for negative binomial CDF, PDF, quantile, and

random number are

nbincdf, nbinpdf, nbininv, and nbinrnd. In WinBUGS the

negative binomial distribution is denoted as

dnegbin(p,r). Note the opposite

order of parameters r and p.

Example 5.14. CASES II. Assume as in Example 5.13 that the prevalence of

“CASES” in a large population is p

= 20%. Subjects are sampled, one by one,

until seven CASES are found and then the sampling is stopped.

(a) What is the probability that the number of CONTROLS among all sam-

pled subjects will be 18?

(b) What is the probability of observing more than the “expected number”

of CONTROLS?

The number of CONTROLS X among all sampled subjects is a negative

binomial, X

∼N B(7, 0.2).

P(X =18) =

25 +7 −1

0.2

(1 −0.2)

8 =0.0310.

Also,

nbinpdf(18,7,0.2)=0.0310.

Thus, with a probability of 0.031 the number of CONTROLS sampled, be-

fore seven CASES are observed, is equal to 18.

(b) The expected number of CONTROLS is

EX =7

0.8

0.2

=28. The probability

of X

> EX is P(X > 28) = 1 −P(X ≤ 28) = 1 −

7+x−1

0.8

0.2

= 0.4328. In

MATLAB

P(X >28) is calculated as 1-nbincdf(28,7,0.20)=0.4328.

The tail probabilities of a negative binomial distribution can be expressed

by binomial probabilities. If X

∼N B(r, p), then

P(X > x) =P(Y < r),

where Y

∼B in(x +r, p). In words, if we have not seen r successes after seing

x failures, then in x

+r experiments the number of successes will be less than

r. In part (b) of the previous example, r

=7, x =28, and p =0.20, so

1 - nbincdf(28, 7, 0.20) % 0.4328

binocdf(7-1, 28+7, 0.20) % 0.4328

5.3 Some Standard Discrete Distributions 155

5.3.7 Multinomial Distribution

The binomial distribution was developed by counting the occurrences two com-

plementary events, A and A

, in n independent trials. Suppose, instead, that

each trial results in one of k

> 2 mutually exclusive events, A

,... A

, so that

S = A

∪···∪ A

. One can deﬁne the vector of random variables (X

,... , X

)

where a component X

counts how many times A

appeared in n trials. The

deﬁned random vector is called multinomial.

The probability mass function for (X

,... , X

) is

,...,X

,..., x

) =

!···x

···p

where p

+···+ p

=1 and x

+···+x

= n. Since p

=1 − p

−···− p

k−1

, there

are k

−1 free parameters to characterize the multinomial distribution, which

is denoted by X

=(X

,... , X

) ∼M n(n, p

,... , p

The mean and variance of the component X

are the same as in the bino-

mial case. It is easy to see that the marginal distribution for a component X

is binomial since the events A

,... , A

can be grouped as A

, A

. Therefore,

E(X

) = np

, Var (X

) = np

(1 − p

). The components X

are dependent since

they sum up to n. For i

6= j, the covariance between X

and X

Cov(X

, X

) =EX

−EX

=−np

. (5.7)

This is easy to verify if X

and X

are represented as the sums of Bernoullis

+···+Y

and Y

+···+Y

, respectively. Since Y

=0 (in a

single trial A

and A

cannot occur simultaneously), it follows that

=(n

−n)p

Since

= n

, the covariance in (5.7) follows.

If X

= (X

, X

,... , X

) ∼ M n(n, p

, p

,... , p

), then X

= (X

+ X

,... , X

)

∼M n(n, p

+ p

,... , p

). This is called the fusing property of the multinomial

distribution.

If X

∼ P oi(λ

), X

∼ P oi(λ

), ..., X

∼ P oi(λ

) are n independent Pois-

son random variables with parameters

,... , λ

, then the conditional distri-

bution of X

, X

,... , X

, given that X

+···+X

= n, is M n(n, p

,... , p

where p

=λ

/(λ

+λ

+···+λ

). This fact is used in modeling contingency ta-

bles with a ﬁxed total and will be discussed in Chap. 14.

In MATLAB, the multinomial PMF is calculated by

mnpdf(x,p), where x

is a 1

×k vector of values, such that

= n, and p is a 1 × k vector of

probabilities, such that

=1.

For example,

%If n=2, Multinomial is Binomial

mnpdf([5 15],[0.6 0.4])

%ans = 0.0013