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

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

6.8 Exercises 227

The original data distinguished ﬁve different strains of trypanosomes, but it

seems that the summary data set, as shown in the table, can be well approx-

imated by a mixture of two normal distributions, p

N µ

,σ

)+p

N µ

,σ

Using MATLAB’s

gmdistribution.fit identify the means of the two normal

components, as well as their weights in the mixture, p

and p

. Plot the

normalized histogram and superimpose the density of the mixture. Data

can be found in

glossina.mat.

6.18. Stabilizing Variance. In Example 6.12 it was stated that the variance

stabilizing transformations for exponential

E (λ) and binomial B in(n, p)

distributions are g(x)

= log(x) and g(x) = arcsin

, respectively. Prove

these statements.

6.19. From Normal to Lognormal. Derive the density of a lognormal distribu-

tion by transforming X

∼N (0,1) into Y =exp{X }.

6.20. The Square of a Standard Normal. If X

∼ N (0,1), show that Y = X

has a density of

(y) =

2Γ

1/2−1

−y/2

, y ≥0,

which is

with 1 degree of freedom.

MATLAB FILES AND DATA SETS USED IN THIS CHAPTER

http://springer.bme.gatech.edu/Ch6.Norm/

acid.m, aviocompany.m, boxcox.m, ch2itf.m, cltdemo.m, glossina.m,

histn.m, ige.m, meanvarind.m, norm2chi2.m, piston.m, plot2dnormal.m,

plotnct.m, quetelet.m, renner.m, tsetse.m

aplysia.odc

glossina.mat, renner.dat|mat

228 6 Normal Distribution

CHAPTER REFERENCES

Altman, D. G. (1980). Statistics and ethics in medical research: misuse of statistics is uneth-

ical. Br. Med. J., 281, 1182–1184.

Casella, G. and Berger, R. (2002). Statistical Inference. Duxbury Press, Belmont

Crow E. L. and Shimizu K. Eds., (1988). Lognormal Distributions: Theory and Application.

Dekker, New York.

Eggert, R. J. (2005). Engineering Design. Pearson - Prentice Hall, Englewood Cliffs

Jung, J. M., Baek, G. H., Kim, J. H., Lee, Y. H., and Chung, M. S. (2001). Changes in ulnar

variance in relation to forearm rotation and grip. J. Bone Joint Surg. Br., 83, 7, 1029–

1033, PubMed PMID: 11603517.

Koike, H. (1987). The extensibility of Aplysia nerve and the determination of true axon

length. J. Physiol., 390, 469–487.

Pearson, K. (1914). On the probability that two independent distributions of frequency are

really samples of the same population, with special reference to recent work on the

identity of trypanosome strains. Biometrika, 10, 1, 85–143.

Renner E. (1970). Mathematisch-statistische Methoden in der praktischen Anwendung.

Parey, Hamburg.

Sautter, C. (1993). Development of a microtargeting device for particle bombardment of

plant meristems. Plant Cell Tiss. Org., 33, 251–257.

Chapter 7

Point and Interval Estimators

A grade is an inadequate report of an inaccurate judgment by a biased and variable

judge of the extent to which a student has attained an undeﬁned level of mastery of

an unknown proportion of an indeﬁnite amount of material. [Dressel, P. L. (1957).

Facts and fancy in assigning grades. Basic College Quarterly, 2, 6–12]

WHAT IS COVERED IN THIS CHAPTER

• Moment Matching and Maximum Likelihood Estimators

• Unbiased and Consistent Estimators

• Estimation of Mean and Variance

• Conﬁdence Intervals

• Estimation of Population Proportions

• Sample Size Design by Length of Conﬁdence Intervals

• Prediction and Tolerance Intervals

• Intervals for the Poisson Rate

7.1 Introduction

One of the primary objectives of inferential statistics is estimation of popula-

tion characteristics (descriptors) on the basis of limited information contained

in a sample. The population descriptors are formalized by a statistical model,

which can be postulated at various levels of speciﬁcity: a broad class of models,

B. Vidakovic, Statistics for Bioengineering Sciences: With MATLAB and WinBUGS Support, 229

