1 / 69

Técnicas de Análisis aplicadas al estudio de la variabilidad del sistema climático

Técnicas de Análisis aplicadas al estudio de la variabilidad del sistema climático. Rosa Hilda Compagnucci 2008. Análisis de Componentes Principales (ACP).

kovit
Download Presentation

Técnicas de Análisis aplicadas al estudio de la variabilidad del sistema climático

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. Técnicas de Análisis aplicadas al estudio de la variabilidad del sistema climático Rosa Hilda Compagnucci 2008

  2. Análisis de Componentes Principales (ACP) • La forma más popular del Análisis de Factores es el ACP y consiste tan sólo el análisis de autovalores y autovectores de la matriz de varianza ó covarianza de una matriz de datos. El esquema moderno del análisis de factores emplea ACP como punto de partida el ACP • Provee una significativa visión de la estructura de la matriz y puede ser frecuentemente interpretado de la misma manera que los factores. • Analiza la interdependencia entre los elementos de una matriz intacta (no particionada) • Procedimiento matemático que transforma un conjunto de variables correlacionadas en un conjunto menor de variables no correlacionadas llamadas componentes principales. Al observar cuidadosamente este nuevo conjunto de variables no correlacionadas, se pueden obtener respuestas a muchas preguntas importantes. • Se obtienen las composiciones lineares de las variables originales llamadas “scores” (componentes principales, rasgos latentes, etc) que contienen ciertas propiedades deseables como máxima varianza y sin correlación con las previamente calculadas. • El procedimiento se focaliza en las variables de la matriz. • El propósito del ACP es reducir el espacio de operación simplificando los datos originales mediante la representación de los mismos objetos en un número menor de variables (dimensiones). • Por otra parte el análisis de “clusters” es el problema inverso, simplifica la matriz de datos reduciendo el número de objetos.

  3. ACP • Pocas • Conceptualmente significantes • Relativamente independientes • Diversas • Difíciles de interpretar • Correlacionadas COMPONENTES VARIABLES

  4. Objetivos del ACP: • “Desenredar” patrones complejos de asociaciones intervariables en datos multivariados. • Investigación exploratoria en la identificación de características “latentes” para futuros experimentos. • Desarrollar tipologías empíricas de las variables. • Reducir la dimensión del grupo de datos multivariados. • Desarrollar un índice unidimensional que máximamente separe a los individuos de la base de datos. • Testear la hipótesis de relación entre ciertas variables. • Transformar la matriz de variables predictoras previamente a aplicarle técnicas como regresión múltiple, o correlación canónica. • Escalamiento y representación espacial de los datos percibidos y preferenciales. • Determinar los “outliers”.

  5. Pasos del ACP: • Rotación de la configuración inicial de puntos (objetos) en nuevas orientaciones, de la misma dimensionalidad, que se caracterizan por ser mutuamente ortogonales y secuencialmente con máxima varianza. Esto es que la proyección de los puntos en la primera dimensión tiene la máxima varianza, la segunda la siguiente mayor varianza. • Reducir la dimensionalidad del espacio transformado descartando las altas dimensiones que exhiben la menor varianza en la proyección de los puntos. • Luego encontrar nuevas orientaciones del espacio reducido (rotación) que hagan que las dimensiones retenidas sean mejor interpretables desde un determinado punto de vista • Interpretación substancial de las dimensiones re-orientadas en términos de las variables que están altamente asociadas con cada una de dichas dimensiones.

  6. Procedimiento del ACP • Consideremos la matriz [X] mxn de datos que contiene m observaciones de dimensión n (número de variables potencialmente inter-correlacionadas) • Las matrices de inter-relación entre las n variables son: [B] Suma en “bruto” de los cuadrados y producto cruzado [S] Correlación-promedio de la suma de los cuadrados y producto cruzado [C] Varianzas y covarianzas [R] Correlaciones • Las matrices difieren en diferencias en medias y/ó desviación estándar entre las variables de interés. • Frecuentemente las investigaciones involucran variables cuyas medias y escala de unidades son arbitrarias con lo cual se utilizan la matriz [Xs] de datos estandarizada (por variables = columnas) ó lo que es igual la matriz de correlación [R] de nxn.

  7. Varianzas y covarianzas como coordenadas [C] : Matriz de varianza / co-varianza Dos variables (n=2) medidas sobre un número de objetos (m=25) Varianza de X1 = 20,3 y de X2 = 24,1 Covarianza (X1,X2) = 15,6 1º Ejemplo numérico Diagrama de dispersión del dato bi-variado en el sistema original La estructura de la matriz varianza-covarianza nxn de los datos puede mostrarse como dos vectores

  8. Los elementos de la matriz nxn pueden considerarse puntos ubicados en una elipse n-dimensional • Los autovectores de la matriz son los ejes principales del elipsoide • Los autovalores son el largo de dichos ejes • Como la matriz de varianza-covarianza es simétrica los autovalores son reales. También son soluciones los vectores multiplicados por una constante  en especial = -1 (aplicado en este caso) La varianza total del grupo de datos es igual a la suma de las varianzas individuales. Equivalente a la traza de [C]. VarTot.= 20,28 + 24,06 = 44,34 X1: 20,3/44,34 = 46%X2: 24,1/44,34=54% Los nuevos ejes (componentes) explican el valor del autovalor dividido por la traza de [C]: 1ºu : 37,8/44,34 = 86% y 2ºu : 6,47/44,34 = 14% Elipse definida por los vectores de varianza y covarianza de la matriz [C].

  9. Componentes Principales, Multiplicadores ó Factor “Score” y “Loadings” La transformación Z1 = 1 X1 + 2 X2 donde (1 ; 2) son los elementos del primer autovector (1ºu) Z2 = 1 X1 + 2 X2 donde (1 ; 2) son los elementos del segundo autovector (2ºu) Crea un nuevo grupo de datos (dos nuevas variables), medidos a lo largo de los ejes ortogonales entre sí y que tienen respectivamente varianza 37,9 y 6,5 Las variables originales del ejemplo (X1 ; X2) son convertidas en nuevas variables (Z1 ; Z2) llamadas “scores”, proyectando las observacionesoriginales sobre los ejes principales mediante: Z1i = 0,66 X1i + 0,75 X2i proyección de las observaciones sobre el 1º eje Z2i = 0,75 X1i – 0,66 X2i proyección de las observaciones sobre el 2º eje Los elementos de los autovalores utilizados para calcular los “scores” son llamados “loadings”[Z]=[X]*[U] Los “loadings” son simplemente los coeficientes de la ecuación linear definida por el autovector. En términos del Análisis de Factor sería “ el loading de la variable X en el 1º Factor y el loading en el 2º Factor”

  10. Tabla con los CPscore del ejemplo1 calculado al proyectar los datos originales en los ejes principales (autovectores) Diagrama de dispersión de los nuevos datos bi-variados: proyecciones de los datos en las CPscores

  11. La matriz de varianza/co-varianza de los CPscores es: • Los CPscores son ortogonales • La suma de varianza de los datos proyectados sobre las CPs es igual a la varianza a lo largo de X1 y X2 • Varianzas de los CPscores: • Z1 = 37,9  37,9/44,34  85,36% • Z2 = 6,5  6,5/44,34  14,64% • Si se descartaban inicialmente X1 ó X2 se perdían 46% ó 54% de la varianza. • Ahora si se descarta Z2 sólo se pierde el 14,64% de la varianza. Al eliminar el segundo eje (autovector) resulta una inevitable pérdida de varianza.

  12. Al ordenar por rango las variables X1 y X2 del ejemplo1 las nuevas X1 y X2 tienen menor dispersión en el diagrama, la covarianza de 15,59 pasa a 21,92 aumentando, sin embargo sus varianzas permanecen inalterables. La nueva matriz [C] es: La nueva matriz [C] sin embargo tiene casi los mismos autovectores que la anterior Sin embargo los dos nuevos autovalores son radicalmente diferentes. El 2º autovalor es casi nulo y la elipse se angosta concentrándose sobre el eje mayor que aumenta su longitud (casi una línea) La 1ºCP (eje principal) ahora explica 99,61% mientras que la 2ºCP sólo explica 0,39% de la varianza total de los datos. Obviamente podemos descartar la 2ºCP reduciendo la dimensionalidad sin perder varianza.

  13. Si se randomizan las variables X1 y X2 originales las varianzas e las variables se mantienen, pero disminuye la co-varianza Los dos autovalores son casi idénticos  Elipse casi circular La correlación entre las dos variables es casi cero, luego los dos ejes originales son casi ortogonales y casi iguales en magnitud, los ejes definen casi una elipse. Ningún otro par de ejes, ni siquiera los obtenidos por el ACP serán significativamente mejores que las dos variables originales. En esta situación no existe transformación de los datos que pueda reducir las dimensiones sin perder información.

  14. Datos estandarizados  Matriz [R] de correlación Conjunto de m datos bi-dimensionales : dos variables estandarizadas El centro de gravedad de los puntos será el origen de coordenadas - el valor medio de ambas variables es cero- y la escala en ambos ejes será la misma - las variables son adimensionales, con varianza unitaria Diagrama bidimensional en el que se representa una situación de dos variables estandarizadas con m puntos de observacion. El eje Z1 maximiza la varianza de los puntos proyectados sobre él. contiene las proyecciones de las observaciones sobre el eje Z1 (nueva variable que extracta la mayor cantidad de información contenida en xs1 y xs2 )

  15. Consecuencia del ACP es la rotación del espacio de datos original en nuevas orientaciones que son las componentes principales Si las variables originales no tuvieran media nula, el centro de gravedad de los puntos observacionales no sería el origen de coordenadas. En este caso, el eje principal Z1 (y/o eventualmente el Z2) no necesariamente debería pasar por el origen. En adelante supondremos que las variables están estandarizadas (o al menos tienen media nula), con lo cual nuestra matriz de datos será siempre la matriz estandarizada por columnas Xs (o al menos la matriz de anomalías Xd). Al considerar sólo dos variables estandarizadas es fácil ver que, de existir, los ejes Z1 y Z2 formarán con los ejes coordenados X1 y X2 un ángulo de 45º. Esto se debe a que las varianzas de los vectores observación proyectados sobre X1 y X2 son iguales. Si las variables no estuvieran estandarizadas estas varianzas no serían iguales y por lo tanto, la mayor o menor varianza sobre uno de éstos ejes haría que los ejes Z1 y Z2 no necesariamente deban formar ángulos de 45º con X1 y X2.

  16. (b) X X (a) 2 2 Z 1 Z 1 Z 2 X X 1 1 (d) X (c) X 2 2 Z 1 Z 2 X X 1 1 Z 3 X 3 • dos variables perfectamente correlacionadas • (b) dos variables no perfectamente correlacionadas. • (c) tres variables no perfectamente correlacionadas • (d) dos variables no correlacionadas. (Reproducido de Green 1978)

  17. La varianza explicada por cada nueva variable z1, z2 y z3 decrece sucesivamente y su suma resulta igual a la varianza total del conjunto de variables originales xs1, xs2 y xs3. • Como las variables están estandarizadas la varianza total es igual al número de variables = n. • Cuando los datos están en tres dimensiones la varianza = 3. Si el rango de la matriz de datos no fuese máximo, significaría que existe alguna dependencia lineal perfecta entre las variables originales y en consecuencia, el número de variables z necesario para explicar la variabilidad contenida en el conjunto de datos se reduce. • Las nuevas variables: componentes principales con (1jn) están formadas por las proyecciones de los puntos de observacion oi (1im) sobre los ejes Zj . Las proyecciones correspondientes z1j,z2j,...,zmj se llaman factores de puntaje. • El Análisis de Componentes Principales (ACP) consiste en calcular los ejes principales Zj para definir las nuevas variables Zj, reteniendo sólo aquellas que expliquen la mayor porción de varianza contenida en el sistema y desechando las restantes.

  18. Cálculo de las Componentes Principales basadas en variables estandarizadas: matriz de correlación [R] zi1proyección de la observación oi = (xi1, xi2) sobre el vector u de norma 1 que es la dirección del eje principal Z1 ángulo  entre oi y u vector u que maximice f : el promedio cuadrático sobre todo i del producto  Teorema del multiplicador de Lagrange: máximaf(u)implica que se verifique que usea autovector unitario de [R] Esquema bidimensional que muestra la proyección de una observación sobre el primer eje principal Z1. El vector ues unitario y colineal con el eje Z1. Se calculan los autovalores por el Método de Jacobi Ru = u

  19. El ACP de variables estandarizadas involucra el problema de diagonalización de la matriz de correlaciones [R] de [X] mxn con m observacionesy n variables. R = U D U´ (1) Donde U es una matriz cuyas columnas son los autovectores de [R] y cumple que U´U = UU´= I Las columnas de U pueden verse como un grupo de direcciones coseno que rotan la matriz [Xs] D es una matriz diagonal cuyos elementos son los autovalores de [R] Los elemontos no nulos de la diagonal de D pueden interpretarse como las varianzas de los puntos proyectados sobre los nuevos ejes, componentes principales (CPs), de [Xs] que son las columnas de [Z] Las componentes principales (CP) son las columnas de [Z] Z = XsU (2)

  20. Distintas matrices de similitud 2º Ejemplo numérico matriz de datos : [ X ] matriz de desvios:[ Xd ] matriz estandarizada: [Xs ] [S]: XdT*Xd producto cruzado de matriz de desvíos [B]: XT*X producto cruzado de matriz de datos [C]:matriz de covarianza [R]:matriz de correlación

  21. xs3 xs2 z1 z1 xs1 xs1 xs3 z1 xs2 ACP de la matriz de correlación[R], matriz de producto cruzado asociada a la matriz dedatos estandarizados Xs El diagrama de dispersión de los pares de variables Xs1, Xs2 y Xs3 de los datos del Ejemplo2 estandarizados exhiben alta correlación positiva. ¿Puede existir sólo un eje -composición linear - cuyos “scores” derivados (proyecciones de los puntos sobre él) máximamente separe s los 12 individuos?. Z1 (1º “score”)primer eje principal que maximiza la varianza de los puntos proyectados sobre él. Se puede obtener un segundo eje Z2 que extrae la máxima varianza remanente y que resulta ortogonal a Z1 y finalmente un tercero Z3 ortogonal a los anteriores que explica el total del remanente.

  22. DESCOMPOSICION EN VALORES SINGULARES (SVD) de [Xs] SDV de una matriz [A] es expresarla como producto de otras tres matrices A = P  Q’ mxn mxr rxr rxn Py Q ortonormalesy  (rxr)es una matriz diagonal de valores ordenados positivos P(mxm)  autovectores de AA’ Q (nxn)  autovectores de A’A  (rxr)  los elementos de la diagonal son las raíces cuadradas de 2, la matriz de los autovalores de AA’ ó de A’A  2 = D Por lo cual AA’=PDP’ A’A = QDQ’ El rango de A está dado por r : número de valores positivos de  (valores singulares) r (A)  min (m , n) Se cumple que  P’P=I y si A tiene rango igual al número de columnas r=n  Q es cuadrada yortogonal Se cumple que  Q’Q=QQ’=I Para la matriz estandarizada de datos [Xs]  Xs = P  Q’ (3)

  23. Xs P  Q’  Varianza de [Z] = ACP / Descomposición en Valores Singulares (SVD) de [Xs] Xs = P  Q’ Multiplicando a ambos lados por Q : la matriz de rotación cuyas columnas dan las direcciones coseno (autovectores)  Xs Q = PQ’Q = P = Z  se producen las tres componentes principales, nuevos ejes, componentes “scores”, columnas de [Z] no-correlacionadas, todas mutuamente ortogonales. Las varianzas de las columnas de [Z] declinan sucesivamente Matriz de covarianza de [Z]: C(Z)

  24. Una forma directa de encontrar las componentes principales Z de la matriz Xs -con variables en las columnas - es encontrar la descomposición en valores singulares de Xs Si sólo se retiene la primer componente z1 se aproxima Xs como: Donde 1 = 5,434 es el primer elemento de la diagonal de  y p1 es la primer columna de P y q1 la primer fila de Q’ Aunque sea de 12x3 tiene rango 1 y es la aproximación más cercana de Xs con rango 1 que pueda ser encontrada Con las componentes z2 y z3 se puede aproximar Xs por y tal que Ejemplo2 las covarianzas de [Z] son cero (las columnas de Z son ortogonales) las varianzas de las tres columnas declinan sucesivamente la suma de las varianzas es: 2,685 + 0,250 + 0,065 = 3 Que es también la traza de la matriz de correlación [R]

  25. RUDU’ Componentes Principales obtenidas de la matriz de correlación [R] Promedio del cross producto de la matriz de datos estandarizada - SIMETRICA Matriz de datos X (mxn)  Descomposición en Valores Singulares de R U matriz de autovectores de R D matriz diagonal con autovalores de R (4)

  26. La matriz de autovalores de R es igual a matriz de autovalores de Xs al cuadrado dividido por (m-1) • La matriz U de autovectores de R es igual a Q los autovectores de Xs’Xs • La matriz de rotación U puede ser encontrada por SVD tanto de R como de Xs • SVD de Xs directamente produce las CPs  (1) Z=P y la varianza de las CPs dada por (2) 2/(m-1) =D Donde i son los autovalores y los ui son los autovectores (columnas de U) con i=1….n R puede ser representada por de rango r ( n) dependiendo r del número de dimensiones retenidas r < n  Uno es cuadrada sino de nxr  U’U=I pero U U’  I (Ur :ortonormal sección de U) útil cuando es posible reducir las dimensiones con poca pérdida de información. La ecuación (4) puede escribirse también como: (5)

  27. Varianza de [Z] = Componentes “scores” y Componentes “loadings” Z=P=XsU ZU’=Xs Al escalar las componentes principales por su varianza el nuevo grupo de componentes principales Zstienen varianza unitaria (6) Estandarizar las componentes principales a varianza unitaria tiene el efecto de: 1º rotar Xs vía U a la orientación de las componentes principales (direcciones de los ejes del hiperelipsoide de puntos) 2º acortar los ejes largos y viceversa llevándolos a todos al largo unitario, transformando el hiperelipsoide en una hiperesfera, vía multiplicar el largo de los ejes por D-1/2 .

  28. La i-esima columna (fi) de la matriz [F] de componentes “loadings” es la correlación entre la matriz de datos estandarizada y la i-esima componente principal estandarizada (Zsi) Las componentes “loadings” significan la correlación entre alguna de las variables originales estandarizadas Xs con una componente principal. (7) A partir del SVD de Xs ó de R (8) Asumiendo que todas las componentes son retenidas tal que UU’=I y también U’U=I, entonces la matriz original [Xs] como producto de la matriz de “scores” estandarizados [Zs] y de la transpuesta de la matriz de “loadings” [F’]. Si sólo son retenidas r componentes Ur’Ur=I, como ya se vio,puede obtenerse la mejor aproximación de [Xs] al rengo r mediante [Ur] y [Dr] • U=Q D=2/(m-1) • valores de  P=XsU

  29. (2) RESUMIENDO: Matriz Xs con m=3 observaciones oi y n=2 variables xsj (1) Factores son los “pesos” utilizados para expresar las variables en Xs como combinación linear de las CPs Análogamente podemos escribir las observaciones como combinación linear de los “loadings” utilizando los “scores” como “pesos” Si bien (2) es una relación válida, no significa que las CPs de las observaciones sean las componentes “loadings” de las variables. La estandarización efectuada sobre las variables x no es válida para las observaciones. Si se desea obtener las CPs de las observaciones se debe transponer la matriz de entrada X tomando ahora las antiguas observaciones como variables y estandarizar sobre ellas. Al estandarizar las variables (columnas de X) las observaciones originales sufren modificaciones  la combinación lineal (2) no es una expresión válida para las observaciones originales.

  30. (Richman, 1986)

  31. Ejemplo2 Componentes “score” Componente “score” sin estandarizar estandarizadas Z = P Zs = Xs U D-1/2 Componentes “loadings” F = UD1/2 Var. Componentes Var. Componentes Var. Componentes i 2,685 0,250 0,065 Var Acum. 0,89 0,98 1,00 A pesar de U=Q (U del SVD de R y Q del SVD de Xs), hay diferencia en signo entre Z y Zs para la 1º y 2º componente debido a que u1 y u2 tienen signo contrario a q1 y q2. Por la propiedad flip flop de las componentes permite multiplicar por (-1) a u1y u2y/ó q1y q2sin alterar el resultado.

  32. Comparando del dos grupos de “scores” cada elemento de la componente j de Z dividido por el correspondiente de Zs es igual a la raíz cuadrada del autovalor correspondiente. • El “loading (i,j)” muestra la correlación de la componente principal j con la variable i • La suma al cuadrado de los “loadings” para una variable (de una fila de F), es igual a 1 si se retienen todas las componentes. Si se retienen r < min(m.n) componentes es la varianza de la variable explicada por las componentes retenidas y se llama comunalidad • La suma del cuadrado de los “loadings” de una componente es igual a la varianza explicada por dicha componente

  33. Estas propiedades se utilizan para calcular la varianza de un patrón dado por una componente principal zj tanto en caso del modo directo al que le corresponden “loadings” positivos (correlación positiva con la variable), como el inverso con “loadings” negativos (correlación negativa con la variables) El patrón dado por la 1ºCP sólo existe en modo negativo ya que las correlaciones con cualquiera de las tres variables es alta y negativa Varianza explicada por el modo directo = 0% y varianza por el modo inverso = 89% Al patrón dado por la 2ºCP le corresponden “loadings” tanto negativos como positivos  Varianza explicada por el modo directo = 4,2% y por el modo inverso = 4,1% Al patrón dado por la 3ºCP también le corresponden “loadings” positivos y negativos Varianza explicada por el modo directo = 1,4% y por el modo inverso = 1,38% Generalmente las cantidades de varianza explicadas por el modo directo difieren de las explicadas por el inverso.

  34. Análisis de Componentes Principales de campos de presión diarios de Junio de 1972 : Matriz de correlación entre los campos diarios (n=31; m= número de puntos de enrejado) Porcentaje de varianza total explicada por la 2ºCP =(  fi22 /n)*100=19,06% Porcentaje de varianza explicada por el patrón con bajas presiones en el área sombreada (“loadings” positivos; modo directo) = 8,15% Porcentaje de varianza explicada por el patrón con altas presiones en el área sombreada (“loadings” negativos; modo inverso) = 10,80%

  35. Vargas y Compagnucci 1983

  36. Compagnucci R. y Vargas (1985 a y b) Objetivo: a) Encontrar regiones homogeneas  Modo-S b) Encontrar los patrones espaciales Modo-T Datos: presión diaria de superficie de Julio 1972-1977 (186 campos) definidos por un enrejado de 81 estaciones Matriz de similitud: correlación (covarianza normalizada) Modo-S RS de 81x81  estandarizar cada serie temporal Modo-T  RT de 186x186  estandarizar cada campo diario Metodología: Método de Jacobi para obtener los autovalores y autovectores R ui = i uidonde son los autovalores yuson los autovectores F = UD1/2Xs= ZsF’  Zs = Xs U D-1/2 donde F es la matriz de factor “loadings” , U la matriz de autovectores, D la matriz diagonal de autovalores, Zs la matriz de factores “scores” y Xs la matriz de datos estandarizada

  37. Modo - T Modo - S

  38. 15,57% 14,94% 8,77% 6,79%

  39. 4,34% 4,11% 3,64% 3,56%

  40. 2,31% 2,51% 1,94% 2,07%

  41. Modo - T

  42. Compagnucci R. and M. Salles (1997):“Surface Pressure patterns during the year over southern South America” International Journal of Climatology, 17, 6, 635-654. OBJETIVO: encontrar los tipos sinópticos de superficie de 12:00 UTC mediante el ACP de los campos diarios de presión a nivel del mar METODO: ACP no-rotadas usando [R] en Modo-T (Correlación entre campos diarios) DATOS: 4383 campos diarios de 1972-1983 observados en 81 estaciones. Notar : enrejado irregular. 12 matrices una para cada mes (ej. [XEnero] de 372x81 y [REnero] de 372x372 )

  43. 2 de Mayo de 1976 16 de Julio de 1975 18 de Julio de 1978 Patrón A Patrón B En la serie de tiempo de las componente “loadings” se eligen los valores máximos y mínimos para encontrar ejemplos en la matriz de datos que se correspondan a los patrones de las CPs. Fotos de los mapas sinópticos de superficie de 21:00 UTC del Servicio Meteorológico Nacional

  44. Patrón C 6-Julio-1978 26-Julio-1980 Patrón D 30-Agosto-1972 25-Julio-1980

  45. Patrón E 27-Junio-1983 21-Junio-1981 Patrón F 28-Mayo-1976 19-Mayo-1977

More Related