Implementación del càlculo de polinomios zonales y aplicaciones en análisis multivariante

Q¨ UESTII´ O,vol. 26, 3, p. 429-441, 2002 IMPLEMENTACI ´ ON DEL C ´ ALCULO DE POLINOMIOS ZONALES Y APLICACIONES EN AN´ ALISIS MULTIVARIANTE J. RODR´ IGUEZ-AVI A. J. S ´ AEZ-CASTILLO A. CONDE-S´ ANCHEZ Universidad de Ja´en∗ En este trabajo se describe la implementaci´ on de un algoritmo para el c´ alculo de polinomios zonales, as´ ı como dos aplicaciones expl´ ıcitas de ´ estos en el ´ ambito del an´ alisis multivariante. Concretamente, esta implementaci´ on permite obtener resultados de sumaci´ on aproximados para funciones hipergeom´ etricas de argumento matricial que, a su vez, pueden utilizarse en la g´ enesis de distribuciones multivariantes discretas con frecuencias sim´ etricas. De igual forma, se pone en pr´ actica un conocido resultado te´ orico que caracteriza la distribuci´ on de la menor ra´ ız caracter´ ıstica de una matriz aleatoria con distribuci´ on de Wishart. Implementing calculus of zonal polynomials and applications in Multivariate Analysis Palabras clave: Polinomioszonales, distribucionesdiscretas,distribuci´on deWishart,funcioneshipergeom´etricasdeargumentomatricial Clasificaci´ on AMS (MSC 2000): 60E05, 62E15 *Dep. de Estad´ıstica e I.O. Despacho 147, Ed. D-3, Paraje Las Lagunillas, Universidad de Ja´en, 23071 Ja´en, Espa˜na. Telf: + 34 953 012.207; fax: + 34 953 012.222; email: [email protected]. –Recibido en diciembre de 2000. –Aceptado en abril de 2002. 429 1. INTRODUCCI ´ ON Los polinomios zonales son una extensi´on multivariante de las funciones potenciales. La familiade polinomiossim´etricos y homog´eneosqueconstituyenha sidoampliamente utilizada dentro de la Estad´ıstica Matem´atica cl´asica en la expresi´on de densidades y distribuciones de formas cuadr´aticas en poblaciones normales, para extender resultados conocidos de la estad´ıstica univariante a un ambiente multivariante (Muirhead, 1982). No obstante, no se conocen f´ormulas generales para estos polinomios, ya que su definici´on, comovamos a ver a continuaci´on,los caracteriza como autofuncionesde un cierto operador diferencial. En concreto, dada una matriz sim´etrica Xm×mcon ra´ıces caracter´ısticas x1,...,xmy κ=(k1,...,km)una partici´on de ken no m´as de mpartes, el polinomio zonal de X correspondiente a κ, denotado por Cκ(X), es un polinomio sim´etrico homog´eneo de grado ken las ra´ıces caracter´ısticas x1,...,xmque debe verificar las siguientes tres condiciones: 1. El t´ermino de mayor peso es, salvo constante, xk1 1···xkm m. 2. Verifica la ecuaci´on diferencial (1) ∆XCκ(X)=ακCκ(X) donde ∆X= m ∑ i=1x2 i∂2 ∂x2 i + m ∑ i=1 m ∑ 1=j=i x2 i xi−xj ∂ ∂xi y ακ= m ∑ i=1 ki(ki−i)+k(m−1). 3. (trX)k=(x1+···+xm)k=∑κCκ(X), donde por ∑κdenotamos la suma en todas las particiones del grado k. Alternativamente a esta ´ultima condici´on, que se utiliza para la normalizaci´on de los polinomios zonales, puede considerarse la evaluaci´on de estos polinomios en un punto concreto. En este sentido (2) Cκ(Im)=22kk!1 2mκ p ∏ i<j(2ki−2kj−i+j) p ∏ i<j(2ki+p−i)! 430 donde pes el n´umero de partes no nulas de κy, en general, (3) (a)κ=(a)k1a−1 2k2 ·····a−1 2(m−1)km , con (a)k=a(a+1)···(a+k−1),(a)0=1. 2. ALGORITMO DE C ´ ALCULO Como se ha comentado, de la definici´on no puede deducirse una f´ormula expl´ıcita, aunque s´ı un algoritmo de c´alculo (James, 1968), que describimos seguidamente. Sea κ=(k1,...,km)una partici´on del entero k. Entonces, el monomio sim´etrico de una matriz Xm×mcon ra´ıces caracter´ısticas x1,...,xmcorrespondiente a κse define como Mκ(X)=xk1 1···xkm m+t´erminos sim´etricos. El algoritmo de James establece que los polinomioszonales son combinacioneslineales de monomios sim´etricos. Concretamente: (4) Cκ(X)= ∑ λ≤κ cκ,λMλ(X), donde la sumatoria se hace sobre todas las particiones λde ktales que λ≤κy (5) cκ,λ=∑ λ<µ≤κ [(li+t)−(lj−t)] ρκ−ρλ cκ,µ. En esta expresi´on, (6) ρκ= m ∑ i=1 ki(ki−i); adem´as, λ=(l1,...,lm)yµ=(l1,...,li+t,...,lj−t,...,lm),para t=1,...,lj,son tales que cuando las partes de la partici´on µse organizan en orden descendente, µes mayor que λy menor o igual que κen ordenlexicogr´afico. La sumatoria en (5) es sobre todas las posibles µ, y una suma vac´ıa se toma como cero. Mediante este algoritmo se determinan todos los coeficientes cκ,λsalvo el de mayor peso; es decir, se determina el polinomio zonal salvo normalizaci´on. Desde el punto de vista computacional, el algoritmo puede inicializarse considerando cκ,κ=1 y calculando los coeficientes cκ,λsalvo constante multiplicativa, que se obtiene mediante (2). 431 La implementaci´on del algoritmo se ha realizado empleando el programa MATLAB. Los correspondientes archivos fuente y las instrucciones para su ejecuci´on no aparecen aqu´ı por brevedad, aunque pueden encontrarse en la direcci´on http://www.ujaen.es/dep/estinv/pdi/ajsaez/software.html. 3. APLICACIONES En este apartadose presentansendas aplicaciones enel ´ambito dela estad´ıstica multivariante que pueden llevarse a cabo gracias al conocimiento expl´ıcito de un buen n´umero de polinomios zonales calculados mediante la metodolog´ıa expuesta en la secci´on anterior. 3.1. Aproximaci´ onde funciones hipergeom´ etricas de argumento matricial y g´ enesis de distribuciones multivariantes Se define la funci´on hipergeom´etrica de argumentomatricial de par´ametros a1,...,apy b1,...,bqasociada a la matriz sim´etrica Xm×mcomo (7) pFq(a1,...,ap;b1,...,bq;X)= ∞ ∑ k=0∑ κ (a1)κ···(ap)κ (b1)κ···(bq)κ Cκ(X) k!, donde ∑κdenota la suma sobre todas las particiones κ=(k1,...,km)de kyCκ(X)es el polinomio zonal de Xcorrespondiente a κ. Ning´un par´ametro bjdel denominador puede ser cero o un medio de un entero menor o igual que 1 2(m−1),ya que, en caso contrario, alguno de los denominadoresde la serie se anular´ıa. De igual forma, si alg´un par´ametro del numerador es un entero negativo, esto es, ai=−n, entonces la funci´on es un polinomio de grado mn. La serie convergepara todo Xsi p≤q; converge para X<1sip=q+1, donde por Xnotamos el m´aximo de los valores absolutos de las ra´ıces caracter´ısticas de X;y salvo que la serie tenga un n´umero finito de t´erminos no nulos, divergepara todo X=0 si p>q+1. En el caso m=1,la serie (7) se reduce la funci´on hipergeom´etrica generalizada cl´asica pFq(a1,...,ap;b1,...,bq;x)= ∞ ∑ k=0 (a1)k·····(ap)k (b1)k·····(bq)k xk k!. Podemos subrayar dos casos especiales de (7): los dados por la funci´on 0F0(X)= etr(X), generalizaci´onde la serie exponencial,y la funci´on 1F0(a;X)=det(Im−X)−a, 432 generalizaci´onde la serie binomial(1−x)−a.M´as all´a de estos dos casos no existenexpresiones generales de las funciones hipergeom´etricas de argumento matricial. Es por ello que para obtener resultados de sumaci´on se suele acudir a aproximaciones num´ericas mediante t´ecnicas de an´alisis que proporcionanel valor de las funciones en algunos casos concretos. Nosotros proponemos desarrollarestas funciones como series de monomios sim´etricos, utilizando la implementaci´on del algoritmo de James. De esta forma, pueden obtenerse resultados aproximadosde sumaci´on de dichas funcionespara un ampliorangode valores de los par´ametros. Asimismo, este desarrollo en serie es el adecuado para utilizarlas como generatrices de distribuciones multivariantes discretas de probabilidad, distribuciones que cuentan con una interesante propiedad de simetr´ıa en sus frecuencias. Tabla 1. Desarrollo truncado de la funci´on 2F1(a,b;c;I2). a b c 2F1(a,b;c;I2)∑20 k=0∑ρfρMρ(I2)Diferencia 2 2 10 3.0569 3.0547 0.0022 3 3 10 20.0909 19.0750 1.0159 5 5 25 13.8359 13.8299 0.006 4 5 20 15.7342 15.7046 0.0296 8 10 50 55.1096 55.0856 0.024 10 10 50 172.7901 172.2222 0.5679 12 12 65 260.1934 259.5012 0.6922 15 15 85 699.2147 695.7047 3.51 Concretamente, si en el desarrollo de la serie hipergeom´etrica dada en (7), expresamos los polinomios zonales en t´erminos de los monomios sim´etricos, se tiene que pFq(a1,...,ap;b1,...,bq;X)= ∞ ∑ k=0∑ κ (a1)κ···(ap)κ (b1)κ···(bq)κ Cκ(X) k! (8) = ∞ ∑ k=0∑ κ (a1)κ···(ap)κ (b1)κ···(bq)κ ∑λ≤κcκ,λMλ(X) k! = ∞ ∑ k=0∑ κ∑ λ≤κ 1 k! (a1)κ···(ap)κ (b1)κ···(bq)κ cκ,λMλ(X) = ∞ ∑ k=0∑ ρ fρMρ(X), donde de nuevo la suma ∑ρes en todas las particiones ρde cada entero k. Una vez calculadoslos coeficientes fρ, la serie puede aproximarsemediantesus sumas parciales. 433 En orden a considerar la bondad de estas aproximaciones, hay que tener en cuenta que la cola de la serie truncada en el grado r, contiene tan s´olo t´erminos en las ra´ıces caracter´ısticas de Xde grado superior a r. En la Tabla 1 se muestran resultados aproximados de sumaci´on referidos a la funci´on bivariante 2F1,para distintos valores de los par´ametros cuando X=I2. Todos ellos se comparan con el valor exacto de la suma, conocido gracias al Teorema de Sumaci´on de Gauss en su versi´on multivariante. En general, podemosestablecer una cota del error, ε, cometido mediante la suma de los M+1 primeros t´erminos cuando a1,a2,b1>0,b1> a1+a2yX≤1, dada por (9) ∞ ∑ k=M+1∑ κ (a1)κ(a2)κ (b1)κ Cκ(X) k!≤ ∞ ∑ k=M+1∑ κ Cκ(X) k! = ∞ ∑ k=M+1 traza(X)k k! ≤ ∞ ∑ k=M+1 mk k! =em− M ∑ k=0 mk k! Para una discusi´on m´as detallada de la metodolog´ıa de este m´etodo de sumaci´on puede verse Guti´errez et al. (2000). Por otra parte, las funciones hipergeom´etricas cl´asicas de argumento escalar, fundamentalmente la funci´on de Gauss 2F1y sus extensiones y particularizaciones, pueden emplearse paragenerar las distribuciones discretas m´as usuales (Johnson,Kotz yKemp, 1992;Guti´errez y Rodr´ıguez,1997;Rodr´ıguez et al.,2001a). Tambi´en se han introducidoextensionesbivariantesdealgunasfuncioneshipergeom´etricasdeargumentoescalar, como la F1olaF3,para generar distribuciones bivariantes discretas. En este sentido, las funciones hipergeom´etricasde argumentomatricial tambi´en pueden utilizarse para generar distribuciones discretas, en este caso multivariantes, que generalicen sus hom´ologas univariantes. Para ello tan s´olo es necesario considerar el desarrollo en serie dado en (8), de manera que las probabilidades de estas distribuciones vienen dadas por las constantes fρconvenientementenormalizadas. La principal caracter´ıstica de estas distribuciones es que sus frecuencias son invariantes frente a permutaciones de las variables que la forman; es decir, si Xm×1es un vector aleatorio generado por una funci´on de este tipo, se verifica: (10) P[X=(i1,...,im)] = P[X=σ(i1,...,im)], donde (i1,...,im)es cualquier vector posible para Xyσes cualquier permutaci´on. Esta propiedad hace que estas distribuciones sean modelos adecuados para fen´omenos 434 multivariantes discretos cuyas componentes no son necesariamente independientes pero s´ıid´enticamente distribuidas, de manera que las probabilidades de ocurrencia no dependan del orden de las variables. As´ı, se han utilizado con ´exito, por ejemplo, en la modelizaci´on del n´umero de ingresos diarios en dos boxes del Servicio de Urgencias del Hospital de San Agust´ın (Linares, Espa˜na) (Rodr´ıguez et al, 2001b). Otros ´ambitos adecuados de aplicaci´on aparecen, por ejemplo, en la modelizaci´on del n´umero de fallos en varias m´aquinas indistinguibles en paralelo o el n´umero de llamadas a varios servidores id´enticos dispuestos tambi´en en paralelo. En la Figura 1 se muestran ejemplos de poligonales de frecuencias generadas mediante funciones 1F1y2F1bivariantes. 0 5 10 15 20 0 5 10 15 20 0 0.005 0.01 0.015 0.02 0.025 0 5 10 15 20 0 5 10 15 20 0 0.005 0.01 0.015 0.02 0.025 0 5 10 15 20 0 5 10 15 20 0 1 2 3 4 5 6 7 8 x 10−3 0 5 10 15 20 0 5 10 15 20 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 Figura 1. Poligonal de frecuencias para distribuciones generadas por funciones 1F1(a.yb.)y 2F1(c.yd.). 3.2. Distribuci´ on de la menor ra´ız caracter´ıstica de una matriz de Wishart en t´ erminos de polinomios zonales Dentro de la Estad´ıstica Matem´atica, los polinomios zonales y las funciones hipergeom´etricas de argumento matricial han sido utilizados, entre otras aplicaciones, en la expresi´onde la distribuci´on de formas cuadr´aticas asociadas a poblaciones normales. As´ı, por ejemplo, expresiones que involucran funciones hipergeom´etricas aparecen en 435 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 00.5 11.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 2. Funciones de distribuci´on te´orica y emp´ırica para lamenor ra´ız caracter´ıstica de A1 W3(10,Σ1)(a.)yA2W3(10,Σ2)(b.). la distribuci´on χ2(Pearson, 1900), de Wishart central y no central (Wishart, 1928), de las ra´ıces caracter´ısticas de la matriz de covarianza de una distribuci´on gaussiana (James, 1964), en la ratio de la varianza F(Fisher, 1924), de la Fno central con pyn grados de libertad (Fisher, 1924), etc. En este apartado describimos, precisamente, una aplicaci´on de la implementaci´on realizada que nos permite evaluar la distribuci´on de probabilidad de la menor ra´ız caracter´ıstica de una matriz de Wishart. Concretamente (Muirhead, 1982), si lmes la menor ra´ız caracter´ıstica de S, donde A= nS esWm(n,Σ),yr=1 2(n−m−1)es un entero positivo, entonces (11) PΣ(lm>x)=etr−1 2nxΣ−1mr ∑ k=0∑ κ ∗Cκ1 2nxΣ−1 k!, donde ∑κ∗denota la suma sobre todas las particiones κ=(k1,...,km)de ktales que k1≤r. La puesta en pr´actica de este teorema tan s´olo implica el desarrollo de una serie finita en t´erminos de polinomios zonales, as´ı que puede realizarse de manera exacta hasta distintos valores de myn,aunque no muy elevados. Nosotros hemos considerado una aplicaci´on en donde m=3yn=10. Se han simulado mediante el programa matem´atico MATLAB 1.000 valores de matrices aleatorias Xn×mN(0,In⊗Σ). A partir de estas 1.000 matrices, x1 n×m,...,x1,000 n×m,se calculan 1.000 valores simulados de sendas distribuciones de Wishart, mediante la expresi´on Ai=xi·xi,i=1,...,1,000,de modo que la distribuci´on te´orica de estas matrices es la de AW3(10,Σ). Posteriormente se han dividido estas matrices por n=10 y a cada una de estas nuevas matrices se les ha calculado su menor ra´ız caracter´ıstica. Esto genera una muestra aleatoria simple, l1,...,l1,000,de 1.000 valores de la variable aleatoria lm,menor ra´ız caracter´ıstica de A/n. En la Figura 2 aparecen superpuestas la funci´on de distribuci´on emp´ırica, calculada a partir de sendas muestras aleatorias simples, y la funci´on de distribuci´on te´orica, dada 436 por (11), para matrices Σ1=⎛ ⎝ 10,58 0,61 0,5810,58 0,61 0,58 1 ⎞ ⎠,Σ2=⎛ ⎝ 6,83 6,40 2,84 6,40 8,02 2,79 2,84 2,79 3,98 ⎞ ⎠ Como puede comprobarse, son pr´acticamente id´enticas, confirmando que los datos emp´ıricos quedan perfectamente ajustados mediante la distribuci´on te´orica. 4. DISCUSI ´ ON La metodolog´ıa descrita en este trabajo, que involucra el c´alculo expl´ıcito de los polinomios zonales, proporciona la evaluaci´on, exacta en algunos casos y aproximada en otros, de cualquier serie hipergeom´etrica de argumentomatricial. Esto permite la puesta en pr´actica de numerosos resultados que involucran a estas series y que hasta ahora, o no se hab´ıan aplicado a casos reales, o si se hab´ıan aplicado, ten´ıan que limitarse a casos donde t´ecnicas de an´alisis num´erico permitieran alg´un resultado concreto. A su vez, tambi´en permite abrir nuevas l´ıneas de investigaci´on, como la de la g´enesis de distribuciones multivariantes de probabilidad mediante funciones hipergeom´etricas de argumento matricial. Hay que decir, no obstante, que el ´exito de esta metodolog´ıa est´a supeditado a las limitaciones computacionales en el c´alculo de los polinomios zonales y, sobre todo, de los monomios sim´etricos en t´erminos de los cuales vienen definidos. 5. AGRADECIMIENTOS Deseamos agradecer a los referees sus comentarios y sugerencias, los cuales han sido de gran ayuda a la hora de elaborar este manuscrito. REFERENCIAS Guti´errez, R. and Rodr´ıguez, J. (1997). «Family of Pearson Discrete Distributions Generated by the Univariate Hypergeometricfunction 3F2(α1,α2,α3;γ1,γ2;λ)».Appl. Stochastics Models and Data Anal, 13, 115-125. Guti´errez R., Rodr´ıguez J. and S´aez A. J. (2000). «Approximation of hypergeometric functions of matricial argument through their development in series of zonal polynomials».Electronic Transactions in Numeric Analysis, 11, 121-130. 437