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

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

166 5 Random Variables

(x) =

Γ(r)

r−1

−λx

, x ≥0,

0, else .

The parameter r

> 0 is called the shape parameter, while λ > 0 is the rate

parameter. Figure 5.11b shows gamma densities for (r,

λ) =(1,1/3),(2,2/3), and

(20,2).

0 1 2 3 4 5 6 7 8

Γ(x)

0 5 10 15

0.05

0.1

0.15

0.2

0.25

0.3

PDF

Ga(1, 1/3)

Ga(2, 2/3)

Ga(20,2)

(a) (b)

Fig. 5.11 (a) Gamma function,

Γ(x). The red dots are values of the gamma function at inte-

gers,

Γ(n) =(n −1)!; (b) Gamma densities: G a(1, 1/3), G a(2,2/3), and G a(20, 2).

The moment-generating function is m(t) =

(

λ/(λ −t)

)

, so in the case r =1,

the gamma distribution becomes the exponential distribution. From m(t) we

have

EX = r/λ and Var X = r/λ

If X

,... , X

are generated from an exponential distribution with (rate) pa-

rameter

λ, it follows from m(t) that Y = X

+···+ X

is distributed as gamma

with parameters

λ and n; that is, Y ∼G a(n,λ). More generally, if X

∼G a(r

,λ)

are independent, then Y

= X

+···+ X

is distributed as gamma with param-

eters

λ and r = r

+···+r

,; that is, Y ∼G a(r,λ) (Exercise 5.16).

Often, the gamma distribution is parameterized with 1/

λ in place of λ,

and this alternative parametrization is used in MATLAB deﬁnitions. The CDF

in MATLAB is

gamcdf(x,r,1/lambda), and the PDF is gampdf(x,r,1/lambda).

The function

gaminv(p,r,1/lambda) computes the pth quantile of the G a(r,λ)

random variable. In WinBUGS,

G a(n, λ) is coded as dgamma(n,lambda).

5.5.5 Inverse Gamma Distribution

Random variable X is said to have an inverse gamma I G (r,λ) distribution

with parameters r

>0 and λ >0 if its density is given by

5.5 Some Standard Continuous Distributions 167

(x) =