Springer Texts in Statistics, DOI 10.1007/978-1-4614-0394-4_7,

230 7 Point and Interval Estimators

a parametric family, or a fully speciﬁc unique model. Often, a functional or dis-

tributional form is fully speciﬁed but dependent on one or more parameters.

Such a model is called parametric. When the model is parametric, the task of

estimation is to ﬁnd the best possible sample counterparts as estimators for

the parameters and to assess the accuracy of the estimators.

The estimation procedure follows standard rules. Usually, a sample is

taken and a statistic (a function of observations) is calculated. The value of

the statistic serves as a point estimator for the unknown population parame-

ter. For example, responses in political pools observed as sample proportions

are used to estimate the population proportion of voters in favor of a particu-

lar candidate. The associated model is binomial and the parameter of interest

is the binomial proportion in the population.

The estimators for a parameter can be given as a single value – point es-

timators or as interval estimators. For example, the sample mean is a point

estimator of the population mean. Conﬁdence intervals and credible sets in a

Bayesian context are examples of interval estimators.

In this chapter we ﬁrst discuss general methods for ﬁnding estimators and

then focus on estimation of speciﬁc population parameters: means, variances,

proportions, rates, etc. Some estimators are universal, that is, they are not con-

nected with any speciﬁc distribution. Universal estimators are a sample mean

for the population mean and a sample variance for the population variance.

However, for interval estimators and for Bayesian estimators, a knowledge of

sampling distribution is critical.

In Chap. 2 we learned about many sample summaries that are good es-

timators for their population counterparts; these will be discussed further in

this chapter. We have also seen some robust competitors based on order statis-

tics and ranks; these will be discussed further in Chap. 12.

The methods for how to propose an estimator for a population parameter

are discussed next. The methods will use knowledge of the form of popula-

tion distribution or, equivalently, distribution of sample summaries treated as

random variables.

7.2 Moment Matching and Maximum Likelihood

Estimators

We describe two approaches for devising point estimators: moment matching

and maximum likelihood.

Matching Estimation. Moment matching is a natural way to propose

an estimator. Theoretical moments of a random variable X with a density

speciﬁed up to a parameter, f (x

|θ), are functions of that parameter:

= h(θ).

7.2 Moment Matching and Maximum Likelihood Estimators 231

For example, if the measurements have a Poisson distribution P oi(λ), the

second moment

is λ +λ

, which is a function of λ. Here, h(x) = x +x

Suppose we obtained a sample X

, X

,... , X

from f (x|θ). The empirical

counterparts for theoretical moments

are sample moments

i=1

By matching the theoretical and empirical moments, an estimator

θ is found

as a solution of the equation

= h(θ).

For example, for the exponential distribution

E (λ), the ﬁrst theoretical

moment is

EX = 1/λ. An estimator for λ is obtained by solving the moment-

matching equation

X = 1/λ, resulting in

= 1/X . Moment-matching es-

timators are not unique; different theoretical and sample moments can be

matched. In the context of an exponential model, the second theoretical mo-

ment is

=2/λ

, leading to an alternative matching equation,

=2/λ

with the solution

The following simple MATLAB code simulates a sample of size 10

from an ex-

ponential distribution with rate parameter

λ =3 and then calculates moment-

matching estimators based on the ﬁrst two moments.

Y = exprnd(1/3, 10e6, 1);

%parametrization in MATLAB is 1/lambda

1/mean(Y) %matching the first moment

ans = 2.9981

sqrt(2/mean(Y.^2)) %matching the second moment

ans = 2.9984

Example 7.1. Consider a sample from a gamma distribution with parameters

r and

λ. It is known that for X ∼G a(r,λ), E(X ) =

, and Var X =EX

−(EX )

. It is easy to see that

(EX )

−(EX )

and λ =

−(EX )

232 7 Point and Interval Estimators

Thus, the moment matching estimators are

(X )

−(X )

and

−(X )

Matching estimation uses mostly moments, but any other statistic that

is (i) easily calculated from a sample and (ii) whose population counterpart

depends on parameter(s) of interest can be used in matching. For example,

the sample/population quantiles can be used.

