Symbolic Cholesky decomposition of the variance—covariance matrix of the negative multinomial distribution

Statistics & Probability North-Holland

Letters

28 September

15 (1992) 103-108

1992

Symbolic Cholesky decomposition of the variance-covariance matrix of the negative multinomial distribution Masahiko

Sagae

Department of Information Sciences, Faculty of Science and Technology, Science University of Tokyo, Japan

Kunio Tanabe The Institute of Statistical Mathematics, Tokyo, Japan Received December 1991 Revised February 1992

Abstract: This note shows a symbolic formula for the square-root-free Cholesky decomposition of the variance-covariance matrix of the negative multinomial distribution. A similar decomposition was given for the multinomial case by Tanabe and Sagae (1984). The evaluation of the symbolic Cholesky factors requires much less arithmetic operations than those with the general Cholesky algorithm. It is applied to obtain a recursive algorithm for generating multivariate normal random numbers which simulate samples from a negative multinomial population, which is similar to Pederson’s procedure for sampling from multinomial populations, An explicit formula of a multivariate normal density approximation to negative multinomial distribution is also given. Keywords: Symbolic number.

Cholesky

decomposition;

negative

multinomial

distribution;

multivariate

normal

approximation;

random

1. Introduction This note shows the square-root-free LDL’= of a positive

symbolic

Cholesky

decomposition,

[r] +rr’,

definite

matrix