(

Γ(r)x

r+1

−λ/x

, x ≥0,

0, else .

The mean and variance of X are

EX = λ/(r −1), r > 1, and Var X = λ

/[(r −

(r − 2)], r > 2, respectively. If X ∼ G a(r,λ), then its reciprocal X

−1

I G (r, λ) distributed. We will see that in the Bayesian context, the inverse

gamma is a natural prior distribution for a scale parameter.

5.5.6 Beta Distribution

We ﬁrst deﬁne two special functions: beta and incomplete beta. The beta func-

tion is deﬁned as B(a, b)

a−1

(1 −t)

b−1

dt =Γ(a)Γ(b)/Γ(a +b). In MATLAB:

beta(a,b). An incomplete beta is B(x,a, b) =

a−1

(1 − t)

b−1

dt, 0 ≤ x ≤ 1. In

MATLAB,

betainc(x,a,b) represents the normalized incomplete beta deﬁned

as I

(a, b) =B(x,a, b)/B(a, b).

The density function for a beta random variable is

(x) =

B(a,b)

a−1

(1 −x)

b−1

, 0 ≤ x ≤1,

0, else ,

and B is the beta function. Because X is deﬁned only in the interval [0,1],

the beta distribution is useful in describing uncertainty or randomness in pro-

portions or probabilities. A beta-distributed random variable is denoted by

∼ B e(a, b). The standard uniform distribution U (0,1) serves as a special

case with (a, b)

=(1,1). The moments of beta distribution are

Γ(a +k)Γ(a +b)

Γ(a)Γ(a +b +k)

a(a +1)...(a +k −1)

(a +b)(a +b +1)...(a +b +k −1)

so that

E(X ) = a/(a +b) and Var X = ab/[(a +b)

(a +b +1)].

In MATLAB, the CDF for a beta random variable (at x

∈ (0,1)) is com-

puted as

betacdf(x,a,b), and the PDF is computed as betapdf(x,a,b). The

pth percentile is

betainv(p,a,b). In WinBUGS, the beta distribution is coded

dbeta(a,b).

To emphasize the modeling diversity of beta distributions, we depict den-

sities for several choices of (a,b), as in Fig. 5.12.

If U

,... ,U

is a sample from a uniform U (0,1) distribution, then

the distribution of the kth component in the ordered sample is beta, U

(k)

∼

e(k, n − k +1), for 1 ≤ k ≤ n. Also, if X ∼ G (m,λ) and Y ∼ G (n,λ), then

X /(X

+Y ) ∼B e(m, n).

168 5 Random Variables

0 0.5 1

0.5

0 0.5 1

Fig. 5.12 Beta densities for (a, b) as (1/2, 1,2), (1,1), (2,2), (10,10), (1,5), (1, 0.4), (3,5), (50,

30), and (5000, 3000).

5.5.7 Double Exponential Distribution

A random variable X has double exponential DE (µ,λ) distribution if its den-

sity is given by

(x) =

−λ|x−µ|

, −∞< x <∞,λ >0.

The expectation of X is

EX =µ, and the variance is Var X =2/λ

. The moment-

generating function for the double exponential distribution is

m(t)

µt

−t

, |t|<λ.

The double exponential distribution is also known as the Laplace distribution.

If X

and X

are independent exponential E (λ), then X

− X

is distributed

DE (0,λ). Also, if X ∼ DE (0,λ), then |X | ∼ E (λ). In MATLAB the double

exponential distribution is not implemented since it can be readily obtained

by folding the exponential distribution. In WinBUGS,

DE (0,λ) is coded as

ddexp(lambda).

5.5 Some Standard Continuous Distributions 169

5.5.8 Logistic Distribution

The logistic distribution is used for growth models in pharmacokinetics, re-

gression with binary responses, and neural networks, to list just a few model-

ing applications.

The logistic random variable can be introduced by a property of its CDF

expressed by a differential equation. Let F(x)

=P(X ≤ x) be the CDF for which

(x) = F(x) ×(1 −F(x)). One interpretation of this differential equation is as

follows. For a Bernoulli random variable

1(X ≤ x) =

1, X ≤ x

0, X

> x

, the change in

probability of 1 as a function of x is equal to its variance. The solution in the

class of CDFs is

F(x)

1 +e

−x

1 +e

which is called the logistic distribution. Its density is

f (x)

(1 +e

)

Graphs of f (x) and F(x) are shown in Fig. 5.13. The mean of the distribution

−6 −4 −2 0 2 4 6

0.05

0.1

0.15

0.2

0.25

PDF

logistic

normal

−6 −4 −2 0 2 4 6

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

CDF

logistic

normal

(a) (b)

Fig. 5.13 (a) Density and (b) CDF of logistic distribution. Superimposed (dotted red) is the

normal distribution with matching mean and variance, 0 and

/3, respectively.

is 0 and the variance is π

/3. For a more general logistic distribution given by

the CDF

F(x)

1 +e

−(x−µ)/σ

the mean is

µ, variance π

/3, skewness 0, and kurtosis 21/5. For the higher

moments one can use the moment generating function

170 5 Random Variables

m(t) =exp{µt}B

(

1 −σt, 1 +σt

)

where B is the beta function. In WinBUGS the logistic distribution is coded as

dlogis(mu,tau), where tau is the reciprocal of σ.

If X has a logistic distribution, then log(X ) has a log-logistic distribution

(also known as the Fisk distribution). The log-logistic distribution is used in

economics (population wealth distribution) and reliability.

The logistic distribution will be revisited in Chap. 17 where we deal with

logistic regression.

5.5.9 Weibull Distribution

The Weibull distribution is one of the most important distributions in survival

theory and engineering reliability. It is named after Swedish engineer and

scientist Waloddi Weibull (Fig. 5.14) after his publication in the early 1950s

(Weibull, 1951).

Fig. 5.14 Ernst Hjalmar Waloddi Weibull (1887–1979).

The density of the two-parameter Weibull random variable X ∼ W ei(r,λ)

is given as

f (x)

=λrx

r−1

−λx

, x >0.

The CDF is given as F(x)

=1 −e

−λx

. Parameter r is the shape parameter and

λ is the rate parameter. Both parameters are strictly positive.

In this form, Weibull X

∼W ei(r, λ) is a distribution of Y

for Y exponential

E (λ).

In MATLAB the Weibull distribution is parameterized by a and r, as in

f (x)

−r

r−1

−(x/a)

, x >0. Note that in this parameterization, a is the scale

parameter and relates to

λ as λ = a

−r

. So when a =λ

−1/r

, the CDF in MATLAB

wblcdf(x,a,r), and the PDF is wblpdf(x,a,r). The function wblinv(p,a,r)

computes the pth quantile of the W ei(r,λ) random variable. In WinBUGS,

W ei(r,λ) is coded as dweib(r,lambda).

5.5 Some Standard Continuous Distributions 171

The Weibull distribution generalizes the exponential distribution (r = 1)

and Rayleigh distribution (r

=2).

The mean of a Weibull random variable X is

EX =

Γ(1+1/r)

1/r

=aΓ

1 +

, and

the variance is

Var X =

Γ(1+2/r)−Γ

(1+1/r)

2/r

(Γ

1 +

−

1 +

). The kth mo-

ment is

Γ(1+k/r)

k/r

1 +

Figure 5.15 shows the densities of the Weibull distribution for r

=2 (blue),

= 1 (red), and r = 1/2 (black). In all three cases, λ = 1/2. The values for the

scale parameter a

=λ

−1/r

are

2,2, and 4, respectively.

0 1 2 3 4 5 6

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

r=2

r=1

r=1/2

Fig. 5.15 Densities of Weibull distribution with r =2 (blue), r =1 (red), and r =2 (black). In

all three cases,

λ =1/2.

5.5.10 Pareto Distribution

The Pareto distribution is named after the Italian economist Vilfredo Pareto.

Some examples in which the Pareto distribution provides an exemplary model

include wealth distribution in individuals, sizes of human settlements, visits

to encyclopedia pages, and ﬁle size distribution of Internet trafﬁc that uses

the TCP protocol. A random variable X has a Pareto

P a(c,α) distribution

with parameters 0

< c <∞ and α >0 if its density is given by

(x) =

α+1

, x ≥ c,

0, else .

The CDF is

(x) =

0, x < c,

−

, x ≥ c.

172 5 Random Variables

The mean and variance of X are EX = αc/(α −1) and Var X = αc

/[(α −

(α −2)]. The median is m = c ·2

1/α

If X

,... , X

are independent P a(c, α), then Y = 2c

ln(X

) ∼ χ

with

2n degrees of freedom.

In MATLAB one can specify the generalized Pareto distribution, which for

some selection of its parameters is equivalent to the aforementioned Pareto

distribution. In WinBUGS, the code is

dpar(alpha,c) (note the permuted order

of parameters).

5.5.11 Dirichlet Distribution

The Dirichlet distribution is a multivariate version of the beta distribution in

the same way that the multinomial distribution is a multivariate extension of

the binomial. A random variable X

=(X

,... , X

) with a Dirichlet distribution

of (X

∼D ir(a

,... , a

)) has a PDF of

f (x

,... , x

) =

Γ(A)

Γ(a

)

i=1

−1

where A

, and x =(x

,... , x

) ≥0 is deﬁned on the simplex x

+···+x

=1.

Then

E(X

) =

Var (X

) =

(A −a

)

(A +1)

, and

Cov(X

, X

) =−

(A +1)

The Dirichlet random variable can be generated from gamma random vari-

ables Y

,... , Y

∼ G a(a, b) as X

= Y

, i = 1,..., k, where S

. Obvi-

ously, the marginal distribution of a component X

is B e(a

, A −a

). This is

illustrated in the following MATLAB m-ﬁle that generates random Dirichlet

vectors.

function drand = dirichletrnd(a,n)

% function drand = dirichletrnd(a,n)

% a - vector of parameters 1 x m

% n - number of random realizations

% drand - matrix m x n, each column one realization.

%---------------------------------------------------

a=a(:);

m=size(a,1);

a1=zeros(m,n);

for i = 1:m

a1(i,:) = gamrnd(a(i,1),1,1,n);

end

for i=1:m

drand(i, 1:n )= a1(i, 1:n ) ./ sum(a1);

end

5.6 Random Numbers and Probability Tables 173

5.6 Random Numbers and Probability Tables

In older introductory statistics texts, many backend pages have been devoted

to various statistical tables. Several decades ago, many books of statistical

tables were published. Also, the most respected of statistical journals occa-

sionally published articles consisting of statistical tables.

In 1947 the RAND Corporation published the monograph A Million Ran-

dom Digits with 100,000 Normal Deviates, which at the time was a state-of-

the-art resource for simulation and Monte Carlo methods. The random dig-

its can be found at

http://www.rand.org/monograph

reports/MR1418.

html

These days, much larger tables of random numbers can be produced by

a single line of code, resulting in a set of random numbers that can pass a

battery of stringent randomness tests. With MATLAB and many other widely

available software packages, statistical tables and tables of random numbers

are now obsolete. For example, tables of binomial CDF and PDF for a speciﬁc

n and p can be reproduced by

disp(’binocdf(0:12, 12, 0.7)’);

binocdf(0:12, 12, 0.7)

disp(’binopdf(0:12, 12, 0.7)’);

binopdf(0:12, 12, 0.7)

We will show how to sample and simulate from a few distributions in MAT-

LAB and compare empirical means and variances with their theoretical coun-

terparts. The following annotated MATLAB code simulates from binomial,

Poisson, and geometric distributions and compares theoretical and empirical

means and variances.

%various

simulations.m

simu = binornd(12, 0.7, [1,100000]);

% simu is 10000 observations from Bin(12,0.7)

disp(’simu = binornd(12, 0.7, [1,100000]); 12

0.7 - mean(simu)’);

0.7 - mean(simu) %0.001069

%should be small since the theoretical mean is n

disp(’simu = binornd(12, 0.7, [1,100000]); ...

0.7

0.3 - var(simu)’);

0.7

0.3 - var(simu) %-0.008350

%should be small since the theoretical variance is n

(1-p)

%% Simulations from Poisson(2)

poi = poissrnd(2, [1, 100000]);

disp(’poi = poissrnd(2, [1, 100000]); mean(poi)’);

mean(poi) %1.9976

disp(’poi = poissrnd(2, [1, 100000]); var(poi)’);

var(poi) %2.01501

174 5 Random Variables

%%% Simulations from Geometric(0.2)

geo = geornd(0.2, [1, 100000]);

disp(’geo = geornd(0.2, [1, 100000]); mean(geo)’);

mean(geo) %4.00281

disp(’geo = geornd(0.2, [1, 100000]); var(geo)’);

var(geo) %20.11996

5.7 Transformations of Random Variables*

When a random variable with known density is transformed, the result is

a random variable as well. The question is how to ﬁnd its distribution. The

general theory for distributions of functions of random variables is beyond the

scope of this text, and the reader can ﬁnd comprehensive coverage in Ross

(2010a, b).

We have already seen that for a discrete random variable X , the PMF of a

function Y

= g(X ) is simply the table

g(X )

g(x

) g(x

) ··· g(x

) ···

Prob p

··· p

···

in which only realizations of X are transformed while the probabilities are

kept unchanged.

For continuous random variables the distribution of a function is more com-

plex. In some cases, however, looking at the CDF is sufﬁcient.

In this section we will discuss two topics: (i) how to ﬁnd the distribution for

a transformation of a single continuous random variable and (ii) how to ap-

proximate moments, in particular means and variances, of complex functions

of many random variables.

Suppose that a continuous random variable has a density f

(x) and that

a function g is monotone on the domain of f , with the inverse function

h, h

= g

−1

. Then the random variable Y = g(X ) has a density

(y) = f (h(y))|h

(y)|. (5.9)

If g is not 1–1, but has k 1–1 inverse branches, h

, h

,... , h

, then

(y) =

i=1

f (h

(y))|h

(y)|. (5.10)

An example of a non 1–1 function is g(x)

= x

, for which inverse branches

(y) =

y and h

(y) =−

y are 1–1.

5.7 Transformations of Random Variables* 175

Example 5.20. Let X be a random variable with an exponential E (λ) distribu-

tion, where

λ > 0 is the rate parameter. Find the distribution of the random

variable Y

X .

Here g(x)

x and g

−1

(y) = y

. The Jacobian is |g

−1

(y)

|=2y, y ≥0.

Thus,

(y) =λe

−λy

·2y, y ≥0,λ >0,

which is known as the Rayleigh distribution.

An alternative approach to ﬁnding the distribution of Y is to consider the

CDF:

(y) =P(Y ≤ y) =P(

X ≤ y) =P(X ≤ y

) =1 −e

−λy

since X has the exponential distribution.

The density is now obtained by taking the derivative of F

(y),

(y) =(F

(y))

=2λye

−λy

, y ≥0,λ >0.

The distribution of a function of one or many random variables is an ul-

timate summary. However, the result could be quite messy and sometimes

the distribution lacks a closed form. Moreover, not all facets of the resulting

distribution may be of interest to researchers; sometimes only the mean and

variance are needed.

If X is a random variable with EX =µ and Var X =σ

, then for a function

= g(X ) the following approximation holds:

EY ≈ g(µ)+

(µ)σ

Var Y ≈

(µ)

. (5.11)

If n independent random variables are transformed as Y

g(X

, X

,... , X

), then

EY ≈ g(µ

,µ

,... , µ

) +

i=1

∂

∂x

(µ

,µ

,... , µ

)σ

Var Y ≈

i=1

∂g

∂x

(µ

,µ

,... , µ

)

, (5.12)

where

=µ

and Var X

=σ