Example 7.2. In one study, the highest 5-year survival rate (90%) for women

was for malignant melanoma of the skin. Assume that survival time T has an

exponential distribution with an unknown rate parameter

λ. Using quantiles,

estimate

λ.

From

P(T

>5) =0.90 ⇒ exp{−5 ·λ} =0.90

it follows that

λ =0.0211.

Maximum Likelihood. An alternative method, which uses a functional

form for distributions of measurements, is the maximum likelihood estimation

(MLE).

The MLE was ﬁrst proposed and used by R. A. Fisher (Fig. 7.1) in the

1920s and remains one of the most popular tools in estimation theory and

broader statistical inference. The method can be formulated as an optimiza-

tion problem involving the search for extrema when the model is considered

as a function of parameters.

Suppose that the sample X

,... , X

comes from a population with distribu-

tion f (x

|θ) indexed by θ that could be a scalar or a vector of parameters. Ele-

ments of the sample are independent; thus the joint distribution of X

,... , X

is a product of individual densities:

f (x

,... , x

|θ) =

i=1

f (x

|θ).

When the sample is observed, the joint distribution remains dependent

upon the parameter,

θ|X

,... , X

) =

i=1

f (X

|θ), (7.1)

7.2 Moment Matching and Maximum Likelihood Estimators 233

and, as a function of the parameter, L is called the likelihood. The value

of the parameter

θ that maximizes the likelihood L(θ|X

,... , X

) is the

MLE,

mle

Fig. 7.1 Sir R. A. Fisher (1890–1962).

The problem of ﬁnding the maximum of L and the value

mle

at which L

is maximized is an optimization problem. In some cases the maximum can be

found directly or with the help of the log transformation of L. Other times

the procedure must be iterative and the solution is an approximation. In some

cases, depending on the model and sample size, the maximum is not unique or

does not exist.

In the most common cases, maximizing the logarithm of likelihood, log-

likelihood, is simpler than maximizing the likelihood directly. This is because

the product in L becomes the sum when a logarithm is applied:

`(θ|X

,... , X

) =log L(θ|X

,... , X

) =

i=1

log f (X

|θ),

and ﬁnding an extremum of a sum is simpler. Since the logarithm is a mono-

tonically increasing function, the maxima of L and

` are achieved at the same

value

mle

(see Fig. 7.2 for an illustration).

Analytically,

mle

=argmax

`(θ|X

,... , X

and it can be found as a solution of

∂`(θ|X

,... , X

)

∂θ

=0 subject to

∂

`(θ|X

,... , X

)

∂θ

<0.

234 7 Point and Interval Estimators

1 2 3 4 5

−2

−1.5

−1

−0.5

0.5

1.5

log L

Fig. 7.2 Likelihood and log-likelihood of exponential distribution with rate parameter λ

when the sample X =[0.4,0.3,0.1,0.5] is observed. The MLE is 1/X =3.077.

In simple terms, the MLE makes the ﬁrst derivative (with respect to θ) of

the log-likelihood equal to 0 and the second derivative negative, which is a

condition for a maximum.

As an illustration, consider the MLE of

λ in the exponential model, E (λ).

After X

,... , X

is observed, the likelihood becomes

λ|X

,... , X

) =

i=1

λe

−λX

=λ

exp

(

−

i=1

)

The likelihood L is obtained as a product of densities f (x

|λ) where the argu-

ments x

s are ﬁxed observations X

. The product is taken over all observa-

tions, as in (7.1). We can search for the maximum of L directly, but since it is

a product of two terms involving

λ, it is beneﬁcial to look at the log-likelihood

instead.

The log-likelihood is

`(λ|X

,... , X

) = n logλ −λ

i=1

The equation to be solved is

∂`

∂λ

−

i=1

=0,

and the solution is

mle

= 1/X . The second derivative of the log-

likelihood,

∂

∂λ

=−

, is always negative; thus the solution

mle

maximizes `,

and consequently L. Figure 7.2 shows the likelihood and log-likelihood as func-

tions of

λ. For sample X = [0.4,0.3,0.1,0.5] the maximizing λ is 1/X = 3.0769.

