$V\PSWRWLFDOO\(I¿FLHQW

(VWLPDWLRQRID%LYDULDWH

*DXVVLDQ–:HLEXOO

'LVWULEXWLRQ DQG DQ

,QWURGXFWLRQ WR WKH $VVRFLDWHG 3VHXGR

WUXQFDWHG :HLEXOO

6WHYH 3 9HUULOO

-DPHV : (YDQV

'DYLG ( .UHWVFKPDQQ

&KHULO\Q $ +DW¿HOG

United States

Department of

Agriculture

Forest Service

Forest

Products

Laboratory

Research

Paper

FPL–RP–666

1RUPDOZLWKPHDQ&9 

:HLEXOOZLWKPHGLDQ&9 

Hí

Hí

Hí

Hí

Hí

Hí

Hí

Hí

Hí



    

    
July 2012

9HUULOO6WHYH3(YDQV-DPHV:.UHWVFKPDQQ'DYLG(+DW¿HOG

&KHULO\Q$$V\PSWRWLFDOO\(I¿FLHQW(VWLPDWLRQRID%LYDULDWH

Gaussian–Weibull Distribution and an Introduction to the Associated 

Pseudo-truncated Weibull . Research Paper FPL-RP-666. Madison, 

WI: U.S. Department of Agriculture, Forest Service, Forest Products 

Laboratory. 76 p.

A limited number of free copies of this publication are available to the 

public from the Forest Products Laboratory, One Gifford Pinchot Drive, 

Madison, WI 53726–2398. This publication is also available online at 

www.fpl.fs.fed.us. Laboratory publications are sent to hundreds of libraries 

in the United States and elsewhere.

The Forest Products Laboratory is maintained in cooperation with the 

University of Wisconsin. 

7KHXVHRIWUDGHRU¿UPQDPHVLQWKLVSXEOLFDWLRQLVIRUUHDGHULQIRUPDWLRQ

and does not imply endorsement by the United States Department of 

Agriculture (USDA) of any product or service.

The USDA prohibits discrimination in all its programs and activities on the 

basis of race, color, national origin, age, disability, and where applicable, 

sex, marital status, familial status, parental status, religion, sexual orientation, genetic information, political beliefs, reprisal, or because all or a part 

of an individual’s income is derived from any public assistance program. 

(Not all prohibited bases apply to all programs.) Persons with disabilities 

who require alternative means for communication of program informaWLRQ%UDLOOHODUJHSULQWDXGLRWDSHHWFVKRXOGFRQWDFW86'$¶V7$5*(7

&HQWHUDW±YRLFHDQG7''7R¿OHDFRPSODLQWRIGLVFULPLQDWLRQZULWHWR86'$'LUHFWRU2I¿FHRI&LYLO5LJKWV,QGHSHQGHQFH

$YHQXH6::DVKLQJWRQ'&±RUFDOO±

(voice) or (202) 720–6382 (TDD). USDA is an equal opportunity provider 

and employer.

$EVWUDFW

Two important wood properties are stiffness (modulus 

of elasticity or MOE) and bending strength (modulus of 

rupture or MOR). In the past, MOE has often been modeled as a Gaussian and MOR as a lognormal or a two- or 

three-parameter Weibull. It is well-known that MOE and 

MOR are positively correlated. To model the simultaneous behavior of MOE and MOR for the purposes of wood 

system reliability calculations, we introduce a bivariate 

Gaussian–Weibull distribution and the associated pseudo-truncated Weibull. We use asymptotically efficient 

likelihood methods to obtain an estimator of the parameter 

vector of the bivariate Gaussian--Weibull, and then obtain 

the asymptotic distribution of this estimator.

Keywords: Reliability, modulus of rupture, modulus of 

elasticity, one-step Newton estimator, Gaussian copula, 

normal distribution, Weibull distribution, likelihood 

methods

$FNQRZOHGJPHQW

We thank the three referees for their very valuable reviews.

&RQWHQWV

1 Introduction .................................................................. 1

2 A bivariate Gaussian–Weibull distribution .................. 2

3 Asymptotic distribution of the estimated parameter 

 vector of the bivariate Gaussian–Weibull 

 distribution ................................................................... 3

$WUXQFDWHGELYDULDWH*DXVVLDQ±:HLEXOO

 distribution ................................................................... 

5 The pseudo-truncated Weibull distribution .................. 

6 Summary ...................................................................... 6

7 References .................................................................... 6

$SSHQGL[$2%LYDULDWH*DXVVLDQ±:HLEXOOGHQVLW\ ..... 8

$SSHQGL[%2¥n-Consistent Initial Estimators of the 

 Parameters ................................................................... 9

10 Appendix C—Partial Derivatives of log (f (x, w))....... 18

11 Appendix D—Partial Derivatives of y ........................ 

12 Appendix E1—Expectations of the first partials of 

 the likelihood .............................................................. 28

13 Appendix E2—Expectations of the second partials of 

 the likelihood .............................................................. 28

$SSHQGL[(2¶V ...................... 31

15 Appendix F—Miscellaneous useful integrals ............. 

16 Appendix G—Three bivariate normal results ............. 52

17 Appendix H—Existence (and finiteness) of the 

 expectations involving partial derivatives of y ........... 53

18 Appendix I—Positive definite information matrix ..... 55

19 Appendix J—Lehmann’s condition D ......................... 57

$SSHQGL[.2%HKDYLRURISVHXGRWUXQFDWHG:HLEXOO

DVȡĺ ........................................................................ 63

21 Appendix L—Proof that pseudo-truncated Weibulls are 

 not Weibulls ................................................................ 

 ̧

 ̧

1

·

 ̈

 ̈

©

§ w w u w w

I

I

I

I ( M L T T
Asymptotically Efficient Estimation of a Bivariate Gaussian–Weibull

Distribution and an Introduction to the Associated Pseudo-truncated

Weibull

Steve P. Verrill, Mathematical Statistician

James W. Evans, Mathematical Statistician

David E. Kretschmann, Research General Engineer

Cherilyn A. Hatfield, Statistician

Forest Products Laboratory, Madison, Wisconsin

1 Introduction

Two important wood properties are stiffness (modulus of elasticity or MOE) and bending strength

(modulus of rupture or MOR). In the past, MOE has often been modeled as a Gaussian and MOR

as a lognormal or a two- or three-parameter Weibull. (See, for example, ASTM 2010a, Evans and

Green 1988, and Green and Evans 1988.)

Design engineers must ensure that the loads to which wood systems are subjected rarely exceed

the systems’ strengths. To this end ASTM D 2915 (ASTM 2010a), and ASTM D 143 or ASTM D

1990 (ASTM 2010b,c) describe the manner in which “allowable properties” are assigned to populations of structural lumber. In essence, an allowable strength property is calculated by estimating

a fifth percentile of a population (actually a 95% content, one-sided lower 75% tolerance bound)

and then dividing that value by “duration of load” (aging) and safety factors. The intent is that

the population can only be used in applications in which the load does not exceed the allowable

property. Of course there are stochastic issues associated with variable loads, uncertainty in estimation, and the division of a percentile with no consideration of population variability. Thus, from

a statistician’s perspective, this is not an ideal approach to ensuring reliability of wood systems.

However, it is the currently codified approach.

To apply this approach, one must obtain estimates of the fifth percentiles of MOR distributions.

Currently, one method for obtaining estimates involves fitting a two-parameter Weibull distribution

to a sample of MORs. To obtain this fit, either a maximum likelihood approach or a linear regression

approach based on order statistics is permitted under ASTM D 5457 (ASTM 2010d).

Unfortunately, these methods are often applied to populations that are not really distributed

as two-parameter Weibulls. For example, in the United States, construction grade 2 by 4’s are

often classified into visual categories—select structural, number 1, number 2—or into machine

stress-rated (MSR) grades. In the case of MSR grades, MOE boundaries are selected, MOE is

measured nondestructively, and boards are placed into categories based upon the MOE bins into

which the boards fall. Because MOE and MOR are correlated, bins with higher MOE boundaries

also tend to contain board populations with higher MOR values. The fifth percentiles of these

MOR populations are sometimes estimated by fitting Weibull distributions to these populations.

Statisticians recognize that this poses a problem. Even if the full population of lumber strengths

were distributed as a Weibull, we would not expect that subpopulations formed by visual grades

or MOE binning would continue to be distributed as Weibulls.

In fact, such a subpopulation is not distributed as a Weibull. Instead, if the full joint MOE–

MOR population were distributed as a bivariate Gaussian–Weibull, the subpopulation would be

1
distributed as a “pseudo-truncated Weibull” (PTW). In this paper, we obtain the distribution of

a PTW and show how to obtain estimates of its parameters and its quantiles by fitting a bivariate

Gaussian–Weibull to the full MOE–MOR distribution. To do this, we first define a particular

form of a bivariate Gaussian–Weibull distribution. In Sections 2 and 3 of this paper, we describe

this form and establish that it can be fit by asymptotically efficient likelihood methods in the full

MOE–MOR case. In Sections 4 and 5, we discuss the truncated case and derive the density of a

PTW.

In a subsequent paper, we will describe a Web-based computer program that we have developed

to perform the asymptotically efficient fit. We will use a related program to fit simulated MOE–

MOR data and evaluate the small sample quality of the fits. In a third paper, we will show that

Weibull fits to PTW data can yield poor estimates of probabilities of failure.

As an aside, we remark that the bivariate Gaussian–Weibull distribution has uses other than as a

generator of pseudo-truncated Weibulls. For example, engineers who are interested in simulating the

performance of wood systems must begin with a model for the joint stiffness, strength distribution

of the members of the system. Provided that we are considering the full population, a Gaussian–

Weibull is one possible model for this joint distribution.

Bivariate Gaussian–Weibull distributions have not yet appeared in the literature. However,

Gumbel (1960), Freund (1961), Marshall and Olkin (1967), Block and Basu (1974), Clayton (1978),

Lee (1979), Hougaard (1986), Sarker (1987), Lu and Bhattacharyya (1990), Patra and Dey (1999),

Johnson et al. (1999), and others have previously investigated bivariate Weibulls.

We note that the bivariate Gaussian–Weibull distribution that we investigate in the current

paper is not the only possible bivariate distribution with Gaussian and Weibull marginals. In

essence we begin with a “Gaussian copula”—a bivariate uniform distribution generated by starting

with a bivariate normal distribution and then applying normal cumulative distribution functions to

its marginals. However, there is a large literature on alternative copulas (multivariate distributions

with uniform marginals). See, for example, Nelsen (1999) and Jaworski (2010). (Also see Wang

et al. (2008) for an application of copulas to joint models of tree heights and diameters.) These

alternatives would lead to alternative bivariate Gaussian–Weibulls. Ultimately, the test of the

usefulness of our proposed version of a Gaussian–Weibull for a particular application will depend

on the match between the theoretical distribution and data. Still, we believe that the analysis of our

proposed version in the current paper represents a useful step in the construction and evaluation

of bivariate Gaussian–Weibull distributions.

2 A bivariate Gaussian–Weibull distribution

To generate a bivariate Gaussian–Weibull distribution, we follow Johnson and Kotz (1972). (Taylor

and Bender (1988, 1989) introduced this technique in a lumber context.) That is, let X1, X2 be

distributed as independent N(0,1)’s. Define X = μ + σX1 and Y = ρX1 + 1 − ρ2X2. Then

X is distributed as a N(μ, σ2), Y is distributed as a N(0,1), and their correlation is ρ. Now let

U = Φ(Y ). Then U is a Uniform(0,1) random variable that is correlated with X. Finally, let

W = (− log(1 − U))1/β/γ. Then W is distributed as a Weibull with shape parameter β and scale

parameter 1/γ, and the pair X,W have our joint “bivariate Gaussian–Weibull” distribution. In

this paper, we require that β > 1. Given this generating process, it is straightforward to show (see

Appendix A) that the joint density is given by

f(x, w) = bivnorm(x, y) × weib(w)/φ(y)

2
where φ is the N(0,1) probability density function, Φ is the N(0,1) cumulative distribution function,

y = Φ−1



1 − exp −(γ × w)

β



weib(w) = γββwβ−1 exp 

−(γ × w)

β



and

bivnorm(x, y) = 1

2π

1

σ

1 − ρ2 exp(−arg)

where

arg = (x − μ)

2/σ2 − 2ρ(x − μ)y/σ + y2 

/(2(1 − ρ2))

In Figures 1–9 we provide contour plots of the bivariate Gaussian–Weibull distribution for

coefficients of variation (CV’s) equal to 0.35, 0.25, and 0.15, and generating correlations equal to

0.5, 0.7, and 0.9. The CV of a univariate distribution is its standard deviation divided by its mean.

Note that as the CV declines from 0.35 to 0.25 to 0.15 (as the Weibull shape parameter increases

from 3.13 to 4.54 to 7.91) the density contours become much less elliptical. That is, the distribution

diverges from a bivariate normal. We would expect this as a Weibull is “like a normal” for shape

near 3.6 (skewness equals 0.00056, excess kurtosis equals −0.28), and a Weibull becomes skewed to

the left and leptokurtic as the shape increases.

3 Asymptotic distribution of the estimated parameter vector of

the bivariate Gaussian–Weibull distribution

Now assume that we have have n independent pairs of observations, (x1, w1),...,(xn, wn) from the

bivariate Gaussian–Weibull distribution. Then we have the following theorem.

Theorem 1

√n

⎛

⎜⎜⎜⎜⎝

⎛

⎜⎜⎜⎜⎝

μˆ

σˆ

ρˆ

γˆ

βˆ

⎞

⎟⎟⎟⎟⎠ −

⎛

⎜⎜⎜⎜⎝

μ

σ

ρ

γ

β

⎞

⎟⎟⎟⎟⎠

⎞

⎟⎟⎟⎟⎠

D

→ N(0, I(θ)

−1)

where θ ≡ (μ, σ, ρ, γ, β)T , ˆμ, ˆσ are one-step Newton estimators based on the standard univariate

normal maximum likelihood estimators of the mean and standard deviation of a Gaussian, ˆγ, βˆ are

one-step Newton estimators based on the standard maximum likelihood estimators of 1/scale and

shape for a Weibull, ˆρ is a one-step Newton estimator based on the √n-consistent estimator of ρ

introduced in Appendix B, and the elements of I(θ) are listed in Appendix E2 (and E3).

Proof

The proof is an application of theorem 4.2 of chapter 6 of Lehmann (1983). To invoke Lehmann’s

theorem we must establish a series of conditions.

Lehmann’s conditions (A0)–(A2)

That conditions (A0)–(A2) hold is clear.

Lehmann’s condition (A)

That condition (A) holds (the existence of third partials of the density) is clear.

3
Lehmann’s condition (B)

Lehmann’s condition (B)(8) is established in Appendix E1. Lehmann’s condition (B)(9) is established in Appendices E2 and E3.

Lehmann’s condition (C)

The fact that the information matrix is positive definite is established in Appendix I.

Lehmann’s condition (D)

Lehmann’s condition (D) is established in Appendix J. 

A Web program that obtains the asymptotically efficient likelihood estimates of the five parameters and confidence bounds on these estimates will be described in a subsequent paper.

4 A truncated bivariate Gaussian–Weibull distribution

In wood engineering applications, it is often the case that we do not have data from a full bivariate

Gaussian–Weibull distribution. Instead, we have data from the subpopulation that is formed by

considering lumber whose MOE values lie between two pre-determined limits, cl and cu (that is,

we have machine stress-rated lumber). It is clear that the joint density in this case is

f(x, w)/ (Φ((cu − μ)/σ) − Φ((cl − μ)/σ)) (1)

for x between cl and cu and 0 elsewhere.

5 The pseudo-truncated Weibull distribution

The pseudo-truncated Weibull distribution function at w is given by integrating the truncated

bivariate Gaussian–Weibull density (1) over the region [cl, cu] × [0, w]. That is (from result (11) in

Appendix A)

FPTW(w) = 
 w

0

F1(s) × F2(s)/ (Φ((cu − μ)/σ) − Φ((cl − μ)/σ)) ds (2)

where

F1(s) ≡ γββsβ−1 exp 

−(γs)

β



(3)

and

F2(s) ≡


 cu

cl

1

√2π

1

σ

1 − ρ2 exp 

− ((x − μ)/σ − ρy)

2 /(2(1 − ρ2))

dx (4)

where

y = Φ−1



1 − exp −(γs)

β



4
We have

F2 =


 (cu−μ)/σ

(cl−μ)/σ

1

√2π

1

1 − ρ2 exp 

−



x/1 − ρ2 − ρy/1 − ρ2

2

/2



dx

=


 (cu−μ)/



σ

√1−ρ2



(cl−μ)/



σ

√1−ρ2



1

√2π

exp −



x − ρy/1 − ρ2

2

/2



dx (5)

=


 (cu−μ)/



σ

√1−ρ2



−ρy/√1−ρ2

(cl−μ)/



σ

√1−ρ2



−ρy/√1−ρ2

1

√2π

exp −x2/2



dx

= Φ (cu − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



− Φ



(cl − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



From results (2), (3), and (5), the pseudo-truncated Weibull density is given by

fPTW(w) = d

dwFPTW(w) (6)

= γββwβ−1 exp 

−(γw)

β



×



Φ



(cu − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



− Φ



(cl − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



÷ (Φ((cu − μ)/σ) − Φ((cl − μ)/σ))

where

y = Φ−1



1 − exp −(γw)

β



Thus, as we would expect, for ρ = 0, the pseudo-truncated Weibull density is simply the Weibull

density, γββwβ−1 exp 

−(γw)β

. In Appendix K, we show that as ρ → 1, the pseudo-truncated

Weibull density converges to 0 for w below wl or w above wu where wl is defined by

Φ

cl − μ

σ



= 1 − exp −(γwl)

β



(7)

and wu is defined by

Φ

cu − μ

σ



= 1 − exp −(γwu)

β



(8)

In Appendix K, we also show that for w ∈ (wl, wu), as ρ → 1, the pseudo-truncated Weibull density

converges to

γββwβ−1 exp 

−(γw)

β



/



exp −(γwl)

β



− exp −(γwu)

β



That is, as ρ → 1, the density of a pseudo-truncated Weibull converges to the density of a truncated

Weibull.

Figures 10 and 11 are Weibull probability plots of PTW data. That is, we plot the ordered

data from a PTW sample against the predicted ordered data from the best Weibull fit to the data.

If the data really were Weibull, then the plots would be approximately linear. In Figure 10, the

generating X, Y correlation was 0, so the data actually was Weibull and the plot is approximately

linear. In Figure 11, the generating X, Y correlation was 0.99, so the data was “far from Weibull”

and the plot is quite nonlinear. For both data sets, the Weibull coefficient of variation was 0.25

and cl and cu corresponded to the 0.2 and 0.8 quantiles of the Gaussian distribution.

In Appendix L, we formally establish that pseudo-truncated Weibull distributions are not

Weibull distributions.

5
6 Summary

In the context of wood strength modeling, we have introduced a bivariate Gaussian–Weibull distribution and the associated pseudo-truncated Weibull distribution. In this paper, we obtain the

asymptotic distribution of the estimated parameter vector for a bivariate Gaussian–Weibull distribution. In a subsequent paper, we will discuss a Web-based computer program that implements

this theory to obtain estimates of the parameters of a bivariate Gaussian–Weibull distribution. In

a third paper, we will investigate the question of whether allowable property estimates based on a

Weibull assumption can be poor if the strength population is actually a pseudo-truncated Weibull

population.

7 References

ASTM. 2010a. Standard practice for evaluating allowable properties for grades of structural

lumber. D 2915-10. West Conshohocken, PA: ASTM International.

ASTM. 2010b. Standard test methods for small clear specimens of timber. D 143-94. West

Conshohocken, PA: ASTM International.

ASTM. 2010c. Standard practice for establishing allowable properties for visually-graded dimension lumber from in-grade tests of full-size specimens. D 1990-07. West Conshohocken, PA:

ASTM International.

ASTM. 2010d. Standard specification for computing reference resistence of wood-based materials and structural connections for load and resistance factor design. D 5457-04a. West

Conshohocken, PA: ASTM International.

Block, H.W. and Basu, A.P. (1974). “A Continuous Bivariate Exponential Extension,” Journal

of the American Statistical Association, 69, pages 1031–1037.

Clayton, D.G. (1978). “A Model for Association in Bivariate Life Tables and its Application in

Epidemiological Studies of Familial Tendency in Chronic Disease Incidence,” Biometrika, 65,

pages 141–151.

Evans, J.W. and Green, D.W. (1988). Mechanical Properties of Visually Graded Dimension Lumber: Vol 2. Douglas Fir–Larch. Pub. PB-88-159-397. Springfield, VA: National Technical

Information Service. 476 p.

Freund, J.E. (1961). “A Bivariate Extension of the Exponential Distribution,” Journal of the

American Statistical Association, 56, pages 971–977.

Gordon, R.D. (1941). “Values of Mills’ Ratio of Area to Bounding Ordinate and of the Normal

Probability Integral for Large Values of the Argument,” Annals of Mathematical Statistics,

12, pages 364–366.

Green, D.W. and Evans, J.W. (1988). Mechanical Properties of Visually Graded Dimension Lumber: Vol 4. Southern Pine. Pub. PB-88-159-413. Springfield, VA: National Technical Information Service. 412 p.

Gumbel, E.J. (1960). “Bivariate Exponential Distributions,” Journal of the American Statistical

Association, 55, pages 698–707.

6
Hougaard, P. (1986). “A Class of Multivariate Failure Time Distributions,” Biometrika, 73, pages

671–678.

Jaworski, P. (2010). Copula Theory and Its Applications. Springer, Dordrecht.

Johnson, N.L. and Kotz, S. (1970). Continuous Univariate Distributions – 2. Houghton Mifflin

Company, Boston.

Johnson, N.L. and Kotz, S. (1972). Distributions in Statistics: Continuous Multivariate Distributions. John Wiley and Sons, New York, New York.

Johnson, R.A., Evans, J.W., and Green, D.W. (1999). “Some Bivariate Distributions for Modeling the Strength Properties of Lumber.” Research Paper FPL-RP-575. Madison, WI: U.S.

Department of Agriculture, Forest Service, Forest Products Laboratory. 11 pages.

Kendall, M. and Stuart, A. (1977). The Advanced Theory of Statistics, Volume 1. MacMillan

Publishing Co., Inc., New York.

Lee, L. (1979). “Multivariate Distributions having Weibull Properties,” Journal of Multivariate

Analysis, 9, pages 267–277.

Lehmann, E.L. (1983). Theory of Point Estimation. John Wiley and Sons, New York, New York.

Lu, J.C. and Bhattacharyya, G.K. (1990). “Some New Constructions of Bivariate Weibull Models,” Annals of the Institute of Statistical Mathematics, 42, Number 3, pages 543–559.

Marshall, A.W. and Olkin, I. (1967). “A Multivariate Exponential Distribution,” Journal of the

American Statistical Association, 62, pages 30–44.

Nelsen, R.B. (1999). An Introduction to Copulas. Springer, New York.

Rudin, W. (1987). Real and Complex Analysis. McGraw-Hill, New York, New York.

Patra, K. and Dey, D.K. (1999). “A Multivariate Mixture of Weibull Distributions in Reliability

Modeling,” Statistics & Probability Letters, 45, pages 225–235.

Sarkar, S.K. (1987). “A Continuous Bivariate Exponential Distribution,” Journal of the American

Statistical Association, 82, pages 667–675.

Taylor, S.E. and Bender, D.A. (1988). “Simulating correlated lumber properties using a modified multivariate normal approach,” Transactions of the American Society of Agricultural

Engineers, 31(1), pages 182–186.

Taylor, S.E. and Bender, D.A. (1989). “A method for simulating multiple correlated lumber

properties,” Forest Products Journal, 39(7/8), pages 71–74.

Verrill, S. and Durst, M. (2005). “The Decline and Fall of Type II Error Rates,” The American

Statistician, 59, No. 4, pages 287–291.

Wang, M., Rennolls, K., and Tang, S. (2008). “Bivariate Distribution Modeling of Tree Diameters

and Heights: Dependency Modeling Using Copulas,” Forest Science, 54, No. 3, pages 284–

293.

7
8 Appendix A—Bivariate Gaussian–Weibull density

Let X, Y have a joint bivariate normal distribution with

X ∼ N(μ, σ2)

Y ∼ N(0, 1)

and correlation(X, Y ) = ρ.

Since Y ∼ N(0, 1), we know that Φ(Y ) is distributed as a Uniform(0,1). (Here, Φ denotes the

N(0,1) cumulative distribution function.) Thus, we know that

W ≡ (− log(1 − Φ(Y )))1/β /γ ∼ Weibull(γ,β) (9)

(a two-parameter Weibull distribution with scale parameter 1/γ and shape parameter β).

We then say that X,W have a bivariate Gaussian–Weibull distribution with parameters μ, σ,

ρ, γ, and β.

Using the multivariate form of the change-of-variables theorem (see, for example, Rudin 1987),

we can calculate the joint density function of X,W. First, we invert Equation (9) to obtain

Y = Φ−1



1 − exp −(γ × W)

β



Thus, the transform that takes (x, w) to (x, y) is

T(x, w) =  T1(x, w)

T2(x, w)



=

 x

Φ−1 

1 − exp(−(γ × w)β)





The corresponding Jacobian matrix is

 ∂T1/∂x ∂T1/∂w

∂T2/∂x ∂T2/∂w 

=

 1 0

0 γββwβ−1 exp(−(γ × w)β)/φ 

Φ−1 

1 − exp(−(γ × w)β)

 

and the absolute value of its determinant is

det = γββwβ−1 exp(−(γ × w)

β)/φ 

Φ−1



1 − exp(−(γ × w)

β)



Thus, the Gaussian–Weibull probability density function (pdf) at x, w is

bivnorm(x, y, μ, σ, ρ) × det (10)

where

y = Φ−1



1 − exp −(γ × w)

β



and

bivnorm(x, y, μ, σ, ρ) = 1

2π

×

1

σ

1 − ρ2

× exp(arg)

where

arg = − (x − μ)

2/σ2 − 2ρ(x − μ)y/σ + y2

/(2(1 − ρ2))

= − (x − μ)

2/σ2 − 2ρ(x − μ)y/σ + ρ2y2 + y2 − ρ2y2

/(2(1 − ρ2))

= − ((x − μ)/σ − ρy)

2 /(2(1 − ρ2)) − y2/2

That is, the Gaussian–Weibull pdf at x, w is given by

gaussweib(x, w; μ, σ, ρ, γ, β) ≡ γββwβ−1 exp 

−(γw)

β



(11)

×

1

√2π

1

σ

1 − ρ2 exp 

− ((x − μ)/σ − ρy)

2 /(2(1 − ρ2))

8
9 Appendix B—√n-Consistent Initial Estimators of the Parameters

We will need the following two lemmas in our development. The first of these lemmas provides

a useful fact about the tail behavior of normal distributions. Versions of this fact have appeared

previously in the statistical literature. See, for example, the discussions of “Mills’ ratio” in Kendall

and Stuart (1977) and Johnson and Kotz (1970). The particular form of the fact described in

Lemma 1 is due to Gordon (1941). A simple proof of Lemma 1 is given in Verrill and Durst (2005).

Lemma 1

For x < 0,

x2/(x2 + 1) < Φ(x)/ (φ(x)/(−x)) < 1 (12)

and for x > 0,

x2/(x2 + 1) < (1 − Φ(x))/(φ(x)/x) < 1 (13)

where Φ(x) is the N(0,1) cumulative distribution function and φ(x) is the N(0,1) probability density

function.

Lemma 2

Let W ∼ Weibull(γ,β). Then

E (log(W))2

< ∞

Proof

We have


 ∞

0

(log(w))2γββwβ−1 exp −(γw)

β



dw ≤ γββ


 1

0

(log(w))2wβ−1dw (14)

+


 ∞

0

w2γββwβ−1 exp −(γw)

β



dw

We know from Weibull theory that the second term on the right hand side in (14) is finite. Consider

the first integral on the right hand side.


 1

0

(log(w))2wβ−1dw =


 0

−∞

w2 exp(w(β − 1)) exp(w)dw (15)

=


 0

−∞

w2 exp(wβ)dw

= (1/β)


 ∞

0

w2β exp(−wβ)dw

= 2/β3 (16)



A referee has noted that

E (log(W)) = − (γE + β log(γ)) /β (17)

and

E (log(W))2

= π2

6β2 + (E (log(W)))2 (18)

9
where γE is the Euler-Mascheroni constant. Result (17) follows from


 ∞

0

log(w)γββwβ−1 exp −(γw)

β



dw =


 ∞

0

log s1/β

γββs(β−1)/β exp 

−(γβs)

 1

β s1/β−1ds

=



γβ/β 
 ∞

0

log(s) exp −(γβs)



ds

= 1

β


 ∞

0

log(t) exp(−t)dt − β log(γ)



= (−γE − β log(γ)) /β (19)

To invoke theorem 4.2 of Lehmann (1983) to establish that our final estimators of the parameters

are asymptotically efficient, we first need to establish that our initial estimates of the parameters

are √n-consistent. (ˆan is a √n-consistent estimator of a if √n(ˆan − a) = Op(1). A sequence of

random variables {Xn} is Op(1) if given any 	 > 0, we can find constants M, N such that n>N

implies that Prob(Xn > M) < 	.) As our initial estimators of μ and σ, we take the standard onevariable estimators  ̄x = xi/n and s = (xi − x ̄)2/(n − 1). As our initial estimators of γ and

β, we take the one-variable maximum likelihood estimators, ˆγ and βˆ. Thus, our initial estimators

of μ, σ, γ, and β are √n-consistent. Our initial estimator of ρ is given by

ρˆ ≡ sˆxy/

sxx × sˆyy (20)

where

sˆxy ≡ n

i=1

(xi − x ̄)(ˆyi − y

ˆ ̄)

sxx ≡ n

i=1

(xi − x ̄)

2

sˆyy ≡ n

i=1

(ˆyi − y

ˆ ̄)

2

y

ˆ ̄ ≡ n

i=1

yˆi/n

and

yˆi ≡ g(wi; ˆγ, βˆ) ≡ Φ−1



1 − exp −(ˆγ × wi)

βˆ  (21)

Theorem 2

√n (ˆρ − ρ) = Op(1)

where ˆρ is defined in Equation (20).

Proof

Define

sxy ≡ n

i=1

(xi − x ̄)(yi − y ̄)

syy ≡ n

i=1

(yi − y ̄)

2

y ̄ ≡ n

i=1

yi/n

10
where

yi ≡ g(wi; γ,β) ≡ Φ−1



1 − exp −(γ × wi)

β

 (22)

(The distinction between the “hatted” variables in definitions (21) and the “unhatted” variables in

definitions (22) is that in the hatted case, γ,β are replaced by their estimates ˆγ, βˆ.)

We know that

r ≡ sxy/

sxx × syy

is a √n-consistent estimator of ρ. (That is, we know that √n(r − ρ) = Op(1).) Thus, we will be

done if we can show that

√n (r − ρˆ) = Op(1) (23)

We have

r − ρˆ = sxy/

sxx × syy − sˆxy/

sxx × sˆyy

=

n

i=1(xi − x ̄)(yi − y ̄)

√sxx × syy

−

n

i=1(xi − x ̄)(ˆyi − y

ˆ ̄)

√sxx × syy

+

n

i=1(xi − x ̄)(ˆyi − y

ˆ ̄)

√sxx × syy

−

n

i=1(xi − x ̄)(ˆyi − y

ˆ ̄)

sxx × sˆyy

≡ D1 + D2 (24)

To show that √nD1 = Op(1), we need to show that

√n

n

i=1

(xi − x ̄)



yi − y ̄ − (ˆyi − y

ˆ ̄)



/n= Op(1) (25)

By the Cauchy–Schwarz inequality and the fact that n

i=1(xi − x ̄)2/n p

→ σ2, we know that we can

establish result (25) by establishing that

n

i=1



yi − y ̄ − (ˆyi − y

ˆ ̄)

2 = Op(1) (26)

and it is clear that result (26) follows if

n

i=1

(yi − yˆi)

2 = Op(1) (27)

(This follows because n

i=1(zi − z ̄)2 ≤ ni=1 z2

i .)

From definitions (21) and (22) we have

n

i=1

(yi − yˆi)

2 = n

i=1



g(wi; γ,β) − g(wi; ˆγ, βˆ)

2

(28)

By Taylor’s theorem this equals

n

i=1

∂g(wi; θ)

∂γ |θ∗,i (ˆγ − γ) + ∂g(wi; θ)∂β |θ∗,i (βˆ − β)

2

(29)

where θ = (γ,β)T and θ∗,i ≡ (γ∗,i, β∗,i)T lies on the line between (γ,β)T and (ˆγ, βˆ)T .

11
Thus, given the Cauchy–Schwarz inequality, to establish result (27), it is sufficient to establish

n

i=1

∂g(wi; θ)

∂γ |θ∗,i2

(ˆγ − γ)

2 = Op(1) (30)

and

n

i=1

∂g(wi; θ)

∂β |θ∗,i2

(βˆ − β)

2 = Op(1) (31)

Because ˆγ and βˆ are the maximum likelihood estimates of γ and β, to establish results (30) and

(31), it is sufficient to establish

n

i=1

∂g(wi; θ)

∂γ |θ∗,i2

/n = Op(1) (32)

and

n

i=1

∂g(wi; θ)

∂β |θ∗,i2

/n = Op(1) (33)

Consider result (32). We have

n

i=1

∂g(wi; θ)

∂γ |θ∗,i2

/n = 

wi<wlow

∂g(wi; θ)

∂γ |θ∗,i2

/n

+ 

wlow≤wi≤wup

∂g(wi; θ)

∂γ |θ∗,i2

/n (34)

+ 

wup<wi

∂g(wi; θ)

∂γ |θ∗,i2

/n

≡ S1 + S2 + S3

where 0 < wlow < wup. We will show that S1, S2, and S3 are Op(1).

First, consider S2. Let 	 > 0 be given. Recall that in this paper we require that β > 1. Because

γˆ p

→ γ, βˆ p→ β > 1, n

i=1 wr

i /n p

→ E(Wr), and (from Lemma 2) n

i=1(log(wi))2/n p

→ E((log(W))2),

given any γlow, γup such that 0 < γlow <γ<γup, γlow < 1 < γup and any βlow, βup such that

1 < βlow <β<βup, we can find an N such that n>N implies

prob (An ∩ Bn ∩ Cn) > 1 − 	 (35)

where An is the set on which

βlow < β<β ˆ up and γlow < γ<γ ˆ up

Bn is the set on which











n

i=1

w

2βup

i /n − E



W2βup











< 1

and Cn is the set on which











n

i=1

(log(wi))2/n − E (log(W))2











< 1

12
Define

w1/4 ≡ (log(4/3))1/β/γ (36)

w3/4 ≡ (log(4))1/β/γ (37)

so

1 − exp −(γ × w1/4)

β



= 1/4 (38)

1 − exp −(γ × w3/4)

β



= 3/4 (39)

Now choose wlow small enough so that

γlow × wlow < γup × wlow < 1, (γup × wlow)

βlow < min 

1/4,(γ × w1/4)

β



(40)

and choose wup large enough so that

wup > 1, γup × wup > γlow × wup > 1, (γlow × wup)

βlow > (γ × w3/4)β (41)

Let

K ≡ min φ



Φ−1



1 − exp(−(γlow × wlow)

βup)

 , φ 

Φ−1



1 − exp(−(γup × wup)

βup)



Then, because γlow × wlow < 1 and γup × wup > 1, on set An we have

φ



Φ−1



1 − exp(−(γ∗,i × wi)

β∗,i )

 ≥ K for all wi in [wlow, wup] (42)

so









∂g(wi; θ)

∂γ |θ∗,i







 = β∗,iγβ∗,i−1

∗,i wβ∗,i

i exp 

−(γ∗,iwi)

β∗,i

/φ Φ−1



1 − exp(−(γ∗,iwi)

β∗,i )



≤ β∗,i × γβ∗,i−1

∗,i × wβ∗,i

i /K (43)

≤ βup × γβup−1

up × wβup

up /K ≡ M1

From results (35) and (43), we know that n>N implies that

prob(S2 < M2

1 ) > 1 − 	

That is,

S2 = Op(1) (44)

Now consider S1. Let 	 > 0 be given. As in the S2 analysis, we can find an N such that n>N

implies that result (35) holds. Then, by our previous choice of wlow (see (40)), on An we have

(γ∗,i × wlow)

β∗,i < (γup × wlow)β∗,i

< (γup × wlow)

βlow (45)

< min 1/4,(γ × w1/4)

β



and thus, for wi < wlow,

0 < (γ∗,i × wi)

β∗,i < (γ∗,i × wlow)β∗,i < min 

1/4,(γ × w1/4)

β



(46)

13
so

Φ−1



1 − exp −(γ∗,i × wi)

β∗,i < Φ−1



1 − exp −(γ × w1/4)

β

 = Φ−1(1/4) < 0 (47)

Thus, on An, for wi < wlow, we can apply Lemma 1 to obtain









∂g(wi; θ)

∂γ |θ∗,i







 = β∗,iγβ∗,i−1

∗,i wβ∗,i

i exp 

−(γ∗,iwi)

β∗,i

/φ Φ−1



1 − exp(−(γ∗,iwi)

β∗,i )



≤ β∗,i × γβ∗,i−1

∗,i × wβ∗,i

i /den (48)

≤ (βup/γlow) × (γ∗,iwi)

β∗,i/den

where

den = −



1 − exp(−(γ∗,iwi)

β∗,i )



×



Φ−1



1 − exp(−(γ∗,iwi)

β∗,i )

 (49)

(We now want to show that den goes to 0 no faster than the (γ∗,iwi)β∗,i in the numerator.)

By results (47)–(49), on An, for wi < wlow,









∂g(wi; θ)

∂γ |θ∗,i









< − (βup/γlow)/Φ−1(1/4)

× (γ∗,iwi)

β∗,i/



1 − exp(−(γ∗,iwi)

β∗,i )



(50)

Now by result (46), on An, for wi < wlow, and x ≡ (γ∗,iwi)β∗,i ,

1 − exp −(γ∗,iwi)

β∗,i

= 1 − 1 − x + x2/2! − x3/3! + x4/4! − ...

= x



1 − x/2! + x2/3! − x3/4! + ...

> x(1 − x/2)

> x(1 − 1/8) = x × 7/8 (51)

or

1/



1 − exp(−(γ∗,iwi)

β∗,i )



< (8/7) × (1/(γ∗,iwi)

β∗,i ) (52)

Thus, by results (50) and (52), on An, for wi < wlow









∂g(wi; θ)

∂γ |θ∗,i









< − (βup/γlow)/Φ−1(1/4)

× (8/7) ≡ M2 (53)

From results (35) and (53), we know that n>N implies that

prob(S1 < M2

2 ) > 1 − 	

That is,

S1 = Op(1) (54)

Now consider S3. Let 	 > 0 be given. As in the S2 analysis, we can find an N such that n>N

implies that result (35) holds. Then, by our previous choice of wup (see (41)), on An, for wi > wup,

we have

(γ∗,i × wi)

β∗,i ≥ (γlow × wup)βlow (55)

≥ (γ × w3/4)

β

so

Φ−1



1 − exp −(γ∗,i × wi)

β∗,i > Φ−1(3/4) > 0 (56)

14
Thus, on An, for wi > wup, we can apply Lemma 1 to obtain









∂g(wi; θ)

∂γ |θ∗,i







 = β∗,iγβ∗,i−1

∗,i wβ∗,i

i exp 

−(γ∗,iwi)

β∗,i

/φ Φ−1



1 − exp(−(γ∗,iwi)

β∗,i )



≤ β∗,i × γβ∗,i−1

∗,i × wβ∗,i

i × exp 

−(γ∗,iwi)

β∗,i

/den (57)

≤ βup × γβup−1

up × wβup

i × exp 

−(γ∗,iwi)

β∗,i

/den

where

den = exp −(γ∗,iwi)

β∗,i

×



Φ−1



1 − exp(−(γ∗,iwi)

β∗,i )

 (58)

By results (56)–(58), on An, for wi > wup,









∂g(wi; θ)

∂γ |θ∗,i









≤ βup × γβup−1

up × wβup

i /Φ−1(3/4) (59)

Thus, on An ∩ Bn,

S3 ≤



βup × γβup−1

up /Φ−1(3/4)2

×n

i=1

w

2βup

i /n (60)

<



βup × γβup−1

up /Φ−1(3/4)2

×



E



W2βup

+ 1≡ M3

(61)

From results (35) and (60), we know that n>N implies that

prob(S3 < M3) > 1 − 	

That is,

S3 = Op(1) (62)

Results (34), (44), (54), and (62) establish result (32). Thus to complete the proof of (27), we

need to establish result (33).

In general, the proof of result (33) is essentially the same as the proof of result (32). However,

there is one distinction that must be addressed. We have

∂g(wi; θ)

∂β |θ∗,i = (γ∗,iwi)

β∗,i × log(γ∗,iwi) × exp 

−(γ∗,iwi)

β∗,i

/φ Φ−1



1 − exp(−(γ∗,iwi)

β∗,i )



(63)

Given this equality, it is clear that the analogues to the results S2 = Op(1) and S3 = Op(1) follow

as they did in the proof of result (32). However, the analogue of S1 = Op(1) requires slightly more

care. In particular, on An, for wi < wlow, we can apply Lemma 1 to obtain









∂g(wi; θ)

∂β |θ∗,i









= (γ∗,iwi)

β∗,i × |log(γ∗,iwi)| × exp 

−(γ∗,iwi)

β∗,i

/φ Φ−1



1 − exp(−(γ∗,iwi)

β∗,i )



≤ (γ∗,iwi)

β∗,i × |log(γ∗,iwi)| × exp 

−(γ∗,iwi)

β∗,i

/den (64)

where

den = −



1 − exp(−(γ∗,iwi)

β∗,i )



×



Φ−1



1 − exp(−(γ∗,iwi)

β∗,i )

 (65)

15
By results (47), (52), (64) and (65), on An, for wi < wlow,









∂g(wi; θ)

∂β |θ∗,i









< − 1/Φ−1(1/4)

× |log(γ∗,iwi)| × (γ∗,iwi)

β∗,i/



1 − exp(−(γ∗,iwi)

β∗,i )



(66)

≤ − (8/7)/Φ−1(1/4)

× |log(γ∗,iwi)|

From results (35) and (66) and the (ubiquitous) Cauchy–Schwarz inequality, we know that n>N

implies that

prob(T1 < M4) > 1 − 	 (67)

where

T1 ≡ 

wi<wlow

∂g(wi; θ)

∂β |θ∗,i2

/n

and

M4 = 

(8/7)/Φ−1(1/4)2×



J2 + 2J E (log(W))2

+ 11/2

+ E (log(W))2

+ 1

where

J ≡ max (| log(γlow)|, | log(γup)|)

That is,

T1 = Op(1) (68)

and result (33) follows.

As noted above, results (32) and (33) establish results (30) and (31), which establish result (27)

which establishes

√n D1 = Op(1) (69)

To complete the proof of the Theorem we now need to show that

√n D2 = Op(1) (70)

To establish (70), we first need to establish a few facts about yi and ˆyi. By the Cauchy–Schwarz

inequality and result (27), we have

√n



y ̄ − y

ˆ ̄



 ≤ √n

n

i=1

|yi − yˆi| /n

≤ √n

n

i=1

(yi − yˆi)

2/n1/2

=

n

i=1

(yi − yˆi)

2

1/2

= Op(1)

Thus,

√n



y

ˆ ̄2 − y ̄2

= √n



y

ˆ ̄ − y ̄

 y

ˆ ̄ +  ̄y



= Op(1) (71)

16
By the Cauchy–Schwarz inequality, we have

n

i=1

(ˆyi + yi)

2/n = n

i=1

(ˆyi − yi + 2yi)

2/n (72)

=











n

i=1

(ˆyi − yi)

2/n + 4n

i=1

(ˆyi − yi)yi/n + 4n

i=1

y2

i /n











≤ n

i=1

(ˆyi − yi)

2/n + 4 n

i=1

(ˆyi − yi)

2/n1/2 n

i=1

y2

i /n1/2

+ 4n

i=1

y2

i /n

By results (27) and (72) and the fact that

n

i=1

y2

i /n p

→ E(Y 2)

we have

n

i=1

(ˆyi + yi)

2/n = Op(1) (73)

By the Cauchy–Schwarz inequality and results (27) and (73) we have

√n











n

i=1



yˆ2

i − y2i



/n









 = √n











n

i=1

(ˆyi − yi) (ˆyi + yi) /n











(74)

≤ √n

n

i=1

(ˆyi − yi)

2 /n1/2 n

i=1

(ˆyi + yi)

2 /n1/2

= Op(1)

By results (71) and (74) we have

√n (ˆsyy/n − syy/n) = √n

n

i=1

yˆ2

i /n − y

ˆ ̄2 −

n

i=1

y2

i /n − y ̄2

 (75)

= √n

n

i=1



yˆ2

i − y2i



/n − y

ˆ ̄2 − y ̄2



= Op(1)

From result (75) we have

√n



sˆyy/n −



syy/n = √n (ˆsyy/n − syy/n) /



sˆyy/n +



syy/n = Op(1) (76)

Now

D2 ≡

n

i=1(xi − x ̄)(ˆyi − y

ˆ ̄)

√sxx × syy

−

n

i=1(xi − x ̄)(ˆyi − y

ˆ ̄)

sxx × sˆyy

(77)

=

n

i=1(xi − x ̄)(ˆyi − y

ˆ ̄)

n

×

sxx × sˆyy/n2 − sxx × syy/n2

sxx × syy × sxx × sˆyy/n4

≡ F1 × F2

17
By the Cauchy–Schwarz inequality and (75)

|F1| ≤ n

i=1

(xi − x ̄)

2/n1/2 n

i=1

(ˆyi − y

ˆ ̄)

2/n1/2

(78)

= sxx/n ×



sˆyy/n p→ σ × 1

By results (76) and (77)

√nF2 =

sxx/n

sxx × syy × sxx × sˆyy/n4

× √n



sˆyy/n −



syy/n = Op(1) (79)

Results (77), (78), and (79) imply that

√n D2 = Op(1) (80)

This completes the proof of the Theorem. 

10 Appendix C—Partial Derivatives of log(f(x, w))

Given result (10) we have

log(f(x, w)) = − log(2π) − log(σ) − log 1 − ρ2



−

 x−μ

σ

2

2(1 − ρ2) +

2ρ

 x−μ

σ



y

2(1 − ρ2) − y22(1 − ρ2) (81)

+ β log(γ) + log(β)+(β − 1) log(w)

− (γw)

β + log √

2π



+ y2/2

Thus, [first partials]

∂ log(f(x, w))

∂μ = x − μσ2(1 − ρ2) − ρyσ(1 − ρ2) (82)

∂ log(f(x, w))

∂σ = −1σ

+ (x − μ)2

σ3(1 − ρ2) − ρ(x − μ)yσ2(1 − ρ2) (83)

∂ log(f(x, w))

∂ρ = ρ1 − ρ2 −

x − μ

σ

2

+ y2



ρ

(1 − ρ2)2 +

x − μ

σ



y

 1 + ρ2

(1 − ρ2)2



(84)

∂ log(f(x, w))

∂γ = ρ

 x−μ

σ

 ∂y

∂γ

1 − ρ2 − y ∂y

∂γ

1 − ρ2 + βγ − wββγβ−1 + y

∂y

∂γ (85)

∂ log(f(x, w))

∂β = ρ

 x−μ

σ

 ∂y

∂β

1 − ρ2 − y ∂y

∂β

1 − ρ2 + log(γ) + 1β + log(w) − (γw)

β log(γw) + y

∂y

∂β (86)

and [second main partials]

18
∂2 log(f(x, w))

∂μ2 = −1σ2(1 − ρ2) (87)

∂2 log(f(x, w))

∂σ2 = 1σ2 − 3(x − μ)2σ4(1 − ρ2) +

2ρ(x − μ)y

σ3(1 − ρ2) (88)

∂2 log(f(x, w))

∂ρ2 = 11 − ρ2 +

2ρ2

(1 − ρ2)2 −

x − μ

σ

2

+ y2

  1

(1 − ρ2)2 +

4ρ2

(1 − ρ2)3



+

x − μ

σ



y

 2ρ

(1 − ρ2)2 + 4(1 + ρ2)ρ(1 − ρ2)3



(89)

∂2 log(f(x, w))

∂γ2 = ρ

 x−μ

σ



1 − ρ2

∂2y

∂γ2 − 11 − ρ2

∂y

∂γ 2

+ y

∂2y

∂γ2



(90)

− β

γ2 − wββ(β − 1)γβ−2 +

∂y

∂γ 2

+ y

∂2y

∂γ2



∂2 log(f(x, w))

∂β2 = ρ

 x−μ

σ



1 − ρ2

∂2y

∂β2 − 11 − ρ2

 ∂y

∂β 2

+ y

∂2y

∂β2



(91)

− 1

β2 − (γw)

β(log(γw))2 +

 ∂y

∂β 2

+ y

∂2y

∂β2



and [second mixed partials, μ]

∂2 log(f(x, w))

∂μ∂σ = −2(x − μ)σ3(1 − ρ2) + ρy1 − ρ2

1

σ2 (92)

∂2 log(f(x, w))

∂μ∂ρ = 2(x − μ)σ2

ρ

(1 − ρ2)2 − y(1 + ρ2)σ(1 − ρ2)2 (93)

∂2 log(f(x, w))

∂μ∂γ = ρ1 − ρ2

∂y

∂γ −1σ



(94)

∂2 log(f(x, w))

∂μ∂β = ρ1 − ρ2

∂y

∂β −1σ



(95)

and [second mixed partials, σ]

∂2 log(f(x, w))

∂σ∂ρ = 2(x − μ)2σ3

ρ

(1 − ρ2)2 − (x − μ)yσ2

(1 + ρ2)

(1 − ρ2)2 (96)

∂2 log(f(x, w))

∂σ∂γ = −(x − μ)σ2 × ρ1 − ρ2 × ∂y∂γ (97)

∂2 log(f(x, w))

∂σ∂β = −(x − μ)σ2 × ρ1 − ρ2 × ∂y∂β (98)

19
and [second mixed partials, ρ]

∂2 log(f(x, w))

∂ρ∂γ =

x − μ

σ

 ∂y

∂γ  11 − ρ2 +

2ρ2

(1 − ρ2)2



− y × ∂y

∂γ ×

2ρ

(1 − ρ2)2 (99)

∂2 log(f(x, w))

∂ρ∂β =

x − μ

σ

 ∂y

∂β  11 − ρ2 +

2ρ2

(1 − ρ2)2



− y × ∂y

∂β ×

2ρ

(1 − ρ2)2 (100)

and [second mixed partials, γ]

∂2 log(f(x, w))

∂γ∂β =

x − μ

σ

 ρ

1 − ρ2

∂2y

∂γ∂β − 11 − ρ2

∂y

∂γ

∂y

∂β + y

∂2y

∂γ∂β 

(101)

+ 1/γ −



wβ log(w)βγβ−1 + wβγβ−1 + wββγβ−1 log(γ)



+ ∂y

∂γ

∂y

∂β + y

∂2y

∂γ∂β

and [third main partials]

1:

∂3 log(f(x, w))

∂μ3 = 0 (102)

2:

∂3 log(f(x, w))

∂σ3 = −2σ3 + 12(x − μ)2σ5(1 − ρ2) − 6ρ(x − μ)yσ4(1 − ρ2) (103)

3:

∂3 log(f(x, w))

∂ρ3 = 6ρ(1 − ρ2)2 +

8ρ3

(1 − ρ2)3 (104)

−

x − μ

σ

2

+ y2

  12ρ

(1 − ρ2)3 +

24ρ3

(1 − ρ2)4



+

x − μ

σ



y

 2

(1 − ρ2)2 +

8ρ2

(1 − ρ2)3 + 4(1 + 3ρ2)(1 − ρ2)3 +

24ρ2(1 + ρ2)

(1 − ρ2)4



4:

∂3 log(f(x, w))

∂γ3 = ρ1 − ρ2

x − μ

σ

 ∂3y

∂γ3 − 11 − ρ2



3

∂y

∂γ

∂2y

∂γ2 + y

∂3y

∂γ3



(105)

+

2β

γ3 − wββ(β − 1)(β − 2)γβ−3 +



3

∂y

∂γ

∂2y

∂γ2 + y

∂3y

∂γ3



5:

∂3 log(f(x, w))

∂β3 = ρ1 − ρ2

x − μ

σ

 ∂3y

∂β3 − 11 − ρ2



3

∂y

∂β

∂2y

∂β2 + y

∂3y

∂β3



(106)

+

2

β3 − (γw)

β (log(γw))3 +



3

∂y

∂β

∂2y

∂β2 + y

∂3y

∂β3



and [third partials, 2 with respect to μ]

20
6:

∂3 log(f(x, w))

∂μ2∂σ = 2σ3(1 − ρ2) (107)

7:

∂3 log(f(x, w))

∂μ2∂ρ = −2ρσ2(1 − ρ2)2 (108)

8:

∂3 log(f(x, w))

∂μ2∂γ = 0 (109)

9:

∂3 log(f(x, w))

∂μ2∂β = 0 (110)

and [third partials, 2 with respect to σ]

10:

∂3 log(f(x, w))

∂σ2∂μ = 6(x − μ)σ4(1 − ρ2) − 2yρσ3(1 − ρ2) (111)

11:

∂3 log(f(x, w))

∂σ2∂ρ = −6(x − μ)2σ4

ρ

(1 − ρ2)2 + 2(x − μ)yσ3

1 + ρ2

(1 − ρ2)2 (112)

12:

∂3 log(f(x, w))

∂σ2∂γ = 2(x − μ)σ3 × ρ1 − ρ2 × ∂y∂γ (113)

13:

∂3 log(f(x, w))

∂σ2∂β = 2(x − μ)σ3 × ρ1 − ρ2 × ∂y∂β (114)

and [third partials, 2 with respect to ρ]

14:

∂3 log(f(x, w))

∂ρ2∂μ = 2(x − μ)σ2

 1

(1 − ρ2)2 +

4ρ2

(1 − ρ2)3



− y

σ

 2ρ

(1 − ρ2)2 +

4ρ(1 + ρ2)

(1 − ρ2)3



(115)

15:

∂3 log(f(x, w))

∂ρ2∂σ = 2(x − μ)2σ3

 1

(1 − ρ2)2 +

4ρ2

(1 − ρ2)3



(116)

− (x − μ)y

σ2

 2ρ

(1 − ρ2)2 +

4ρ(1 + ρ2)

(1 − ρ2)3



16:

∂3 log(f(x, w))

∂ρ2∂γ =

x − μ

σ

 ∂y

∂γ  2ρ(1 − ρ2)2 +

4ρ

(1 − ρ2)2 +

8ρ3

(1 − ρ2)3



(117)

− y

∂y

∂γ  2(1 − ρ2)2 +

8ρ2

(1 − ρ2)3



17:

∂3 log(f(x, w))

∂ρ2∂β =

x − μ

σ

 ∂y

∂β  2ρ(1 − ρ2)2 +

4ρ

(1 − ρ2)2 +

8ρ3

(1 − ρ2)3



(118)

− y

∂y

∂β  2(1 − ρ2)2 +

8ρ2

(1 − ρ2)3



21
and [third partials, 2 with respect to γ]

18:

∂3 log(f(x, w))

∂γ2∂μ = −1σ

× ρ

1 − ρ2 ×

∂2y

∂γ2 (119)

19:

∂3 log(f(x, w))

∂γ2∂σ = −(x − μ)σ2 × ρ1 − ρ2 ×

∂2y

∂γ2 (120)

20:

∂3 log(f(x, w))

∂γ2∂ρ =

x − μ

σ

  1

1 − ρ2 +

2ρ2

(1 − ρ2)2

 ∂2y

∂γ2 − 2ρ(1 − ρ2)2

∂y

∂γ 2

+ y

∂2y

∂γ2



(121)

21:

∂3 log(f(x, w))

∂γ2∂β =

x − μ

σ



× ρ

1 − ρ2 ×

∂3y

∂γ2∂β (122)

− 1

1 − ρ2



2 × ∂y

∂γ ×

∂2y

∂γ∂β + ∂y∂β ×

∂2y

∂γ2 + y ×

∂3y

∂γ2∂β 

− 1

γ2 −



wβ log(w)(β2 − β)γβ−2 + wβ(2β − 1)γβ−2 + wβ(β2 − β)γβ−2 log(γ)



+



2 × ∂y

∂γ ×

∂2y

∂γ∂β + ∂y∂β ×

∂2y

∂γ2 + y ×

∂3y

∂γ2∂β 

and [third partials, 2 with respect to β]

22:

∂3 log(f(x, w))

∂β2∂μ = −1σ

× ρ

1 − ρ2 ×

∂2y

∂β2 (123)

23:

∂3 log(f(x, w))

∂β2∂σ = −(x − μ)σ2 × ρ1 − ρ2 ×

∂2y

∂β2 (124)

24:

∂3 log(f(x, w))

∂β2∂ρ =

x − μ

σ

  1

1 − ρ2 +

2ρ2

(1 − ρ2)2

 ∂2y

∂β2 − 2ρ(1 − ρ2)2

 ∂y

∂β 2

+ y

∂2y

∂β2



(125)

25:

∂3 log(f(x, w))

∂β2∂γ =

x − μ

σ



× ρ

1 − ρ2 ×

∂3y

∂γ∂β2 (126)

− 1

1 − ρ2



2 × ∂y

∂β ×

∂2y

∂γ∂β + ∂y∂γ ×

∂2y

∂β2 + y ×

∂3y

∂γ∂β2



−



wββγβ−1(log(γw))2 + (γw)

β2 log(γw)/γ

+



2 × ∂y

∂β ×

∂2y

∂γ∂β + ∂y∂γ ×

∂2y

∂β2 + y ×

∂3y

∂γ∂β2



and [third partials, μ, σ, ρ]

26:

∂3 log(f(x, w))

∂μ∂σ∂ρ = −4(x − μ)σ3 × ρ(1 − ρ2)2 + yσ2 ×

 1

1 − ρ2 +

2ρ2

(1 − ρ2)2



(127)

22
and [third partials, μ, σ, γ]

27:

∂3 log(f(x, w))

∂μ∂σ∂γ = 1σ2 × ρ1 − ρ2 × ∂y∂γ (128)

and [third partials, μ, σ, β]

28:

∂3 log(f(x, w))

∂μ∂σ∂β = 1σ2 × ρ1 − ρ2 × ∂y∂β (129)

and [third partials, μ, ρ, γ]

29:

∂3 log(f(x, w))

∂μ∂ρ∂γ = −1σ

×

1 + ρ2

(1 − ρ2)2 × ∂y∂γ (130)

and [third partials, μ, ρ, β]

30:

∂3 log(f(x, w))

∂μ∂ρ∂β = −1σ

×

1 + ρ2

(1 − ρ2)2 × ∂y∂β (131)

and [third partials, μ, γ, β]

31:

∂3 log(f(x, w))

∂μ∂γ∂β = −1σ

× ρ

1 − ρ2 ×

∂2y

∂γ∂β (132)

and [third partials, σ, ρ, γ]

32:

∂3 log(f(x, w))

∂σ∂ρ∂γ = −(x − μ)σ2 ×

1 + ρ2

(1 − ρ2)2 × ∂y∂γ (133)

and [third partials, σ, ρ, β]

33:

∂3 log(f(x, w))

∂σ∂ρ∂β = −(x − μ)σ2 ×

1 + ρ2

(1 − ρ2)2 × ∂y∂β (134)

and [third partials, σ, γ, β]

34:

∂3 log(f(x, w))

∂σ∂γ∂β = −(x − μ)σ2 × ρ1 − ρ2 ×

∂2y

∂γ∂β (135)

and [third partials, ρ, γ, β]

35:

∂3 log(f(x, w))

∂ρ∂γ∂β = x − μσ

 1

1 − ρ2 +

2ρ2

(1 − ρ2)2

 ∂2y

∂γ∂β − 2ρ(1 − ρ2)2

∂y

∂γ

∂y

∂β + y

∂2y

∂γ∂β 

(136)

11 Appendix D—Partial Derivatives of y

Recall that

y ≡ Φ−1



1 − exp(−(γw)

β



Thus, we have

∂y

∂γ = exp 

−(γw)β

βγβ−1wβ

φ(y) (137)

= √

2π × βγβ−1 × wβ × exp −(γw)

β



× exp y2/2



23
∂y

∂β = √

2π × (γw)

β log(γw) × exp 

−(γw)

β



× exp y2/2

 (138)

∂2y

∂γ2 = √

2π βwβ



(β − 1)γβ−2 exp 

−(γw)

β



exp y2/2

 (139)

+ γβ−1 exp −(γw)

β



(−1)wββγβ−1 exp 

y2/2



+ γβ−1 exp 

−(γw)

β



exp y2/2



y

∂y

∂γ 

∂2y

∂β2 = √

2π log(γw)



(γw)

β log(γw) exp 

−(γw)

β



exp y2/2

 (140)

+ (γw)

β exp 

−(γw)

β



(−1)(γw)

β log(γw) exp 

y2/2



+ (γw)

β exp 

−(γw)

β



exp y2/2



y

∂y

∂β 

∂2y

∂γ∂β = √

2π



(1/γ)(γw)

β exp 

−(γw)

β



exp y2/2

 (141)

+ log(γw)wββγβ−1 exp 

−(γw)

β



exp y2/2



+ log(γw)(γw)

β exp 

−(γw)

β



(−1)wββγβ−1 exp 

y2/2



+ log(γw)(γw)

β exp 

−(γw)

β



exp y2/2



y

∂y

∂γ 

∂3y

∂γ3 = √

2πβwβ



(β − 1) × T1 − wββ × T2 + T3



(142)

where

T1 = (β − 2)γβ−3 exp 

−(γw)

β



exp y2/2



+ γβ−2 exp 

−(γw)

β



(−1)wββγβ−1 exp 

y2/2



+ γβ−2 exp 

−(γw)

β



exp y2/2



y

∂y

∂γ

T2 = (2β − 2)γ2β−3 exp 

−(γw)

β



exp y2/2



+ γ2β−2 exp 

−(γw)

β



(−1)wββγβ−1 exp 

y2/2



+ γ2β−2 exp 

−(γw)

β



exp y2/2



y

∂y

∂γ

24
T3 = (β − 1)γβ−2 exp 

−(γw)

β



exp y2/2



y

∂y

∂γ

+ γβ−1 exp 

−(γw)

β



(−1)wββγβ−1 exp 

y2/2



y

∂y

∂γ

+ γβ−1 exp 

−(γw)

β



exp y2/2



y2

∂y

∂γ 2

+ γβ−1 exp −(γw)

β



exp y2/2

 ∂y

∂γ 2

+ γβ−1 exp −(γw)

β



exp y2/2



y

∂2y

∂γ2

∂3y

∂β3 = √

2π log(γw) (log(γw) × T1 − log(γw) × T2 + T3) (143)

where

T1 = (γw)

β log(γw) exp 

−(γw)

β



exp y2/2



+ (γw)

β exp 

−(γw)

β



(−1)(γw)

β log(γw) exp 

y2/2



+ (γw)

β exp 

−(γw)

β



exp y2/2



y

∂y

∂β

T2 = (γw)

2β2 log(γw) exp 

−(γw)

β



exp y2/2



+ (γw)

2β exp 

−(γw)

β



(−1)(γw)

β log(γw) exp 

y2/2



+ (γw)

2β exp 

−(γw)

β



exp y2/2



y

∂y

∂β

T3 = (γw)

β log(γw) exp 

−(γw)

β



exp y2/2



y

∂y

∂β

+ (γw)

β exp 

−(γw)

β



(−1)(γw)

β log(γw) exp 

y2/2



y

∂y

∂β

+ (γw)

β exp 

−(γw)

β



exp y2/2



y2

 ∂y

∂β 2

+ (γw)

β exp 

−(γw)

β



exp y2/2

  ∂y

∂β 2

+ (γw)

β exp 

−(γw)

β



exp y2/2



y

∂2y

∂β2

∂3y

∂γ2∂β = √

2π



wβ × T1 + wββ × T2 − wββ × T3 + wβ × T4



(144)

where

T1 = (β − 1)γβ−2 exp 

−(γw)

β



exp y2/2



+ γβ−1 exp 

−(γw)

β



(−1)wββγβ−1 exp 

y2/2



+ γβ−1 exp 

−(γw)

β



exp y2/2



y

∂y

∂γ

25
T2 = (1/γ)γβ−1 exp 

−(γw)

β



exp y2/2



+ log(γw)(β − 1)γβ−2 exp 

−(γw)

β



exp y2/2



+ log(γw)γβ−1 exp 

−(γw)

β



(−1)wββγβ−1 exp 

y2/2



+ log(γw)γβ−1 exp 

−(γw)

β



exp y2/2



y

∂y

∂γ

T3 = (1/γ)(γw)

β exp 

−(γw)

β



γβ−1 exp 

y2/2



+ log(γw)wββγβ−1 exp 

−(γw)

β



γβ−1 exp 

y2/2



+ log(γw)(γw)

β exp 

−(γw)

β



(−1)wββγβ−1γβ−1 exp 

y2/2



+ log(γw)(γw)

β exp 

−(γw)

β



(β − 1)γβ−2 exp 

y2/2



+ log(γw)(γw)

β exp 

−(γw)

β



γβ−1 exp 

y2/2



y

∂y

∂γ

T4 = (1/γ)γβ exp 

−(γw)

β



exp y2/2



y

∂y

∂γ

+ log(γw)βγβ−1 exp 

−(γw)

β



exp y2/2



y

∂y

∂γ

+ log(γw)γβ exp −(γw)

β



(−1)wββγβ−1 exp 

y2/2



y

∂y

∂γ

+ log(γw)γβ exp 

−(γw)

β



exp y2/2



y2

∂y

∂γ 2

+ log(γw)γβ exp −(γw)

β



exp y2/2

 ∂y

∂γ 2

+ log(γw)γβ exp −(γw)

β



exp y2/2



y

∂2y

∂γ2

∂3y

∂γ∂β2 = √

2π (T1 + log(γw) × T2 − log(γw) × T3 + log(γw) × T4) (145)

where

T1 = γβ−1 log(γ)wβ exp 

−(γw)

β



exp y2/2



+ γβ−1wβ log(w) exp 

−(γw)

β



exp y2/2



+ γβ−1wβ exp 

−(γw)

β



(−1)(γw)

β log(γw) exp 

y2/2



+ γβ−1wβ exp 

−(γw)

β



exp y2/2



y

∂y

∂β

26
T2 = wβ log(w)βγβ−1 exp 

−(γw)

β



exp y2/2



+ wβγβ−1 exp 

−(γw)

β



exp y2/2



+ wββγβ−1 log(γ) exp 

−(γw)

β



exp y2/2



+ wββγβ−1 exp 

−(γw)

β



(−1)(γw)

β log(γw) exp 

y2/2



+ wββγβ−1 exp 

−(γw)

β



exp y2/2



y

∂y

∂β

T3 = (γw)

β log(γw) exp 

−(γw)

β



wββγβ−1 exp 

y2/2



+ (γw)

β exp 

−(γw)

β



(−1)(γw)

β log(γw)wββγβ−1 exp 

y2/2



+ (γw)

β exp 

−(γw)

β



wβ log(w)βγβ−1 exp 

y2/2



+ (γw)

β exp 

−(γw)

β



wβγβ−1 exp 

y2/2



+ (γw)

β exp 

−(γw)

β



wββγβ−1 log(γ) exp 

y2/2



+ (γw)

β exp 

−(γw)

β



wββγβ−1 exp 

y2/2



y

∂y

∂β

T4 = (γw)

β log(γw) exp 

−(γw)

β



exp y2/2



y

∂y

∂γ

+ (γw)

β exp 

−(γw)

β



(−1)(γw)

β log(γw) exp 

y2/2



y

∂y

∂γ

+ (γw)

β exp 

−(γw)

β



exp y2/2



y2 ∂y

∂γ

∂y

∂β

+ (γw)

β exp 

−(γw)

β



exp y2/2

 ∂y

∂γ

∂y

∂β

+ (γw)

β exp 

−(γw)

β



exp y2/2



y

∂2y

∂γ∂β

(146)

12 Appendix E1—Expectations of the first partials of the likelihood

Recall that, by construction, x, y are normal random variables with correlation 1, means μ and 0,

respectively, and standard deviations σ and 1, respectively.

Thus, from result (82), we have

E

∂ log(f(x, w))

∂μ 

= 0 − 0 = 0 (147)

From result (83), we have

E

∂ log(f(x, w))

∂σ 

= −1

σ

+

σ2

σ3(1 − ρ2) − ρ2σσ2(1 − ρ2) = −1σ

+

1

σ(1 − ρ2) − ρ2σ(1 − ρ2) = 0 (148)

27
From result (84), we have

E

∂ log(f(x, w))

∂ρ 

= ρ

1 − ρ2 − 2ρ(1 − ρ2)2 + ρ(1 + ρ2)(1 − ρ2)2 = 0 (149)

From results (85), (261), and (277), we have

E

∂ log(f(x, w))

∂γ 

=

ρ2E



y ∂y

∂γ 

1 − ρ2 −

E



y ∂y

∂γ 

1 − ρ2 +0+ E



y

∂y

∂γ 

= 0 (150)

From results (86), (262), and (278), we have

E

∂ log(f(x, w))

∂β 

=

ρ2E



y ∂y

∂β

1 − ρ2 −

E



y ∂y

∂β

1 − ρ2 +0+ E



y

∂y

∂β 

= 0 (151)

13 Appendix E2—Expectations of the second partials of the likelihood

(The existence and finiteness of the expectations involving partial derivatives of y are established

in Appendix H.)

Recall that, by construction, x, y are normal random variables with correlation 1, means μ and

0, respectively, and standard deviations σ and 1, respectively.

Thus, from result (87), we have

E

∂2 log(f(x, w))

∂μ2



= −1

σ2(1 − ρ2) (152)

From result (88), we have

E

∂2 log(f(x, w))

∂σ2



= 1

σ2 − 3σ2(1 − ρ2) +

2ρ2

σ2(1 − ρ2) = −2 + ρ2σ2(1 − ρ2) (153)

From result (89), we have

E

∂2 log(f(x, w))

∂ρ2



= 1

1 − ρ2 +

2ρ2

(1 − ρ2)2 − 2

 1

(1 − ρ2)2 +

4ρ2

(1 − ρ2)3



+ ρ

 2ρ

(1 − ρ2)2 + 4(1 + ρ2)ρ(1 − ρ2)3



= −(1 + ρ2)

(1 − ρ2)2 (154)

From results (90), (256), and (282), we have

E

∂2 log(f(x, w))

∂γ2



= ρ2

1 − ρ2 E



y

∂2y

∂γ2



− ρ2

1 − ρ2 E

∂y

∂γ 2

+ y

∂2y

∂γ2



− β2

γ2

= −ρ2

1 − ρ2 E

∂y

∂γ 2



− β2

γ2 (155)

28
From results (91), (256), and (287), we have

E

∂2 log(f(x, w))

∂β2



= ρ2

1 − ρ2 E



y

∂2y

∂β2



− ρ2

1 − ρ2 E

 ∂y

∂β 2

+ y

∂2y

∂β2



(156)

− 1

β2 − E



(γw)

β(log(γw))2



= −ρ2

1 − ρ2 E

 ∂y

∂β 2



− 1

β2 − E



(γw)

β(log(γw))2



= −ρ2

1 − ρ2 E

 ∂y

∂β 2



− 1

β2

− E



(log(w))2



− 2

β E (log(w))

− 2 log(γ)E (log(w)) − 2 log(γ)

β

− (log(γ))2

From result (92), we have

E

∂2 log(f(x, w))

∂μ∂σ 

= −2E(x − μ)

σ3(1 − ρ2) + ρE(y)1 − ρ2

1

σ2 = 0 (157)

From result (93), we have

E

∂2 log(f(x, w))

∂μ∂ρ 

= 2E(x − μ)

σ2

ρ

(1 − ρ2)2 − E(y)(1 + ρ2)σ(1 − ρ2)2 = 0 (158)

From result (94), we have

E

∂2 log(f(x, w))

∂μ∂γ 

= ρ

1 − ρ2 E

∂y

∂γ  −1σ



(159)

From result (95), we have

E

∂2 log(f(x, w))

∂μ∂β 

= ρ

1 − ρ2 E

 ∂y

∂β  −1σ



(160)

From results (96) and (256), we have

E

∂2 log(f(x, w))

∂σ∂ρ 

= 2ρ

σ(1 − ρ2)2 − ρ(1 + ρ2)σ(1 − ρ2)2

= ρ

σ(1 − ρ2) (161)

From results (97) and (256), we have

E

∂2 log(f(x, w))

∂σ∂γ 

= −ρ2

σ(1 − ρ2)

E



y

∂y

∂γ 

(162)

From results (98) and (256), we have

E

∂2 log(f(x, w))

∂σ∂β 

= −ρ2

σ(1 − ρ2)

E



y

∂y

∂β 

(163)

29
From results (99) and (256), we have

E

∂2 log(f(x, w))

∂ρ∂γ 

=

 ρ

1 − ρ2 +

2ρ3

(1 − ρ2)2 − 2ρ(1 − ρ2)2



E



y

∂y

∂γ 

= −ρ

1 − ρ2 × E



y

∂y

∂γ 

(164)

From results (100) and (256), we have

E

∂2 log(f(x, w))

∂ρ∂β 

=

 ρ

1 − ρ2 +

2ρ3

(1 − ρ2)2 − 2ρ(1 − ρ2)2



E



y

∂y

∂β 

= −ρ

1 − ρ2 × E



y

∂y

∂β 

(165)

From results (101) and (256), we have

E

∂2 log(f(x, w))

∂γ∂β 

= ρ2

1 − ρ2 E



y

∂2y

∂γ∂β 

− ρ2

1 − ρ2 E

∂y

∂γ

∂y

∂β + y

∂2y

∂γ∂β 

+ 1/γ − E



wβ log(w)βγβ−1



− E



wβγβ−1



− E



wββγβ−1 log(γ)



≡ − ρ2

1 − ρ2 E

∂y

∂γ

∂y

∂β 

+ 1/γ − T1 − T2 − T3 (166)

We have

T1 = − log(w)wββγβ−1 exp 

−(γw)

β







∞

0

+


 ∞

0

γβ−1β



βwβ−1 log(w) + wβ−1



exp −(γw)

β



dw

= β

γ


 ∞

0

log(w)γββwβ−1 exp −(γw)

β



dw +

1

γ


 ∞

0

γββwβ−1 exp −(γw)

β



dw

= β

γ

× E (log(w)) + 1

γ (167)

and

T2 = −γβ−1wβ exp 

−(γw)

β







∞

0

+


 ∞

0

γβ−1βwβ−1 exp −(γw)

β



dw

= 1

γ (168)

and

T3 = − log(γ)γβ−1βwβ exp 

−(γw)

β







∞

0

+ β log(γ)

γ


 ∞

0

γββwβ−1 exp −(γw)

β



dw

= β log(γ)

γ (169)

From results (166)–(169), we have

E

∂2 log(f(x, w))

∂γ∂β 

= − ρ2

1 − ρ2 E

∂y

∂γ

∂y

∂β 

− β

γ

× E (log(w))

− 1

γ − β log(γ)γ (170)

30
14 Appendix E3—E

∂f /∂θi

f × ∂f /∂θjf



’s

(The existence and finiteness of the expectations involving partial derivatives of y are established

in Appendix H. The existence and finiteness of E(log(w)) and E (log(w))2

 follow from Lemma

2.)

Recall from (11) that the pdf of the Gaussian-Weibull is given by

f(x, w) ≡ γββwβ−1 exp 

−(γw)

β



(171)

×

1

√2π

×

1

σ

1 − ρ2

× exp −

x − μ

σ − ρy2

/



2



1 − ρ2

where

y = Φ−1



1 − exp −(γ × w)

β



Thus

∂f

∂μ

f = 1σ

 x−μ

σ − ρy

1 − ρ2 (172)

∂f

∂σ

f = −1σ

+

1

σ

 x−μ

σ − ρy

1 − ρ2

x − μ

σ



∂f

∂ρ

f = ρ1 − ρ2 +

 x−μ

σ − ρy

y

1 − ρ2 −

 x−μ

σ − ρy2

ρ

(1 − ρ2)2

∂f

∂γ

f = βγ − wββγβ−1 +

x − μ

σ − ρy ρ1 − ρ2

∂y

∂γ

∂f

∂β

f = log γ +

1

β + log(w) − (γw)

β log(γw) + x − μ

σ − ρy ρ1 − ρ2

∂y

∂β

From (172) and (256), we have

E

⎛

⎝

 ∂f

∂μ

f

2

⎞

⎠ = E

⎛

⎝

 x−μ

σ − ρy 1σ

1 − ρ2

2

⎞

⎠

= 1

σ2(1 − ρ2)2 E

x − μ

σ

2

− 2ρy x − μ

σ



+ ρ2y2



= 1

σ2(1 − ρ2)2



1 − 2ρ2 + ρ2= 1

σ2(1 − ρ2) (173)

31
From (172), (256), (289), and (290), we have

E

⎛

⎝

 ∂f

∂σ

f

2

⎞

⎠ = E

 1

σ2 − 2σ2(1 − ρ2)

x − μ

σ − ρy x − μσ



+

1

σ2(1 − ρ2)2

x − μ

σ − ρy2 x − μσ

2



= 1

σ2 − 2σ2(1 − ρ2)

(1 − ρ2)

+

1

σ2(1 − ρ2)2 E

x − μ

σ

4

− 2ρy x − μ

σ

3

+ ρ2y2

x − μ

σ

2



= −1

σ2 +

1

σ2(1 − ρ2)2



3 − 6ρ2 + ρ2 + 2ρ4

= −1

σ2 +

1

σ2(1 − ρ2)2



3 − 3ρ2 − 2ρ2 + 2ρ4

= −1

σ2 +

1

σ2(1 − ρ2)



1+2 − 2ρ2

= 1

σ2



1 +

1

1 − ρ2



(174)

From (172), (256), and (273)–(276), we have

E

⎛

⎝

 ∂f

∂ρ

f

2

⎞

⎠ = ρ2

(1 − ρ2)2 +

2ρ

(1 − ρ2)2 E

x − μ

σ



y − ρy2



− 2ρ2

(1 − ρ2)3 E

x − μ

σ

2

− 2ρy x − μ

σ



+ ρ2y2



+

1

(1 − ρ2)2 E

x − μ

σ − ρy2

y2



− 2ρ

(1 − ρ2)3 E

x − μ

σ − ρy3

y



+ ρ2

(1 − ρ2)4 E

x − μ

σ − ρy4



= ρ2

(1 − ρ2)2 +

2ρ

(1 − ρ2)2 (ρ − ρ) − 2ρ2(1 − ρ2)3



1 − 2ρ2 + ρ2

+

1

(1 − ρ2)2 (1 − ρ2) − 2ρ(1 − ρ2)3 × 0 + ρ2(1 − ρ2)4 3(1 − ρ2)

2

= 1 + ρ2

(1 − ρ2)2 (175)

32
From (172), we have

E

⎛

⎝

 ∂f

∂γ

f

2

⎞

⎠ = E

β

γ − wββγβ−1

2



+ 2E

β

γ − wββγβ−1

 x − μ

σ − ρy ρ1 − ρ2

∂y

∂γ 

+ E

x − μ

σ − ρy2

ρ2

(1 − ρ2)2

∂y

∂γ 2



≡ T1 + T2 + T3 (176)

T1 =


 ∞

0

β2

γ2 − 2wββ2γβ−2 + w2ββ2γ2β−2



γββwβ−1 exp 

−(γw)

β



dw

= β2

γ2 − 2


 ∞

0

wββ2γβ−2γββwβ−1 exp −(γw)

β



dw

+


 ∞

0

w2ββ2γ2β−2γββwβ−1 exp −(γw)

β



dw

≡ β2

γ2 − 2U1 + U2 (177)

U1 = −wββ2γβ−2 exp 

−(γw)

β







∞

0

+


 ∞

0

γβ−2β3wβ−1 exp −(γw)

β



dw

= β2

γ2


 ∞

0

γββwβ−1 exp −(γw)

β



dw = β2

γ2 (178)

U2 = −w2ββ2γ2β−2 exp 

−(γw)

β







∞

0

+


 ∞

0

2βw2β−1β2γ2β−2 exp −(γw)

β



dw

=


 ∞

0

2β2wβγβ−2γββwβ−1 exp −(γw)

β



dw

= −2β2wβγβ−2 exp 

−(γw)

β







∞

0

+


 ∞

0

2β3wβ−1γβ−2 exp −(γw)

β



dw

= 2β2

γ2 (179)

From results (177)–(179), we have

T1 = β2

γ2 (180)

From result (258), we have

T2 =


 ∞

0

2

β

γ − wββγβ−1



denweib ρ

1 − ρ2

∂y

∂γ 
 ∞

−∞ x − μ

σ − ρy

denin dx dw = 0 (181)

where denin is defined in (254) and denweib is defined in (259).

From result (264), we have

T3 =


 ∞

0

denweib ρ2

(1 − ρ2)2

∂y

∂γ

2 
 ∞

−∞ x − μ

σ − ρy2

denin dx dw = ρ2

1 − ρ2 E

∂y

∂γ 2



(182)

33
From results (176), and (180)–(182), we have

E

⎛

⎝

 ∂f

∂γ

f

2

⎞

⎠ = ρ2

1 − ρ2 E

∂y

∂γ 2



+ β2

γ2 (183)

From (172), we have

E

⎛

⎝

 ∂f

∂β

f

2

⎞

⎠ = E



log γ +

1

β + log(w) − (γw)

β log(γw)

2



+ 2E

log γ +

1

β + log(w) − (γw)

β log(γw)

 x − μ

σ − ρy ρ1 − ρ2

∂y

∂β 

+ E

x − μ

σ − ρy2

ρ2

(1 − ρ2)2

 ∂y

∂β 2



≡ T1 + T2 + T3 (184)

From result (264), we have

T3 =


 ∞

0

denweib ρ2

(1 − ρ2)2

 ∂y

∂β

2 
 ∞

−∞ x − μ

σ − ρy2

denin dx dw = ρ2

1 − ρ2 E

 ∂y

∂β 2



(185)

From result (258), we have

T2 =


 ∞

0

2



log γ +

1

β + log(w) − (γw)

β log(γw)



denweib ρ

1 − ρ2

∂y

∂β 
 ∞

−∞ x − μ

σ − ρy

denin dx dw = 0

(186)

Next, we have

T1 =



log γ +

1

β



E



log γ +

1

β + log(w) − (γw)

β log(γw)



+



log γ +

1

β



E



log(w) − (γw)

β log(γw)



+ E



log(w) − (γw)

β log(γw)

2



≡ U1 + U2 + U3 (187)

By result (281)

U1 = 0 (188)

and

U2 = −



log(γ) + 1

β

2

(189)

Now,

U3 = E



(log(w))2



+ E



(γw)

2β (log(γ))2



+ E



(γw)

2β (log(w))2



− 2E



(γw)

β log(γ) log(w)



− 2E



(γw)

β (log(w))2



+ 2E



(γw)

2β log(γ) log(w)



≡ V1 + V2 + V3 + V4 + V5 + V6 (190)

34
V2 = (log(γ))2



− (γw)

2β exp 

−(γw)

β







∞

0

+


 ∞

0

γ2β2βw2β−1 exp −(γw)

β



dw

= (log(γ))2


 ∞

0

γβ2wβγββwβ−1 exp −(γw)

β



dw

= (log(γ))2



−2γβwβ exp 

−(γw)

β







∞

0

+ 2 
 ∞

0

γββwβ−1 exp −(γw)

β



dw

= 2 (log(γ))2 (191)

V3 = − (γw)

2β (log(w))2 exp 

−(γw)

β







∞

0

+


 ∞

0

γ2β2βw2β−1(log(w))2 exp −(γw)

β



dw

+


 ∞

0

(γw)

2β 2 log(w)

1

w

exp −(γw)

β



dw

=


 ∞

0

γβ2wβ(log(w))2γββwβ−1 exp −(γw)

β



dw

+

1

β


 ∞

0

γβ2wβ log(w)γββwβ−1 exp −(γw)

β



dw

= −2γβwβ(log(w))2 exp 

−(γw)

β







∞

0

+


 ∞

0

2γβ



βwβ−1(log(w))2 + wβ2 log(w)

1

w



exp −(γw)

β



dw

− 1

β γβ2wβ log(w) exp 

−(γw)

β









∞

0

+


 ∞

0

2

β γβ



βwβ−1 log(w) + wβ 1

w



exp −(γw)

β



dw

= 2 
 ∞

0

(log(w))2γββwβ−1 exp −(γw)

β



dw

+

4

β


 ∞

0

log(w)γββwβ−1 exp −(γw)

β



dw

+

2

β


 ∞

0

log(w)γββwβ−1 exp −(γw)

β



dw

+

2

β2


 ∞

0

γββwβ−1 exp −(γw)

β



dw

= 2E (log(w))2

+

6

β E (log(w)) + 2β2 (192)

V4 = −2 log(γ)



−(γw)

β log(w) exp 

−(γw)

β







∞

0

+


 ∞

0

log(w)γββwβ−1 exp −(γw)

β



dw

+

1

β


 ∞

0

γββwβ−1 exp −(γw)

β



dw

= −2 log(γ)



E (log(w)) + 1

β



(193)

35
V5 = −2



−(γw)

β(log(w))2 exp 

−(γw)

β







∞

0

+


 ∞

0

γβ



βwβ−1(log(w))2 + wβ2 log(w)

1

w



exp −(γw)

β



dw

= −2


 ∞

0

(log(w))2γββwβ−1 exp −(γw)

β



dw +

2

β


 ∞

0

log(w)γββwβ−1 exp −(γw)

β



dw

= −2E (log(w))2

− 4

β E (log(w)) (194)

V6 = 2 log(γ)



−(γw)

2β log(w) exp 

−(γw)

β







∞

0

+


 ∞

0

γ2β



2βw2β−1 log(w) + w2β−1



exp −(γw)

β



dw

= 2 log(γ)


 ∞

0

2γβwβ log(w)γββwβ−1 exp −(γw)

β



dw

+


 ∞

0

γβ 1

β wβγββwβ−1 exp 

−(γw)

β



dw

= 2 log(γ)



−2γβwβ log(w) exp 

−(γw)

β







∞

0

+ 2 
 ∞

0

γβ



βwβ−1 log(w) + wβ−1



exp −(γw)

β



dw

− γβ 1

β wβ exp 

−(γw)

β









∞

0

+

1

β


 ∞

0

γββwβ−1 exp −(γw)

β



dw

= 2 log(γ)



2E (log(w)) + 2

β +

1

β



= 4 log(γ)E (log(w)) + 6 log(γ)

β (195)

From results (190)–(195), we have

U3 = E



(log(w))2



+ 2 (log(γ))2

+ 2E (log(w))2

+

6

β E (log(w)) + 2β2

− 2 log(γ)



E (log(w)) + 1

β



− 2E (log(w))2

− 4

β E (log(w))

+ 4 log(γ)E (log(w)) + 6 log(γ)

β

= E



(log(w))2



+

2

β E (log(w)) + 4 log(γ)β

+ 2 log(γ) (E (log(w))) + 2 (log(γ))2 +

2

β2 (196)

36
Results (187)–(189) and (196) yield

T1 = E



(log(w))2



+

2

β E (log(w)) + 2 log(γ)β

+ 2 log(γ) (E (log(w))) + (log(γ))2 +

1

β2 (197)

Results (184)–(186) and (197) yield

E

⎛

⎝

 ∂f

∂β

f

2

⎞

⎠ = ρ2

1 − ρ2 E

 ∂y

∂β 2



+ E



(log(w))2



+

2

β E (log(w)) + 2 log(γ)β

+ 2 log(γ) (E (log(w))) + (log(γ))2 +

1

β2 (198)

From results (172), (258), and (265), we have

E

 ∂f

∂μ

f ×

∂f

∂σ

f



= E



1

σ

 x−μ

σ − ρy

1 − ρ2 ×



−1

σ

+

1

σ

 x−μ

σ − ρy

1 − ρ2

x − μ

σ



= 0 (199)

From results (172), (258), (266), and (267), we have

E

 ∂f

∂μ

f ×

∂f

∂ρ

f



= E



1

σ

 x−μ

σ − ρy

1 − ρ2 ×



ρ

1 − ρ2 +

 x−μ

σ − ρy

y

1 − ρ2 −

 x−μ

σ − ρy2

ρ

(1 − ρ2)2



= 0 (200)

From results (172), (258), and (264), we have

E

 ∂f

∂μ

f ×

∂f

∂γ

f



= E



1

σ

 x−μ

σ − ρy

1 − ρ2 ×

β

γ − wββγβ−1 +

x − μ

σ − ρy ρ1 − ρ2

∂y

∂γ 

= ρ

σ(1 − ρ2)

E

∂y

∂γ 

(201)

From results (172), (258), and (264), we have

E

 ∂f

∂μ

f ×

∂f

∂β

f



= E



1

σ

 x−μ

σ − ρy

1 − ρ2

×



log γ +

1

β + log(w) − (γw)

β log(γw) + x − μ

σ − ρy ρ1 − ρ2

∂y

∂β 

= ρ

σ(1 − ρ2)

E

 ∂y

∂β 

(202)

37
Next, from result (172), we have

E

 ∂f

∂σ

f ×

∂f

∂ρ

f



= E

−1

σ

+

1

σ

 x−μ

σ − ρy

1 − ρ2

x − μ

σ



×



ρ

1 − ρ2 +

 x−μ

σ − ρy

y

1 − ρ2 −

 x−μ

σ − ρy2

ρ

(1 − ρ2)2



= E



−1

σ



ρ

1 − ρ2 +

 x−μ

σ − ρy

y

1 − ρ2 −

 x−μ

σ − ρy2

ρ

(1 − ρ2)2



+ E



1

σ

 x−μ

σ − ρy

1 − ρ2

x − μ

σ

 ρ

1 − ρ2



+ E



1

σ

 x−μ

σ − ρy

1 − ρ2

x − μ

σ

  x−μ

σ − ρy

y

1 − ρ2



− E



1

σ

 x−μ

σ − ρy

1 − ρ2

x − μ

σ

  x−μ

σ − ρy2

ρ

(1 − ρ2)2



≡ T1 + T2 + T3 + T4 (203)

Results (258) and (264) yield

T1 = 0 (204)

The fact that (x − μ)/σ and y are bivariate standard normals with correlation ρ yields

T2 = ρ

σ(1 − ρ2) (205)

By results (289) and (290), we have

T3 = 1

σ(1 − ρ2)2



E

x − μ

σ

3

y



− 2ρE x − μ

σ

2

y2



+ ρ2E

x − μ

σ



y3



= 1

σ(1 − ρ2)2



3ρ − 2ρ(1 + 2ρ2)+3ρ3

= 1

σ(1 − ρ2)2 × ρ(1 − ρ2) = ρσ(1 − ρ2) (206)

By result (268), we have

T4 = −ρ

σ(1 − ρ2)3 E

x − μ

σ − ρy3 x − μσ



= −3ρ

σ(1 − ρ2) (207)

Results (203)–(207) yield

E

 ∂f

∂σ

f ×

∂f

∂ρ

f



= −ρ

σ(1 − ρ2) (208)

38
Next, from result (172), we have

E

 ∂f

∂σ

f ×

∂f

∂γ

f



= E

−1

σ

+

1

σ

 x−μ

σ − ρy

1 − ρ2

x − μ

σ



×

β

γ − wββγβ−1 +

x − μ

σ − ρy ρ1 − ρ2

∂y

∂γ 

= E

−1

σ

β

γ − wββγβ−1



− ρ

σ(1 − ρ2)

E

x − μ

σ − ρy ∂y∂γ 

+ E



1

σ(1 − ρ2)

x − μ

σ

2

− ρ

x − μ

σ



y

 β

γ − wββγβ−1



+ ρ

σ(1 − ρ2)2 E

x − μ

σ − ρy2 x − μσ

 ∂y

∂γ 

≡ T1 + T2 + T3 + T4 (209)

Result (277) yields

T1 = 0 (210)

Result (258) yields

T2 = 0 (211)

Results (256), (263), and (277) yield

T3 = 0 (212)

Result (270) yields

T4 = ρ2

σ(1 − ρ2)

E



y

∂y

∂γ 

(213)

Results (209)–(213) yield

E

 ∂f

∂σ

f ×

∂f

∂γ

f



= ρ2

σ(1 − ρ2)

E



y

∂y

∂γ 

(214)

39
Next, from result (172), we have

E

 ∂f

∂σ

f ×

∂f

∂β

f



= E

−1

σ

+

1

σ

 x−μ

σ − ρy

1 − ρ2

x − μ

σ



×



log γ +

1

β + log(w) − (γw)

β log(γw) + x − μ

σ − ρy ρ1 − ρ2

∂y

∂β 

= E

−1

σ



log γ +

1

β + log(w) − (γw)

β log(γw)



− ρ

σ(1 − ρ2)

E

x − μ

σ − ρy ∂y∂β 

+

1

σ(1 − ρ2)

E

log γ +

1

β + log(w) − (γw)

β log(γw)

 x − μ

σ − ρy x − μσ



+ ρ

σ(1 − ρ2)2 E

x − μ

σ − ρy2 x − μσ

 ∂y

∂β 

≡ T1 + T2 + T3 + T4 (215)

Result (281) yields

T1 = 0 (216)

Result (258) yields

T2 = 0 (217)

Results (256), (263), and (281) yield

T3 = 0 (218)

Result (270) yields

T4 = ρ2

σ(1 − ρ2)

E



y

∂y

∂β 

(219)

Results (215)–(219) yield

E

 ∂f

∂σ

f ×

∂f

∂β

f



= ρ2

σ(1 − ρ2)

E



y

∂y

∂β 

(220)

40
Next, from result (172), we have

E

 ∂f

∂ρ

f ×

∂f

∂γ

f



= E

 ρ

1 − ρ2 +

 x−μ

σ − ρy

y

1 − ρ2 −

 x−μ

σ − ρy2

ρ

(1 − ρ2)2



×

β

γ − wββγβ−1 +

x − μ

σ − ρy ρ1 − ρ2

∂y

∂γ 

= ρ

1 − ρ2 E

β

γ − wββγβ−1 +

x − μ

σ − ρy ρ1 − ρ2

∂y

∂γ 

+ E

β

γ − wββγβ−1

  x−μ

σ − ρy

y

1 − ρ2



+ E



ρy

(1 − ρ2)2

∂y

∂γ x − μσ − ρy2



− ρ

(1 − ρ2)2 E

β

γ − wββγβ−1

 x − μ

σ − ρy2



− ρ2

(1 − ρ2)3 E



∂y

∂γ x − μσ − ρy3



≡ T1 + T2 + T3 + T4 + T5 (221)

Results (258) and (277) yield

T1 = 0 (222)

Result (258) yields

T2 = 0 (223)

Result (264) yields

T3 = ρ

1 − ρ2 E



y

∂y

∂γ 

(224)

Results (264) and (277) yield

T4 = 0 (225)

Result (274) yields

T5 = 0 (226)

Results (221)–(226) yield

E

 ∂f

∂ρ

f ×

∂f

∂γ

f



= ρ

1 − ρ2 E



y

∂y

∂γ 

(227)

41
Next, from result (172), we have

E

 ∂f

∂ρ

f ×

∂f

∂β

f



= E

 ρ

1 − ρ2 +

 x−μ

σ − ρy

y

1 − ρ2 −

 x−μ

σ − ρy2

ρ

(1 − ρ2)2



×



log γ +

1

β + log(w) − (γw)

β log(γw) + x − μ

σ − ρy ρ1 − ρ2

∂y

∂β 

= ρ

1 − ρ2 E



log γ +

1

β + log(w) − (γw)

β log(γw) + x − μ

σ − ρy ρ1 − ρ2

∂y

∂β 

+ E



log γ +

1

β + log(w) − (γw)

β log(γw)

  x−μ

σ − ρy

y

1 − ρ2



+ E



ρy

(1 − ρ2)2

∂y

∂β x − μσ − ρy2



− ρ

(1 − ρ2)2 E



log γ +

1

β + log(w) − (γw)

β log(γw)

 x − μ

σ − ρy2



− ρ2

(1 − ρ2)3 E



∂y

∂β x − μσ − ρy3



≡ T1 + T2 + T3 + T4 + T5 (228)

Results (258) and (281) yield

T1 = 0 (229)

Result (258) yields

T2 = 0 (230)

Result (264) yields

T3 = ρ

1 − ρ2 E



y

∂y

∂β 

(231)

Results (264) and (281) yield

T4 = 0 (232)

Result (274) yields

T5 = 0 (233)

Results (228)–(233) yield

E

 ∂f

∂ρ

f ×

∂f

∂β

f



= ρ

1 − ρ2 E



y

∂y

∂β 

(234)

42
Next, from result (172), we have

E

 ∂f

∂γ

f ×

∂f

∂β

f



= E

β

γ − wββγβ−1 +

x − μ

σ − ρy ρ1 − ρ2

∂y

∂γ 

×



log γ +

1

β + log(w) − (γw)

β log(γw) + x − μ

σ − ρy ρ1 − ρ2

∂y

∂β 

= E

β

γ − wββγβ−1

 x − μ

σ − ρy ρ1 − ρ2

∂y

∂β 

+ E

log γ +

1

β + log(w) − (γw)

β log(γw)

 x − μ

σ − ρy ρ1 − ρ2

∂y

∂γ 

+ ρ2

(1 − ρ2)2 E



∂y

∂γ

∂y

∂β x − μσ − ρy2



+ E

β

γ − wββγβ−1

 log γ +

1

β + log(w) − (γw)

β log(γw)



≡ T1 + T2 + T3 + T4 (235)

Result (258) yields

T1 = 0 (236)

and

T2 = 0 (237)

Result (264) yields

T3 = ρ2

1 − ρ2 E

∂y

∂γ

∂y

∂β 

(238)

Now note that

T4 =


 ∞

0

β

γ − wββγβ−1

 log γ +

1

β + log(w) − (γw)

β log(γw)



γββwβ−1 exp −(γw)

β



dw

= β

γ



log γ +

1

β



−



log γ +

1

β

 
 ∞

0

wββγβ−1γββwβ−1 exp −(γw)

β



dw

+ β

γ


 ∞

0



log(w) − (γw)

β log(γw)



γββwβ−1 exp −(γw)

β



dw

−


 ∞

0

wββγβ−1



log(w) − (γw)

β log(γw)



γββwβ−1 exp −(γw)

β



dw

≡ U1 + U2 + U3 + U4 (239)

We have


 ∞

0

wββγβ−1γββwβ−1 exp −(γw)

β



dw = −wββγβ−1 exp 

−(γw)

β







∞

0

+ β

γ


 ∞

0

γββwβ−1 exp −(γw)

β



dw

= β

γ

43
Thus

U2 = −



log γ +

1

β

 β

γ (240)

By results (279) and (280), we have

U3 = −β

γ



log γ +

1

β



(241)

Next

U4 = −


 ∞

0

wββγβ−1 log(w)γββwβ−1 exp −(γw)

β



dw

+


 ∞

0

wββγβ−1(γw)

β log(γw)γββwβ−1 exp 

−(γw)

β



dw

≡ V1 + V2 (242)

We have

−V1 = −wββγβ−1 log(w) exp 

−(γw)

β







∞

0

+


 ∞

0

βγβ−1



βwβ−1 log(w) + wβ−1



exp −(γw)

β



dw

= β

γ


 ∞

0

log(w)γββwβ−1 exp −(γw)

β



dw

+

1

γ


 ∞

0

γββwβ−1 exp −(γw)

β



dw

= β

γ

E (log(w)) + 1

γ (243)

Next,

V2 =


 ∞

0

γ2β−1w2ββ log(γ)γββwβ−1 exp −(γw)

β



dw

+


 ∞

0

γ2β−1w2ββ log(w)γββwβ−1 exp −(γw)

β



dw

≡ W1 + W2 (244)

We have

W1 = −γ2β−1w2ββ log(γ) exp 

−(γw)

β







∞

0

+ 2β log(γ)


 ∞

0

γβ−1wβγββwβ−1 exp −(γw)

β



dw

= 2β log(γ)



−γβ−1wβ exp 

−(γw)

β







∞

0

+


 ∞

0

γβ−1βwβ−1 exp −(γw)

β



dw

= 2β log(γ)

γ (245)

44
Also,

W2 = −γ2β−1w2ββ log(w) exp 

−(γw)

β







∞

0

+


 ∞

0

γ2β−1β



2βw2β−1 log(w) + w2β−1



exp −(γw)

β



dw

= 2β


 ∞

0

γβ−1wβ log(w)γββwβ−1 exp −(γw)

β



dw

+


 ∞

0

γβ−1wβγββwβ−1 exp −(γw)

β



dw

≡ X1 + X2 (246)

We have

X1 = 2β



−γβ−1wβ log(w) exp −(γw)

β







∞

0

+


 ∞

0

γβ−1



βwβ−1 log(w) + wβ−1



exp −(γw)

β



dw

= 2β



E

1

γ log(w)



+

1

βγ 

= 2 β

γ

E (log(w)) + 1

γ



(247)

and

X2 = −γβ−1wβ exp 

−(γw)

β







∞

0

+


 ∞

0

γβ−1βwβ−1 exp −(γw)

β



dw = 1

γ (248)

Results (246)–(248) yield

W2 = 2 β

γ

E (log(w)) + 1

γ



+

1

γ (249)

Results (244), (245), and (249) yield

V2 = 2 β

γ

E (log(w)) + 1

γ



+

1

γ

+

2β log(γ)

γ (250)

Results (242), (243), and (250) yield

U4 = β

γ

E (log(w)) + 1

γ

+

1

γ

+

2β log(γ)

γ (251)

Results (239), (240), (241), and (251) yield

T4 = −β

γ



log(γ) + 1

β



+ β

γ

E (log(w)) + 1

γ

+

1

γ

+

2β log(γ)

γ

= β

γ

E (log(w)) + 1

γ

+ β log(γ)

γ (252)

Results (235)–(238) and (252) yield

E

 ∂f

∂γ

f ×

∂f

∂β

f



= ρ2

1 − ρ2 E

∂y

∂γ

∂y

∂β 

+ β

γ

E (log(w)) + 1

γ

+ β log(γ)

γ (253)

45
15 Appendix F—Miscellaneous useful integrals

(The existence and finiteness of the expectations involving partial derivatives of y are established

in Appendix H. The existence and finiteness of E(log(w)) and E (log(w))2

 follow from Lemma

2.)

Define

denin ≡ 1

√2π

×

1

σ

1 − ρ2

× exp −

x − μ

σ − ρy2

/



2



1 − ρ2

(254)

Then, because denin is simply the pdf of a Nμ + ρσy, σ2(1 − ρ)2



,


 ∞

−∞

denin dx = 1 (255)

(Note that the demonstrations of some of the results below could also be based on the fact that

denin is the pdf of a Nμ + ρσy, σ2(1 − ρ)2



. However, instead, in most cases we have found it

expedient to explicitly perform some of the integral calculations. Either approach is valid.)

Next


 ∞

−∞

(x − μ)/σ × denin dx =


 ∞

−∞

1

√2π

× (x − μ)/σ

σ

1 − ρ2

× exp −

x − μ

σ − ρy2

/



2



1 − ρ2

dx

=


 ∞

−∞

1

√2π

×

x

1 − ρ2

× exp − (x − ρy)

2 /



2



1 − ρ2

dx

= 1 − ρ2


 ∞

−∞

x

√2π

× exp −



x − ρy/1 − ρ2

2

/2



dx

= 1 − ρ2 × ρy/1 − ρ2 = ρy (256)

and


 ∞

−∞

(x − μ)/σ × y × denin dx = ρy2 (257)

From results (255) and (256) we have


 ∞

−∞ x − μ

σ − ρy

denin dx = ρy − ρy = 0 (258)

Define

denweib ≡ γββwβ−1 exp 

−(γw)

β



(259)

Then from (11)

gaussweib(x, w; μ, σ, ρ, γ, β) = denweib × denin (260)

and by result (256)

E

x − μ

σ

× ∂y

∂γ 

=


 ∞

0

denweib 
 ∞

−∞

x − μ

σ

× ∂y

∂γ × denin dx dw (261)

=


 ∞

0

denweib × ρy

∂y

∂γ dw

= E



ρy

∂y

∂γ 

46
Similarly,

E

x − μ

σ

× ∂y

∂β 

= E



ρy

∂y

∂β 

(262)

Next, as in the development of result (256),


 ∞

−∞ x − μ

σ

2

× denin dx = (1 − ρ2)


 ∞

−∞

x2

√2π

× exp −



x − ρy/1 − ρ2

2

/2



dx

= (1 − ρ2) ×



1 + ρ2y2

1 − ρ2



= 1 − ρ2 + ρ2y2 (263)

From results (256) and (263) we have


 ∞

−∞ x − μ

σ − ρy2

× denin dx =


 ∞

−∞ x − μ

σ

2

− 2ρy x − μ

σ



+ ρ2y2



denin dx

= 1 − ρ2 + ρ2y2 − 2ρ2y2 + ρ2y2 = 1 − ρ2 (264)

By result (288), we have

E

x − μ

σ − ρy2 x − μσ



= E

x − μ

σ

3

− 2ρ

x − μ

σ

2

y + ρ2

x − μ

σ



y2



= 0 (265)

and

E

x − μ

σ − ρy2

y



= E

x − μ

σ

2

y − 2ρ

x − μ

σ



y2 + ρ2y3



= 0 (266)

and

E

x − μ

σ − ρy3



= E

x − μ

σ

3

− 3ρ

x − μ

σ

2

y + 3ρ2

x − μ

σ



y2 − ρ3E 

y3



= 0

(267)

By results (289) and (290), we have

E

x − μ

σ − ρy3 x − μσ



= E

x − μ

σ

4

− 3

x − μ

σ

3

ρy

+ 3 x − μ

σ

2

ρ2y2 −

x − μ

σ



ρ3y3



= 3 − 9ρ2 + 3(1 + 2ρ2)ρ2 − 3ρ4

= 3 − 6ρ2 + 3ρ4 = 3 − 3ρ2 − 3ρ2 + 3ρ4

= 3(1 − ρ2) − 3ρ2(1 − ρ2) = 3(1 − ρ2)

2 (268)

47
Next, as in the development of (256), we have


 ∞

−∞ x − μ

σ

3

denin dx =


 ∞

−∞

1

√2π

×

x3

1 − ρ2

× exp − (x − ρy)

2 /



2



1 − ρ2

dx

= (1 − ρ2)

3/2


 ∞

−∞

x3

√2π

× exp −



x − ρy/1 − ρ2

2

/2



dx

= (1 − ρ2)

3/2


 ∞

−∞ 

x3 +

3ρyx2

1 − ρ2

+

3ρ2y2x

1 − ρ2 + ρ3y3

(1 − ρ2)3/2



φ(x) dx

= (1 − ρ2)

3/2



3ρy

1 − ρ2

+ ρ3y3

(1 − ρ2)3/2



= 3ρ(1 − ρ2)y + ρ3y3 (269)

From results (256), (263), and (269), we have


 ∞

−∞ x − μ

σ − ρy2 x − μσ



× denin dx

=


 ∞

−∞ x − μ

σ

3

− 2ρ

x − μ

σ

2

y + ρ2

x − μ

σ



y2



× denin dx

= 3ρ(1 − ρ2)y + ρ3y3 − 2ρy(1 − ρ2 + ρ2y2) + ρ2y2(ρy)

= ρ(1 − ρ2)y (270)

Next, by result (261),

E

x − μ

σ − ρy ∂y∂γ 

= E

x − μ

σ

∂y

∂γ 

− E



ρy

∂y

∂γ 

= 0 (271)

Similarly,

E

x − μ

σ − ρy ∂y∂β 

= 0 (272)

Next,

E

x − μ

σ − ρy4



=


 ∞

0

denweib 
 ∞

−∞ x − μ

σ − ρy4

denin dx dw

=


 ∞

0

denweib 
 ∞

−∞

(x − ρy)

4

√2π

×

1

1 − ρ2

× exp − (x − ρy)

2 /



2



1 − ρ2

dx dw

=


 ∞

0

denweib 
 ∞

−∞

1 − ρ2 x − ρy4

√2π

× exp −



x − ρy/1 − ρ2

2

/2



dx dw

=


 ∞

0

denweib 
 ∞

−∞

1 − ρ2

x + ρy/1 − ρ2



− ρy4

√2π

× exp −x2/2



dx dw

= 1 − ρ22


 ∞

0

denweib 
 ∞

−∞

x4

√2π

× exp −x2/2



dx dw

= 3 1 − ρ22 (273)

48
Similarly,


 ∞

−∞ x − μ

σ − ρy3

denin dx = (1 − ρ2)

3/2


 ∞

−∞

x3

√2π

× exp −x2/2



dx = 0 (274)

and

E

x − μ

σ − ρy3

y



= 1 − ρ23/2


 ∞

0

denweib 
 ∞

−∞

x3y

√2π

× exp −x2/2



dx dw = 0 (275)

and (the result below can also be obtained by noting that (x − μ)/σ and y are distributed as

standard normals with correlation ρ and then making use of results (289) and (290))

E

x − μ

σ − ρy2

y2



= 1 − ρ2 
 ∞

0

denweib 
 ∞

−∞

x2y2

√2π

× exp −x2/2



dx dw

= 1 − ρ2 
 ∞

0

denweib × y2 dw

= 1 − ρ2 
 ∞

0

γββwβ−1 exp −(γw)

β



×



Φ−1



1 − exp −(γw)

β

2

dw

= 1 − ρ2 
 ∞

0

dFW (w)

dw

× Φ−1 (FW (w))2 dw

= 1 − ρ2 
 1

0

dFW



F −1

W (s)



dw

× Φ−1 (s)

2 dF −1

W (s)

= 1 − ρ2 
 1

0

dFW



F −1

W (s)



dw

× Φ−1 (s)

2 1

dFW



F −1

W (s)



/dw

ds

= 1 − ρ2 
 1

0



Φ−1 (s)

2 ds

= 1 − ρ2 
 ∞

−∞

s2 dΦ(s) = 1 − ρ2 
 ∞

−∞

s2φ(s) ds = 1 − ρ2 (276)

Next

E

β

γ − Wββγβ−1



=


 ∞

0

β

γ − wββγβ−1



γββwβ−1 exp 

−(γw)

β



dw

= β

γ

+



wββγβ−1 exp 

−(γw)

β







∞

0 −


 ∞

0

β2(γw)

β−1 exp 

−(γw)

β



dw

= β

γ − βγ


 ∞

0

γββwβ−1 exp −(γw)

β



dw = 0 (277)

Next, we want to show that

E



log(γ)+1/β + log(W) − (γW)

β log(γW)



= 0 (278)

We have


 ∞

0

− (γw)

β log(γ)



γββwβ−1 exp −(γw)

β

 dw

= (γw)

β log(γ) exp 

−(γw)

β







∞

0 − log(γ)


 ∞

0

γββwβ−1 exp −(γw)

β



dw

= − log(γ) (279)

49
Also,


 ∞

0

− (γw)

β log(w)



γββwβ−1 exp −(γw)

β

 dw

= (γw)

β log(w) exp 

−(γw)

β







∞

0 −


 ∞

0

d

dw



(γw)

β log(w)



exp −(γw)

β



dw

= −


 ∞

0



γββwβ−1 log(w) + γβwβ−1



exp −(γw)

β



dw

= −


 ∞

0

log(w)γββwβ−1 exp −(γw)

β



dw −


 ∞

0

γβwβ−1 exp −(γw)

β



dw

= −


 ∞

0

log(w)γββwβ−1 exp −(γw)

β



dw − 1/β (280)

Results (279) and (280) imply that

E



log(γ)+1/β + log(W) − (γW)

β log(γW)



= log(γ)+1/β + E (log(W)) (281)

− E



(γW)

β log(γ)



− E



(γW)

β log(w)



= log(γ)+1/β + E (log(W))

− log(γ) − E (log(W)) − 1/β = 0

as desired.

Next,


 ∞

0

 β

γ2 + wββ(β − 1)γβ−2



γββwβ−1 exp 

−(γw)

β



dw

= β

γ2 +



−wββ(β − 1)γβ−2 exp −(γw)

β







∞

0

+


 ∞

0

βwβ−1β(β − 1)γβ−2 exp −(γw)

β



dw

= β

γ2 + β(β − 1)γ2


 ∞

0

γββwβ−1 exp −(γw)

β



dw

= β2

γ2 (282)

Next,


 ∞

0

(γw)

β (log(γw))2

γββwβ−1 exp −(γw)

β



dw

=


 ∞

0

γβwβ (log(w))2γββwβ−1 exp −(γw)

β



dw

+


 ∞

0

γβwβ2 log(γ) log(w)γββwβ−1 exp −(γw)

β



dw

+


 ∞

0

γβwβ (log(γ))2γββwβ−1 exp −(γw)

β



dw

≡ T1 + T2 + T3 (283)

50
We have

T1 = −γβwβ (log(w))2

exp −(γw)

β







∞

0

+


 ∞

0

γβ



βwβ−1(log(w))2 + 2 log(w)wβ−1



exp −(γw)

β



dw

= E



(log(w))2



+

2

β E (log(w)) (284)

and

T2 = −2 log(γ)γβwβ log(w) exp 

−(γw)

β







∞

0

+ 2 log(γ)


 ∞

0

γβ



βwβ−1 log(w) + wβ−1



exp −(γw)

β



dw

= 2 log(γ)E (log(w)) + 2 log(γ)

β (285)

and

T3 = − (log(γ))2

γβwβ exp −(γw)

β







∞

0 + (log(γ))2


 ∞

0

γββwβ−1 exp −(γw)

β



dw

= (log(γ))2 (286)

Results (283)–(286) yield


 ∞

0

(γw)

β (log(γw))2

γββwβ−1 exp −(γw)

β



dw = E



(log(w))2



+

2

β E (log(w))

+ 2 log(γ)E (log(w)) + 2 log(γ)

β

+ (log(γ))2 (287)

16 Appendix G—Three bivariate normal results

Assume that

 S

T



∼ N

 0

0



,

 1 ρ

ρ 1



Then

E ST2 =


 ∞

−∞

t

2

√2π

exp −t

2/2

 
 ∞

−∞

1

√2π

s

1 − ρ2 exp 

−(s − ρt)

2/(2(1 − ρ2))

ds dt

=


 ∞

−∞

t

2

√2π

exp −t

2/2



ρt dt = 0 (288)

Similarly,

E S2T2 =


 ∞

−∞

t

2

√2π

exp −t

2/2

 
 ∞

−∞

1

√2π

s2

1 − ρ2 exp 

−(s − ρt)

2/(2(1 − ρ2))

ds dt

=


 ∞

−∞

t

2

√2π

exp −t

2/2



(1 − ρ2 + ρ2t

2) dt =1+2ρ2 (289)

51
and

E ST3 =


 ∞

−∞

t

3

√2π

exp −t

2/2

 
 ∞

−∞

1

√2π

s

1 − ρ2 exp 

−(s − ρt)

2/(2(1 − ρ2))

ds dt

=


 ∞

−∞

t

3

√2π

exp −t

2/2



ρt dt = 3ρ (290)

17 Appendix H—Existence (and finiteness) of the expectations

involving partial derivatives of y

To complete the calculations presented in Appendices E1-E3, we need to prove that a number of

integrals involving y and its partial derivatives are well-defined. The “proofs” are all similar and

are related to arguments made in the proof of Theorem 2. Let wlow and wup be defined as in that

proof. (We assume that results (40) and (41) hold.)

First, consider E( ∂y

∂γ ). From result (137) we have

∂y

∂γ = exp 

−(γw)β

βγβ−1wβ

φ(y) (291)

where

y ≡ Φ−1



1 − exp(−(γw)

β



Our claim, which in essence establishes all the results in this Appendix, is that


 ∞

0



wβ exp 

−(γw)β

φ(y)



γββwβ−1 exp 

−(γw)

β



dw

is well-defined and finite.

We have


 ∞

0



wβ exp 

−(γw)β

φ(y)



γββwβ−1 exp 

−(γw)

β



dw

=


 wlow

0



wβ exp 

−(γw)β

φ(y)



γββwβ−1 exp 

−(γw)

β



dw

+


 wup

wlow 

wβ exp 

−(γw)β

φ(y)



γββwβ−1 exp 

−(γw)

β



dw

+


 ∞

wup 

wβ exp 

−(γw)β

φ(y)



γββwβ−1 exp 

−(γw)

β



dw

≡ T1 + T2 + T3 (292)

It is clear that T2 is well-defined and finite.

Now consider T3. By result (41), for w>wup

(γw)

β > (γwup)β > (γlowwup)β

> (γlowwup)

βlow > (γw3/4)β

52
so

1 − exp −(γw)

β



> 1 − exp −(γw3/4)

β



= 3/4 (293)

By Lemma 1 and (293), we have

T3 <


 ∞

wup

wβ exp −(γw)β

exp (−(γw)β) × Φ−1(3/4) γββwβ−1 exp 

−(γw)

β



dw (294)

which is clearly well-defined and finite.

Now consider T1. By result (40), for w<wlow

(γw)

β < (γwlow)β < (γupwlow)β

< (γupwlow)

βlow < min 

1/4,(γw1/4)

β



(295)

so

1 − exp −(γw)

β



< 1 − exp −(γw1/4)

β



= 1/4 (296)

By Lemma 1 and (296), we have

T1 <


 wlow

0

wβ exp −(γw)β

(1 − exp (−(γw)β)) × (−Φ−1(1/4)) γββwβ−1 exp 

−(γw)

β



dw (297)

Now for w<wlow, by (40) we have (as in result (295))

(γw)

β < (γupwlow)βlow < 1/4 (298)

Thus, for w<wlow and x ≡ (γw)β,

1 − exp(−(γw)

β)=1 − 

1 − x + x2/2! − x3/3! + x4/4! − ...

= x



1 − x/2! + x2/3! − x3/4! + ...

> x(1 − x/2)

> x(1 − 1/8) = x × 7/8 (299)

Results (297) and (299) imply

T1 <


 wlow

0

wβ exp −(γw)β

((γw)β × 7/8) × (−Φ−1(1/4)) γββwβ−1 exp 

−(γw)

β



dw (300)

which is clearly well-defined and finite. Results (291), (292), (294), and (300) imply that E

∂y

∂γ 

is well-defined and finite.

Now consider E(y ∂y

∂γ ). We have


 ∞

0



wβ exp 

−(γw)β

φ(y)



|y|γββwβ−1 exp 

−(γw)

β



dw

=


 wlow

0



wβ exp 

−(γw)β

φ(y)



|y|γββwβ−1 exp 

−(γw)

β



dw

+


 wup

wlow 

wβ exp 

−(γw)β

φ(y)



|y|γββwβ−1 exp 

−(γw)

β



dw

+


 ∞

wup 

wβ exp 

−(γw)β

φ(y)



|y|γββwβ−1 exp 

−(γw)

β



dw

≡ T1 + T2 + T3 (301)

53
It is clear that T2 is well-defined and finite.

Now consider T3. By Lemma 1, we have

T3 <


 ∞

wup

wβ exp −(γw)β

exp (−(γw)β) × y

yγββwβ−1 exp 

−(γw)

β



dw (302)

which is clearly well-defined and finite.

Now consider T1. By Lemma 1, we have

T1 <


 wlow

0

wβ exp −(γw)β

(1 − exp (−(γw)β)) × (−y)

(−y)γββwβ−1 exp 

−(γw)

β



dw (303)

Results (303) and (299) imply

T1 <


 wlow

0

wβ exp −(γw)β

((γw)β × 7/8) γββwβ−1 exp 

−(γw)

β



dw (304)

which is clearly well-defined and finite. Results (301), (302) and (304) imply that E



y ∂y

∂γ 

is

well-defined and finite.

It is clear via an inspection of the remaining partial derivatives in Appendix D that their

expectations can be handled similarly. (We also use the additional fact that for β > 1, log(w)wβ−1

converges to 0 as w decreases to 0.)

18 Appendix I—Positive definite information matrix

To invoke Lehmann’s Theorem 4.2, we need to establish that the information matrix is positive

definite. In Appendices E2 and E3 we establish that

E



−∂2 log(f(x, w))

∂θi∂θj



= E

∂f /∂θi

f × ∂f /∂θjf



(305)

Thus

aT I(θ)a = 

5

i=1



5

j=1

aiajE

∂f /∂θi

f × ∂f /∂θjf



= E

⎛

⎝



5

i=1

ai

∂f /∂θi

f

2

⎞

⎠ ≥ 0 (306)

To complete the proof that I(θ) is positive definite we need to show that



5

i=1

ai

∂f /∂θi

f = 0 a.e. (307)

54
implies a = 0. From result (172) we have



5

i=1

ai

∂f /∂θi

f = a1 ×



1

σ

 x−μ

σ − ρy

1 − ρ2



(308)

+ a2 ×



−1

σ

+

1

σ

 x−μ

σ − ρy

1 − ρ2

x − μ

σ



+ a3 ×



ρ

1 − ρ2 +

 x−μ

σ − ρy

y

1 − ρ2 −

 x−μ

σ − ρy2

ρ

(1 − ρ2)2



+ a4 ×

β

γ − wββγβ−1 +

x − μ

σ − ρy ρ1 − ρ2

∂y

∂γ 

+ a5 ×



log γ +

1

β + log(w) − (γw)

β log(γw) + x − μ

σ − ρy ρ1 − ρ2

∂y

∂β 

From result (137), we have

∂y

∂γ = βγβ−1 × wβ × exp 

−(γw)

β



/φ(y) (309)

From result (138), we have

∂y

∂β = (γw)

β log(γw) × exp 

−(γw)

β



/φ(y) (310)

Recall that

y ≡ Φ−1



1 − exp(−(γw)

β



Now let 	 > 0 be given. Then results (307) through (310) imply that given any w0, we can find an

associated x, w rectangle chosen so that (x − μ)/σ − ρy is small in the rectangle such that









a2 ×

−1

σ



+ a3 ×

 ρ

1 − ρ2



(311)

+ a4 ×

β

γ − wββγβ−1



+ a5 ×



log γ +

1

β + log(w) − (γw)

β log(γw)









< 	/2

for some (x,w) in the rectangle.

A suitable rectangle can be written as [x0 − δ, x0 + δ] × [w0 − δ, w0 + δ] where δ can be made

arbitrarily small, (x0 − μ)/σ − ρy0 = 0, and y0 = Φ−1 

1 − exp(−(γw0)β. By (307), there must be

some (x,w) in the rectangle for which 5

i=1 ai

∂f /∂θi

f = 0.

Taking w0 large enough

|a4 + a5K log(γw)| < 	 (312)

for K fixed and positive and w arbitrarily large. As 	 was arbitrary, this implies that a4 and a5

equal 0.

Now, given results (307) and (308) and the fact that a4 = a5 = 0, given any 	 > 0, we can find

an x, w region of positive measure (chosen so that y is large and (x − μ)/σ is bounded) such that

(taking y large enough)









a3 ×

 −ρ

1 − ρ2 − ρ3(1 − ρ2)2









< 	 (313)

55
This implies that a3 = 0 or ρ = 0. If ρ = 0, then (given that a4 = a5 = 0)



5

i=1

ai

∂f /∂θi

f = a1 ×

 1

σ

x − μ

σ



+ a2 ×



−1

σ

+

1

σ

x − μ

σ

2



+ a3 ×

x − μ

σ

y



(314)

Given results (307) and (314), given any 	 > 0, we can find an x, w region of positive measure

(chosen so that y is large and (x − μ)/σ is bounded above and bounded below away from 0) such

that (taking y large enough)









a3

x − μ

σ









< 	

for arbitrary (x − μ)/σ in the bounded region. Thus, a3 = 0.

Next, given results (307) and (308) and the fact that a3 = a4 = a5 = 0, given any 	 > 0, we can

find an x, w region of positive measure (chosen so that (x − μ)/σ is large and y is bounded) such

that (letting x get large enough)









a2 ×

1

σ(1 − ρ2)









< 	 (315)

This implies that a2 = 0.

Finally, results (307) and (308) and the fact that a2 = a3 = a4 = a5 = 0 imply that a1 = 0, or

a = 0 as needed.

19 Appendix J—Lehmann’s condition D

Lehmann’s condition D requires that there exists an open 5-dimensional set, S, that contains the

point θ0 = (γ0, β0, ρ0, μ0, σ0), and that there exist functions Mijk(x, w) such that









∂3 log(f(x, w; θ))

∂θiθjθk









≤ Mijk(x, w) (316)

for all θ ∈ S, and

Eθ0 (Mijk(x, w)) < ∞ (317)

for all i, j, k.

The third partials of log(f(x, w; θ)) are displayed in Appendix C. There are 35 third partials.

They can all be handled with five techniques. We will present these techniques in detail by using

them to handle some of the 35 partials, and then describe which of them are needed for each of the

remaining 35 third partials.

Let 0 < γlow < γ0 < γup, γlow < 1 < γup, 1 < βlow < β0 < βup, −1 < ρlow < ρ0 < ρup < 1,

μlow < μ0 < μup, and 0 < σlow < σ0 < σup. Then

S ≡ (γlow, γup) × (βlow, βup) × (ρlow, ρup) × (μlow, μup) × (σlow, σup)

is an open neighborhood of θ0 = (γ0, β0, ρ0, μ0, σ0). Let

S ̄ ≡ [γlow, γup] × [βlow, βup] × [ρlow, ρup] × [μlow, μup] × [σlow, σup]

denote the closure of S.

56
Technique 1—portion of a partial is a continuous function of the parameters

alone (not involving w or x)

Because S ̄ is compact, the image of any continuous function on S ̄ is bounded. Thus, in those cases

in which a portion of the third partial is simply a function of γ, β, ρ, μ, σ (and not of x or w) we can

find a constant M that dominates the absolute value of that portion of the partial for all elements

of S.

This technique permits us to fully handle third partials (listed in Appendix C) 1, 6–9, and to

handle those portions of the other third partials that are dependent only on the parameters (and

not on x or w).

Technique 2—portion of a partial is wβ × exp −(γw)β

× y × exp (y2/2)k

We proceed here in a fashion similar to the fashion in which we proceeded in Appendix H. Let wlow

and wup be defined as in the proof of Theorem 2. That is, we assume that results (36)–(41) hold

with γ,β replaced by γ0, β0. We have


 ∞

0











wβ exp 

−(γw)β

y

φ(y)











k

γβ0

0 β0wβ0−1 exp 

−(γ0w)

β0



dw

=


 wlow

0











wβ exp 

−(γw)β

y

φ(y)











k

γβ0

0 β0wβ0−1 exp 

−(γ0w)

β0



dw

+


 wup

wlow











wβ exp 

−(γw)β

y

φ(y)











k

γβ0

0 β0wβ0−1 exp 

−(γ0w)

β0



dw

+


 ∞

wup











wβ exp 

−(γw)β

y

φ(y)











k

γβ0

0 β0wβ0−1 exp 

−(γ0w)

β0



dw

≡ T1 + T2 + T3 (318)

Note that we are now drawing a distinction between general (γ,β) ∈ [γlow, γup]×[βlow, βup] and

the true parmaeter vector (γ0, β0).

It is clear that for (γ, β, ρ, μ, σ) ∈ S ̄, the 

wβ exp 

−(γw)β

y/φ(y)





k in T2 can be dominated

by a constant that does not depend on the parameter values, and that the resulting integral over

(wlow, wup) is finite.

Now consider T3. By Lemma 1 and the facts that wup > 1 and for w>wup,

(γw)

β > (γlowwup)β > (γlowwup)βlow > (γ0w3/4)β0

(assumption (41)), for w>wup we have











wβ exp 

−(γw)β

y

φ(y)











k

<











wβ × exp 

−(γw)β

× y

exp (−(γw)β) × y











k

< wβupk (319)

Clearly, wβupk does not depend upon the parameter values and the resulting integral over (wup, ∞)

is finite.

57
Now consider T1. By Lemma 1, (299), and the facts that γlow < 1 and

(γw)

β < (γupwlow)β < (γupwlow)βlow < min(1/4,(γ0w1/4)β0 )

we have











wβ exp 

−(γw)β

y

φ(y)











k

<











wβ × exp 

−(γw)β

× y

(1 − exp (−(γw)β)) × y











k

<

 wβ

(γw)β × 7/8

k

<



1

γβup

low × 7/8

k

(320)

Clearly, this quantity does not depend upon the parameter values and the resulting integral over

(−∞, wlow) is finite.

This analysis of T1, T2, and T3 establishes that on S,



wβ × exp 

−(γw)

β



× y × exp y2/2

k

can be dominated by a function that does not depend on the parameter values whose expectation

(with respect to the θ0 pdf) is finite.

Techniques 1, 2, and 5 (see below) together permit us to fully handle third partials 18, 22, and

27–31.

Technique 3—handling (x − μ)/σ and ((x − μ)/σ)

2

First note that |x − μ|/σ < ((x − μ)/σ)2 + 1 so we only need to address ((x − μ)/σ)2.

Define

M ≡ max (|μlow|, |μup|)

Then for θ ∈ S

x − μ

σ

2

≤ (|x| + M)2

σ2

low

= |x|

2 + 2|x|M + M2

σ2

low

(321)

Define

y0 ≡ Φ−1



1 − exp −(γ0w)

β0



Then we have


 ∞

−∞

|x| 1

√2π

1

σ0

1 − ρ2

0

exp − ((x − μ0)/σ0 − ρ0y0)

2 /(2(1 − ρ2

0))

dx

=


 ∞

−∞

|x + μ0| 1

√2π

1

σ0

1 − ρ2

0

exp −



x/ σ0



1 − ρ2

0



− ρ0y0/



1 − ρ2

0

2

/2



dx

=


 ∞

−∞









σ0



1 − ρ2

0 x + μ0









1

√2π

exp −



x − ρ0y0/



1 − ρ2

0

2

/2



dx

=


 ∞

−∞









σ0



1 − ρ2

0



x + ρ0y0/



1 − ρ2

0



+ μ0









1

√2π

exp −x2/2



dx

=


 ∞

−∞









σ0



1 − ρ2

0 x + σ0ρ0y0 + μ0









1

√2π

exp −x2/2



dx

≤ σ0



1 − ρ2

0

2

√2π

+ σ0|ρ0||y0| + |μ0| (322)

58
By a similar argument,


 ∞

−∞

|x|

2 1

√2π

1

σ0

1 − ρ2

0

exp − ((x − μ0)/σ0 − ρ0y0)

2 /(2(1 − ρ2

0))

dx

= σ2

0(1 − ρ20) + σ20ρ20y20 + 2μ0σ0ρ0y0 + μ20 (323)

Results (321)–(323) and technique 4 below establish that for θ ∈ S, ((x − μ)/σ)2 (and thus

|x − μ|/σ) can be dominated by a function that does not depend on θ and that is integrable with

respect to f(x, w; θ0).

Technique 4—handling y and y2

Because |y| < y2 + 1, we can focus on y2.

Let wlow and wup be defined as in the proof of Theorem 2. That is, we assume that results

(37)–(41) hold with γ,β replaced by γ0, β0. We have


 ∞

0

y2γβ0

0 β0wβ0−1 exp 

−(γ0w)

β0



dw =


 wlow

0

y2γβ0

0 β0wβ0−1 exp 

−(γ0w)

β0



dw

+


 wup

wlow

y2γβ0

0 β0wβ0−1 exp 

−(γ0w)

β0



dw

+


 ∞

wup

y2γβ0

0 β0wβ0−1 exp 

−(γ0w)

β0



dw

≡ T1 + T2 + T3 (324)

It is clear that for (γ, β, ρ, μ, σ) ∈ S ̄, y2 = 

Φ−1(1 − exp(−(γw)β))2 in T2 can be dominated

by a constant that does not depend on the parameter values, and that the resulting T2 integral is

finite.

Now consider T3.

We claim that for w large enough, y2 is “like” 2(γw)β. From (39) and (41) (with γ,β replaced

by γ0, β0), we know that for w>wup and θ ∈ S,

(γw)

β > (γlowwup)β > (γlowwup)βlow > (γ0w3/4)β0 (325)

so

y = Φ−1



1 − exp −(γw)

β

 > Φ−1(3/4) > 0

Thus, we can apply Lemma 1 to obtain

y2

1 + y2

1

√2π

exp −y2/2



< exp −(γw)

β



y <

1

√2π

exp −y2/2

 (326)

Taking logs, we then obtain

log  y2

1 + y2



+ log  1

√2π



− y2

2 < −(γw)

β + log(y) < log  1

√2π



− y2

2 (327)

or

log  y2

1 + y2



/(−y2/2) + log  1

√2π



/(−y2/2) + 1 > (γw)

β/(y2/2) + log(y)/(−y2/2)

> log  1

√2π



/(−y2/2) + 1 (328)

59
By (325), for w>wup and any θ ∈ S,

y = Φ−1



1 − exp −(γw)

β

 > Φ−1

1 − exp −(γlowwup)

βlow 

so we can make y arbitrarily large by choosing wup large enough. Thus, by result (328), given any

δ > 0, we can choose wup large enough so that w>wup implies

1 − δ < (γw)

β/(y2/2) < 1 + δ

or, taking δ = 1/2,

y2/4 < (γw)

β

or

y2 < 4(γw)

β < 4(γupw)βup (329)

Thus, for (γ, β, ρ, μ, σ) ∈ S, the y2 in T3 can be dominated by a function that does not depend on

the parameter values, and the resulting integral over (wup, ∞) is finite.

Now consider T1.

We claim that for w small enough, y2 is “like” −2β log(γw). From (38) and (40) (with γ,β

replaced by γ0, β0), we know that for w<wlow and θ ∈ S,

(γw)

β < (γupwlow)βlow < min 

1/4,(γ0w1/4)

β0



(330)

so

y = Φ−1(1 − exp 

−(γw)

β



< Φ−1(1/4) < 0

Thus, we can apply Lemma 1 to obtain

y2

1 + y2

1

√2π

exp −y2/2



<



1 − exp −(γw)

β

 (−y) <

1

√2π

exp −y2/2

 (331)

Taking logs, we then obtain

log  y2

1 + y2



+ log  1

√2π



− y2

2 < log 

1 − exp −(γw)

β

 + log(−y) < log  1

√2π



− y2

2 (332)

or

log  y2

1 + y2



/(−y2/2) + log  1

√2π



/(−y2/2) + 1 > − log 

1 − exp −(γw)

β

 /(y2/2)

+ log(−y)/(−y2/2)

> log  1

√2π



/(−y2/2) + 1 (333)

Because γupwlow < 1 (result (40)), for w<wlow and and θ ∈ S,

y = Φ−1



1 − exp −(γw)

β

 < Φ−1

1 − exp −(γupwlow)

βlow 

so we can make −y arbitrarily large by choosing wlow small enough. Thus, by result (333), given

any δ > 0, we can choose wlow small enough so that w<wlow implies

1 − δ < − log 1 − exp −(γw)

β

 /(y2/2) < 1 + δ (334)

60
Now by results (299) and (330), for w<wlow and θ ∈ S,

1 − exp −(γw)

β



> (γw)

β × 7/8

so

log 1 − exp −(γw)

β

 > β log(γw) + log(7/8)

or

− log 1 − exp −(γw)

β

 /(y2/2) < −β log(γw)/(y2/2) − log(7/8)/(y2/2) (335)

Combining results (334) and (335) and the fact that −y can be made arbitrarily large by making

wlow small enough, we know that we can choose wlow small enough so that for w<wlow and θ ∈ S

y2/4 < −β log(γw) (336)

By assumption (40), for w<wlow and θ ∈ S

1 > γupwlow > γw > γloww

so

β log(γw) > β log(γloww) > βup log(γloww) (337)

By results (336) and (337), for w<wlow (and wlow sufficiently small)

y2 < −4βup log(γloww) (338)

Thus, for (γ, β, ρ, μ, σ) ∈ S, the y2 in T1 can be dominated by a function that does not depend on

the parameter values, and the resulting integral over (0, wlow) is finite (because log(w)rwβ0−1 → 0

as w → 0).

Thus, we have established that if the parameter vector lies in S, then y2 (and hence y) is

dominated by an integrable function that does not depend on the parameter values.

Technique 5—handling wkβ and log(w)j

On (wup, ∞), wkβ log(w)j is dominated by wkβup+j which is integrable with respect to

γβ0

0 β0wβ0−1 exp 

−(γ0w)β0



.

For w in the compact interval [wlow, wup], and θ in the compact region S ̄, the continuous function

wkβ log(w)j is dominated by a constant which is integrable with respect to γβ0

0 β0wβ0−1 exp 

−(γ0w)β0



.

For w ∈ (0, wlow), |wkβ log(w)j | is dominated by | log(w)|

j which is integrable with respect to

γβ0

0 β0wβ0−1 exp 

−(γ0w)β0



because | log(w)|

jwβ0−1 → 0 as w → 0.

Thus, if the parameter vector lies in S, then |wkβ log(w)j | is dominated by an integrable function

that does not depend on the parameter values.

Applying the five techniques

We have already used techniques 1, 2, and 5 to establish that third partials (see results (102)

through (136)) 1, 6–9, 18, 22, and 27–31 satisfy Lehmannn’s condition D. Now we simply list the

remaining partials and the techniques needed to handle them. (Note that partials of y are given in

Appendix D.)

2: Techniques 1, 3, 4

3: Techniques 1, 3, 4

61
4: Techniques 1–5

5: Techniques 1–5

10: Techniques 1, 3, 4

11: Techniques 1, 3, 4

12: Techniques 1–4

13: Techniques 1–5

14: Techniques 1, 3, 4

15: Techniques 1, 3, 4

16: Techniques 1–4

17: Techniques 1–5

19: Techniques 1–5

20: Techniques 1–5

21: Techniques 1–5

23: Techniques 1–5

24: Techniques 1–5

25: Techniques 1–5

26: Techniques 1, 3, 4

32: Techniques 1–4

33: Techniques 1–5

34: Techniques 1–5

35: Techniques 1–5

20 Appendix K—Behavior of pseudo-truncated Weibull as ρ → 1

From result (6) it is clear that we need to investigate the behavior of

N/D ≡



Φ



(cu − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



− Φ



(cl − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



÷ (Φ((cu − μ)/σ) − Φ((cl − μ)/σ)) (339)

as ρ → 0 where

y = Φ−1



1 − exp −(γw)

β



We have

N = Φ 

cu − μ

σ

1 − ρ2 − ρy

1 − ρ2



− Φ



cl − μ

σ

1 − ρ2 − ρy

1 − ρ2



(340)

= Φ 

cu − μ

σ

1 − ρ2 − y

1 − ρ2

+ (1 − ρ)y

1 − ρ2



− Φ



cl − μ

σ

1 − ρ2 − y

1 − ρ2

+ (1 − ρ)y

1 − ρ2



≡ Φ (argu) − Φ (argl)

First note that

(1 − ρ)y

1 − ρ2 = (1 − ρ)y√1 − ρ

√1 + ρ

→ 0 (341)

62
as ρ → 1. Next note that by the definition (7) of wl, for w<wl

y = Φ−1



1 − exp −(γw)

β



< Φ−1



1 − exp −(γwl)

β

 = cl − μ

σ

(342)

Thus, from results (340)–(342), for w<wl both argl and argu converge to +∞ as ρ → 1, and,

consequently N → 1 − 1 = 0 as ρ → 1.

By the definition (8) of wu, for w>wu

y = Φ−1



1 − exp −(γw)

β



> Φ−1



1 − exp −(γwu)

β

 = cu − μ

σ

(343)

So from results (340), (341), and (343), for w>wu both argl and argu converge to −∞ as ρ → 1,

and, consequently N → 0 − 0 = 0 as ρ → 1.

For w ∈ (wl, wu)

y = Φ−1



1 − exp −(γw)

β



> Φ−1



1 − exp −(γwl)

β

 = cl − μ

σ

(344)

and

y = Φ−1



1 − exp −(γw)

β



< Φ−1



1 − exp −(γwu)

β

 = cu − μ

σ

(345)

So from results (340), (341), (344), and (345), for w ∈ (wl, wu), argl converges to −∞ and argu

converges to +∞ as ρ → 1, and, consequently N → 1 − 0 = 1 as ρ → 1.

For w = wl

y = Φ−1



1 − exp −(γwl)

β

 = cl − μ

σ

(346)

From results (340), (341), and (346), for w = wl, argl converges to 0 and argu converges to +∞ as

ρ → 1, and, consequently N → 1 − 1/2=1/2 as ρ → 1.

For w = wu

y = Φ−1



1 − exp −(γwu)

β

 = cu − μ

σ

(347)

From results (340), (341), and (347), for w = wu, argl converges to −∞ as ρ → 1 and argu converges

to 0, and, consequently N → 1/2 − 0=1/2 as ρ → 1.

21 Appendix L—Proof that pseudo-truncated Weibulls are not

Weibulls

We need to show that, given γ1 > 0, β1 > 1, there is no γ2 > 0, β2 > 0 such that, for all w > 0,

fPTW(w) = γβ1

1 β1wβ1−1 exp 

−(γ1w)

β1



×



Φ



(cu − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



− Φ



(cl − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



÷ (Φ((cu − μ)/σ) − Φ((cl − μ)/σ))

= γβ2

2 β2wβ2−1 exp 

−(γ2w)

β2



63
where

y = Φ−1



1 − exp −(γ1w)

β1



That is, we need to show that, for at least one w > 0,

R ≡



γβ1

1 /γβ22



(β1/β2) wβ1−β2 exp 

−(γ1w)

β1



exp (γ2w)

β2



(348)

×



Φ



(cu − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



− Φ



(cl − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



÷ (Φ((cu − μ)/σ) − Φ((cl − μ)/σ))

= 1

First, consider the case in which β1 ≥ β2. In this case, as w → 0, wβ1−β2 , exp 

−(γ1w)β1



, and

exp (γ2w)β2



stay bounded while

Φ



(cu − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



− Φ



(cl − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



converges to 0. (We are assuming that ρ = 0.) Thus, result (348) holds.

Now consider the case in which β1 < β2. Let w → ∞. We have

exp −(γ1w)

β1



exp (γ2w)

β2



= exp wβ2



γβ2

2 − γβ11 /wβ2−β1



and it is clear that for w large enough, this is greater than

exp kwβ2



(349)

for some k > 0.

Next, let au ≡ (cu − μ)/



σ

1 − ρ2



, al ≡ (cl − μ)/



σ

1 − ρ2



, and b ≡ ρ/1 − ρ2. Then,

for w large enough,

S ≡ Φ



(cu − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



− Φ



(cl − μ)/



σ

1 − ρ2



− ρy/1 − ρ2



=


 au−by

al−by

φ(s) ds

> (au − al) × φ(al − by)

> (au − al) × φ(2by)

= (au − al) ×



1/

√

2π

 

exp −y2/2

4b2

(350)

Now it is clear from results (348)–(350), that to determine the behavior of R for large w we

need only consider

T ≡ wβ1−β2 exp 

kwβ2



(φ(y))4b2

(351)

From Lemma 1, we know that for large w (and thus for large y) we can replace φ(y) by (1−Φ(y))y

to obtain

T = wβ1−β2 exp 

kwβ2



((1 − Φ(y))y)

4b2

M(w)

= wβ1−β2 exp 

kwβ2

 exp 

−(γ1w)

β1



y

4b2

M(w)

= wβ1−β2 exp 

wβ2



k − 4b2γβ1

1 wβ1−β2

 y4b2

M(w) (352)

64
where M(w) → 1 as w → ∞. For large enough w, this is greater than

w−β2 exp 

wβ2 k/2



y4b2M(w)

which (recall that β2 > 1) clearly converges to ∞ as w → ∞. Thus, result (348) holds.

65
Normal with mean 100, CV = .35

Weibull with median 100, CV = .35

 2e−05 

 4e−05 

 6e−05 

 8e−05 

 1e−04 

 0.00012 

 0.00014 

50 100 150

50 100 150

  	

 

 
     
 

   
Æ 
   
    


 
  

!!
Normal with mean 100, CV = .35

Weibull with median 100, CV = .35

 2e−05 

 4e−05 

 6e−05 

 8e−05 

 1e−04 

 0.00012 

 0.00014 

 0.00016 

50 100 150

50 100 150

  	

 

 
     
 

   
Æ 
   
    


 
 !

"!
Normal with mean 100, CV = .35

Weibull with median 100, CV = .35

 5e−05 

 1e−04 

 0.00015 

 2e−04 

 0.00025 

50 100 150

50 100 150

  	

 

 
     
 

   
Æ 
   
    


 
  

!"
Normal with mean 100, CV = .25

Weibull with median 100, CV = .25

 5e−05 

 1e−04 

 0.00015 

 2e−04 

 0.00025 

60 80 100 120 140

60 80 100 120 140

  	

 

 
     
 

   
Æ 
   
    


 
  

!"
Normal with mean 100, CV = .25

Weibull with median 100, CV = .25

 5e−05 

 1e−04 

 0.00015 

 2e−04 

 0.00025 

 3e−04 

60 80 100 120 140

60 80 100 120 140

  	

 

 
     
 

   
Æ 
   
    


 
  

 
Normal with mean 100, CV = .25

Weibull with median 100, CV = .25

 5e−05 

 5e−05 

 1e−04 

 0.00015 

 2e−04 

 0.00025 

 3e−04 

 0.00035 

 4e−04 

 0.00045 

 5e−04 

 0.00055 

60 80 100 120 140

60 80 100 120 140

  	

 

 
     
 

   
Æ 
   
    


 
 !

"#
Normal with mean 100, CV = .15

Weibull with median 100, CV = .15

 1e−04 

 2e−04 

 3e−04 

 4e−04 

 5e−04 

 6e−04 

 7e−04 

 8e−04 

60 80 100 120 140

60 80 100 120 140

  	

 

 
     
 

   
Æ 
   
    


 
  

!
Normal with mean 100, CV = .15

Weibull with median 100, CV = .15

 1e−04 

 2e−04 

 3e−04 

 4e−04 

 5e−04 

 6e−04 

 7e−04 

 8e−04 

 9e−04 

 0.001 

60 80 100 120 140

60 80 100 120 140

  	

 

 
     
 

   
Æ 
   
    


 
 !

!"
Normal with mean 100, CV = .15

Weibull with median 100, CV = .25

 2e−04 

 4e−04 

 6e−04 

 8e−04 

 0.001 

 0.0012 

 0.0014 

 0.0016 

60 80 100 120 140

60 80 100 120 140

  	

 

 
     
 

   
Æ 
   
    


 
 

!"
50 100 150

50 100 150

Predicted order statistics

Actual order statistics

 	 
 
 
   
 
  

 Æ         

          

 
60 80 100 120

60 80 100 120

Predicted order statistics

Actual order statistics

  	 


 
 
  
 	  

 
Æ 
 
  
    

 


       
  
 

!"