(1) [rl + rrt, where

[rl =diag(~~/~~, PJP~,..-,PJP~), is an IZ x 12 diagonal

matrix,

r = ( P~/P~, is an II x 1 vector tpi< i=O

(2)

p2/p0,.

. . , Pn/Polt,

(3)

such that 1

and

p, >

0,

i=o,

1 ,...> n,

Correspondence to: Masahiko Sagae, Department of Information of Tokyo, 2641 Yamazaki, Nodashi, Chiba-ken 278, Japan. 0167-7152/92/$05.00

0 1992 - Elsevier

Science

Publishers

(4)

Sciences,

Faculty

B.V. All rights reserved

of Science

and Technology,

Science

University

103

Volume

15, Number

2

STATISTICS

& PROBABILITY

LETTERS

L = {Zi,j) is an IZX n lower triangular matrix with unit diagonal D = diag(d,, . . . , d,) is an n x n diagonal matrix.

elements,

28 September

1992

L’ is its transpose

and

Let the random number xi be the number of independent trials corresponding to mutually exclusive events Ei (i = 0, 1,. . . , n) and let pi be the probability of occuring Ei. Then, the joint distribution, Negative multinomial distribution, of x i, . . . , x, which are necessary to obtain m occurrences of an event E,, is defined by

(5) where pO = 1 - Cl_Ipi. Its variance-covariance matrix is given by m([r] + rrt) (e.g. Khatri and Mitra, 19691. The decomposition (1) is needed in sampling by a multivariate Normal approximation to this distribution. A Normal random vector x with mean 0 and the variance-covariance matrix m([r] + rr’) is generated by the transformation, x = &LD’/‘e, of a random vector E whose elements are independently identically normally distributed with mean 0 and variance 1. The numerical algorithm is generally used for the Cholesky decomposition of a positive definite matrix (e.g. Forsythe and Moler, 1967). If some of the pi’s are very small, then [r] + rrt is an ill-conditioned matrix. If this is the case, it is difficult to obtain an accurate Cholesky decomposition of the matrix with a finite precision calculation. However, a symbolic formula can be obtained for the Cholesky decomposition. The evaluation of the symbolic Cholesky factor matrices requires much less arithmetic operations than those with the general Cholesky algorithm. Besides, it is not affected by ill-condition of the matrix. The symbolic formula and related identities are given in Section 2. The decomposition is used to obtain a simple recursive algorithm for generating a Normal random vector which simulates a negative multinomial sample, which is similar to Pederson’s algorithm of sampling a multinomial distribution. A multivariate Normal density approximation to the negative Multinomial distribution is given in Section 3.

2. Symbolic Cholesky decomposition

of [r-l + t-r-*

Let a positive quantity qi be defined by

k=O

then we have pi = qi - qi_ 1, i = 1, 2,. . . , n, where pO = qo. The following theorem gives the exact square-root-free symbolic Cholesky decomposition [r] + rrt. Theorem

of the matrix

1. Given a vector p, there exists a decomposition, LDL’ = [r] + rrf,

(7)

where L = {li,j) is an n X n lower triangular matrk with i >.i,

i=j, i
(f-9

Volume

15, Number

STATISTICS

2

& PROBABILITY

LETTERS

28 September

1992

d,, . . . , d,) is an n x n diagonal matrix with

and D = diagcd,, di=

*,

i=1,2

n.

,...,

(9)

PO4i-1

The following Lemma

lemmas

are easily proved.

2. The inverse L ~ ’ of the matrix L is given by \

1 1

-32 L-1

=

-sx

-s3

-s,

-s,

1

(10)

7

-s,

...

-s,

I/

where L is defined by (8), and si’s are defined by Pi

s-= 1 4i-1

i= 1,2 ,...,

=lii_l, ’

n.

0

(11)

Lemma 3. The equality (i.e.,

r=Ls

L-‘r

= s)

(12)

holds where s is an n X 1 vector defined by s = (sl, s2,. . . , sJ. Proof of Theorem

1. By Lemma

L-‘([r]

where

L-’

3, we have

+ rr’)L-‘=

denotes

CL’)-‘.

q

L-l[r]L-’

+ L-‘rr’L_‘=

The (i, i)th element

‘p;l{pi

L-l[r]L-‘+ss’,

of this matrix

+sfqL_l}

is given by

= di.

= = PO4i-

Similarly,

if i > j, then the (i, j)th element

-I PO

sisj i

c I} 111

+

pk - s.p. L

S.S. I I

=p

‘PO

Since the matrix L-‘([r]

which implies

is symmetric, +rrt)L-’

the equality

of the matrix

(13)

(14)

1

is given by

tT’{‘i’j( -SiPi) PO+ l
lsi(

sjqj-

*

-pi,,.>

=

0.

(15)

we have = D,

(7).

(16)

0

The number of multiplicative operations required for evaluating the symbolic n2/2 + O(n) as compared with n3/6 + O(n*) with the general Cholesky algorithm.

Cholesky

factors

is

105

Volume

15, Number

2

STATISTICS

The following corollary gives the bounds be too large unless pa is too small. Corollary

& PROBABILITY

for the elements

LETTERS

of the matrices

28 September

1992

L, Lp ’ and D. They cannot

4. The following inequalities hold: i>j,

O
+si)


(17) +p;‘p;),

i= 1, 2,...,n,

i = 1, 2,. . .,n.


(18) (19)

Proof. Since, qj +pi =pO + Ci=,pk +pi < 1, i >j, we have (17). The inequalities seen by noting (17), and 0
(18) and (19) are easily

The diagonal element of the Cholesky factor mD of m([r] + rrt) is given the following interpretation: let pi be the probability at birth of a new born to die at age i, and m is the number of babies who dies at the age 0. The ith diagonal element of mD is given by

which is the variance of the number who die at age i, until the number who die under the age i attains its variance of Xi given X0 +X1 +X, expected value mq,_ 1/p,,, that is, its md, is the conditional + ... +x,-r =mq;_,/p, as follows, mdi=Var(X,]X,+X,+X,+ *.. +Xi_,=mq,_I/po) for i= 1>. . . 3n. The ith principal minor Mi of [r] + rrt, is given by M; =p;i-l( In particular,

p1p2 . . . pi)qi.

if the equality det([r]

+rr’)

(21)

holds in (4), then

=p;n-1(p,p2

. ..p.).

(22)

Pederson’s (1973, pp. 815-816) recursive method for generating random numbers of multinomial distribution is based essentially on the Cholesky decomposition of the corresponding variance-covariante matrix by Tanabe and Sagae (1984). A similar formula to Pederson’s can be obtained for the negative multinomial distribution. 5 (Recursive random number generator). An n-dimensional Normal random vector x = x,)~ with mean mr and variance-covariance matrix m([r] + rt) is generated by the recursive (Xl? x2,. . . , formulas Corollary

x1 = ms, + {md, e1

(23)

and xi=

(m+

icIxk)s1+i;;;;7ei,

i=2,3

,...,

n,

(24)

where E;‘s, i = 1, 2,. . . , n, are independent successive random choice of the standard Normal distribution, di and si are defined respectively in (9) and (11). 106

Volume

15, Number

STATISTICS

2

& PROBABILITY

28 September

LETTERS

1992

Proof. The desired random vector x is given by

x = mp + &LD’/2~, where D112 = diag(&, (25) by L-‘, we have

&,

L-lx = mL-‘p

(25) . . . , &>

+ &D1/2E

and E = (or, Ed,. . . , E,Y. Premultiplying

both sides of the equality

= ms + &D’12~,

by Lemma 3. Hence, x=L-‘x-

(L-1 -1)x=

ms + fiD”2~

- (L-’

- Z)x,

which implies (23) and (24), since L-’ is of the form (10). The corollary gives another interpretation of di and S, as follows; si is the expected value of the conditional mean of x,/Cm + CL:lIxk), and di is the expected value of the conditional variance of xi/ 6, given x1, x2,. . . , x,_~. 0

3. Normal approximation

of negative

multinomial

distribution

The inverse matrix of [r] + rr’ is a positive definite Stieltjes matrix given by ([r]

+rr’)_’ = [r]_’ -

PO

pll’,

YLOPk

where 1 is an n x 1 vector defined by 1 = (1, 1,. . . , 1)‘. In particular, if the equality holds in (4), then ([r] + I;r’)-’ = [r]-r -poll’. Now we give an explicit representation of Normal approximation of the negative multinomial distribution. Proposition 6. A multivariate Normal density approximation (Khatri and Mitra, 1969) of the negative multinomial distribution (5) is given by

_ (x-mr)‘([r]

+r’)-l(x-mr) 2m

i

= ( j&)1’2[p,r4”.1.

p, i1’2 ( xk - mp, 1pk)2

_

-mp;‘)’

(m + C;=lxk -1

mp; lpk

mpo

Proof. The result is easily seen from the fact that det([r] + rr’) =p1p2

11. . . . p,/pz+‘,

(27)

and the following

equality holds: (x-mr)‘([r]

+rrt)-‘(x-mr)/m

= (x-mr)t([r]-l = k k=l

(xk-mp,‘pk)’ mp0 lpk

-pollt)(x-mr)/m -

(m+E;=lxk-mp;1)2 mpo

-1

107

Volume

15, Number

2

STATISTICS

& PROBABILITY

LETTERS

28 September

1992

The numerical difficulty with the general Cholesky algorithm for computing the decompositions (7) is measured by the condition number given in the following corollary. The matrix [rl + mt is ill-conditioned when some of the elements pi, i = 0, 1, 2,. . . , n, are very small. q Corollary

7. We have

cond,([rl +rrt) = where cond,(A) = IIA II1 II A-’

IKiXjC~=~I ai,j I.

-$(y,y

pk

II 1 and

G=OPk +n .

_2 7

mlnk,lPk

IIA II1 is the matrix norm of A = Iai,j) defined

by

IIA

II 1=

0

By the corollary we have 1

1 -Po II

i

.

mlnk,lpk

+-

n-2

< cond,([r]

PO

+rrt)

Gi

1 . mmk,lpk

(29)

References Forsythe, G.E. and C.M. Moler (19671, Computer Solution of Linear Algebraic Systems (Prentice-Hall, Englewood Cliffs, NJ). Khatri, C.G. and S.K. Mitra (19691, Some identities and approximations concerning positive and negative multinomial distributions, in: P.R. Krishnaiah, ed., Multiuariute Analysis II, Proc. of the Second Internut. Symp. (Academic Press, New York) pp. 241-260. Pederson, D.G. (19731, An approximation method of sampling a multinomial population, Biometrics 29, 814-821. Tanabe, K. (19891, Symbolic LDL’ Decomposition of the

108

symmetric matrix D, + rrt, Res. Mem. No. 366, Inst. Statist. Math. (Tokyo). Tanabe, K. and M. Sagae (1984), The exact Cholesky decomposition and the generalized inverse of the variance-covariance matrix of the multinomial distribution and their applications, to appear in: J. Roy. Statist. Sot. Ser. B. Tanabe, K. and M. Sagae (19891, Symbolic Cholesky decomposition of the variance-covariance matrix of the negative multinomial distribution, Res. Mem. No. 365, Inst. Statist. Math. (Tokyo).