7.2 Moment Matching and Maximum Likelihood Estimators 235

Note that both the likelihood and log-likelihood are maximized at the same

value.

For the alternative parametrization of exponentials via a scale parameter,

as in MATLAB, f (x

|λ) =

−x/λ

, the estimator is, of course,

mle

= X .

An important property of MLE is their invariance property.

Invariance Property of MLEs. Let

mle

be an MLE of θ and let η =

g(θ), where g is an arbitrary function. Then

mle

= g(

mle

) is an MLE of

η.

For example, if the MLE for

λ in the exponential distribution was 1/X , then

for a function of the parameter

η =λ

−sin(λ) the MLE is (1/X )

−sin

1/X

In MATLAB, the function

mle ﬁnds the MLE when inputs are data and

the name of a distribution with a list of options. The normal distribution is

the default. For example,

parhat = mle(data) calculates the MLE for µ and σ

of a normal distribution, evaluated at vector data. One of the outputs is the

conﬁdence interval. For example,

[parhat, parci] = mle(data) returns MLEs

and 95% conﬁdence intervals for the parameters. The conﬁdence intervals,

as interval estimators, will be discussed later in this chapter. The command

[...] = mle(data,’distribution’,dist) computes parameter estimations for

the distribution speciﬁed by

dist. Acceptable strings for dist are as follows:

’beta’ ’bernoulli’ ’binomial’

’discrete uniform’ ’exponential’ ’extreme value’

’gamma’ ’generalized extreme value’ ’generalized pareto’

’geometric’ ’lognormal’ ’negative binomial’

’normal’ ’poisson’ ’rayleigh’

’uniform’ ’weibull’

Example 7.3. MLE of Beta in MATLAB. The following MATLAB commands

show how to estimate parameters a and b in a beta distribution. We will simu-

late a sample of size 1000 from a beta

B e(2, 3) distribution and then from the

sample ﬁnd MLEs of a and b.

a = betarnd( 2, 3,[1, 1000]);

thetahat = mle(a,’distribution’, ’beta’)

%thetahat = 1.9991 3.0267

It is possible to ﬁnd the MLE using MATLAB’s mle command for distribu-

tions that are not on the list. The code is given at the end of Example 7.4in

which moment-matching estimators and MLEs for parameters in a Maxwell

distribution are compared.

236 7 Point and Interval Estimators

Example 7.4. Moment-Matching Estimators and MLEs in a Maxwell

Distribution. The Maxwell distribution models random speeds of molecules

in thermal equilibrium as given by statistical mechanics. A random variable

X with a Maxwell distribution is given by the probability density function

f (x

|θ) =

3/2

−θx

, θ >0, x >0.

Assume that we observed velocities X

,... , X

and want to estimate the

unknown parameter

θ.

The following theoretical moments for the Maxwell distribution are avail-

able: the expectation

EX = 2

πθ

, the second moment EX

= 3/θ, and the

fourth moment

= 15/θ

. To ﬁnd moment-matching estimators for θ the

theoretical moments are “matched” with their empirical counterparts

X , X

, and X

, and the resulting equations are solved with

respect to

θ :

X =2

πθ

⇒

π(X )

i=1

⇒

i=1

⇒

15n

To ﬁnd the MLE of

θ we show that the log-likelihood has the form

logθ−

+ part free of θ. The maximum of the log-likelihood is achieved at

MLE

, which is the same as the moment-matching estimator

Speciﬁcally, if X

= 1.4, X

= 3.1, and X

= 2.5 are observed, the MLE

θ is

MLE

17.82

= 0.5051. The other two moment-matching estimators are

=0.4677 and

=0.5768.

In MATLAB, the Maxwell distribution can be custom-deﬁned using the

function

maxwell = @(x,theta) sqrt(2/pi)

...

theta.^(3/2)

x.^2 .

exp( - theta

x.^2/2);

mle([1.4 3.1 2.5], ’pdf’, maxwell, ’start’, rand)

%ans = 0.5051

In most cases, taking the log of likelihood simpliﬁes ﬁnding the MLE. Here

is an example in which the maximization of likelihood was done without the

use of derivatives.