miércoles, 8 de agosto de 2012

Calculando coeficientes de correlación parcial con la HP50G

Estoy repasando mis apuntes de estadística (por afición, creo que me puede ser útil en el futuro), y me he decidido a hacer un programa para mi calculadora (una maravillosa HP50G, que es el orgullo de todo ingeniero) para ayudarme a hallar coeficientes de correlación y coeficientes de correlación parcial entre variables marginales de muestras de variables n-dimensionales.

Antes de que los no iniciados a la estadística se asusten, se puede leer sobre el tema en el libro: "Elementos de estadística aplicada: Cálculo de probabilidades y teoría de variable aleatoria" de José Javier Muruzábal Irigoyen , que es el que estoy usando yo ahora mismo.

Además, se puede consultar la wikipedia, que es una auténtica fuente de sabiduría: http://en.wikipedia.org/wiki/Partial_correlation - en inglés.

El problema es el siguiente: tenemos varias listas de datos y nos interesa saber si hay alguna relación considerable entre unas variables y otras. Para ello, se calcula un coeficiente de correlación para cada pareja de variables. Cada coeficiente de correlación irá de 0 a 1, y puede ser de signo positivo o negativo si la relación entre las variables es directa o inversa.

Además, nos interesa calcular la correlación entre cada pareja de variables, pero esta vez sin tener en cuenta el efecto de las demás variables. Para ello calcularemos los coeficientes de correlación parcial, que también valdrán entre 0 y 1, y, de nuevo, su signo será positivo si la relación entre la pareja de variables es directa, y negativo si la relación es inversa.

También se pueden entender estos coeficientes como la medida en que la varianza de una variable puede explicar la varianza de otra, considerando la afección de las demás variables en el caso de los coeficientes de correlación, y eliminando el efecto de las demás variables en el caso de los coeficientes de correlación parcial.

Se entiende mucho mejor con un ejemplo (un clásico): Se mide la estatura, el coeficiente intelectual y la edad de un grupo de estudiantes. Estos son los resultados:

Altura (cm)            IQ                   Edad  

   188                    121                    30
   183                    126                    35
   170                    106                    16
   150                    101                     9
   190                    119                    26
   191                    124                    38
   160                    106                    12
   165                    107                    15
   185                    116                    20
   175                    111                    19


Lo que queremos saber es si hay alguna relación notable entre la altura y el coeficiente intelectual, o entre la altura y la edad, o entre la edad y el coeficiente intelectual. Por ejemplo... ¿Son mas listos los mas altos? ¿y los mas jóvenes?

Bueno, pues calcular estos coeficientes (el método viene en los enlaces anteriores) requiere hacer un gran número de cálculos, que es la tarea perfecta para la calculadora hp50G. He escrito el siguiente código (codificado en ASCII, pero que nadie se asuste. Al final del texto pondré un enlace para poder descargarlo directamente a la calculadora.):


<< { ∑DAT ∑mat P2mat r2mat r2mADJ } DUP PURGE { 1. 1. 1. 1. 1. } SWAP STO "MAT. COVARIANZAS" { { "MUESTRA:" "VARIABLES POR COLUMNAS" } } { 1. 0. } { } {
[[ 1. 2. ]
 [ 3. 5. ]
 [ 6. 3. ]
 [ 5. 8. ]] } INFORM
  IF NOT
  THEN KILL
  END EVAL DUP SIZE EVAL
  IF DUP 2. <
  THEN "NO HAY SUFICIENTES COLUMNAS DE DATOS PARA RALACIONAR" MSGBOX KILL
  END -> GSAMP NFILAS NCOLS
  << GSAMP DUP "DATA" ->TAG SWAP CL∑ STO∑ 1. NCOLS
    FOR J 1. NCOLS
      FOR I I J COL∑ ∑XY N∑ / ∑X N∑ / ∑Y N∑ / * -
      NEXT NCOLS COL->
    NEXT NCOLS ROW-> '∑mat' DUP PURGE STO CLLCD "∑mat - Matriz de covarianzas" MSGBOX ∑mat "∑mat" ->TAG 1. NCOLS
    FOR I I XCOL ∑X2 N∑ / ∑X N∑ / 2. ^ - v/ INV
    NEXT NCOLS ROW-> NCOLS DIAG-> DUP ∑mat SWAP * * 1. NCOLS
    FOR J 1. NCOLS
      FOR I DUP I J 2. ->LIST GET
        IF DUP 0. >
        THEN 2. ^
        ELSE 2. ^ NEG
        END I J 2. ->LIST SWAP 4. RND PUT
      NEXT
    NEXT DUP "P2mat" ->TAG SWAP 'P2mat' DUP PURGE STO "P2mat - matriz correlaciones lineales simples elementos Gr^2" MSGBOX 'r2mat' DUP PURGE NCOLS IDN SWAP STO { MATCOV ∑mat P2mat r2mat r2mADJ rADJ } ORDER 1. NCOLS 1. -
    FOR I I 1. + NCOLS
      FOR J I J 1. NCOLS
        FOR AUXV1
          IF AUXV1 I =/ AUXV1 J =/ AND
          THEN AUXV1
          END
        NEXT NCOLS
        IF DUP 3. <
        THEN DROPN KILL
        END ->LIST rAUX r2mat SWAP I J 2. ->LIST SWAP
        IF DUP 0. >
        THEN 2. ^
        ELSE 2. ^ NEG
        END 4. RND PUT 'r2mat' STO
      NEXT
    NEXT "r2mat - elemetos r^2 correlaciones parciales, por parejas sin las demas variables" MSGBOX r2mat 'r2mat' ->TAG
  >>
>>


El que quiera puede escribirlo en la calculadora, y almacenarlo en una variable que se llame RUN o MATCOV: estará listo para ser ejecutado. El programa requiere una variable auxiliar que debe llamarse rAUX, y que debe contener el siguiente código:


<< DUP HEAD SWAP TAIL DUP HEAD SWAP TAIL DUP SIZE SWAP DUP REVLIST HEAD SWAP REVLIST TAIL REVLIST "<< -> IAUX JAUX NCONT"
  IF NCONT TYPE 0. SAME  THEN NCONT ->STR +  END " COLAL COLAS
  <<
    IF NCONT" +  IF NCONT TYPE 0. SAME  THEN NCONT ->STR +  END " 1. >
    THEN IAUX COLAL COLAS + + rAUX DUP 2. ^ 1. SWAP - √ JAUX COLAL COLAS + + rAUX
DUP 2. ^ 1. SWAP - √ SWAP 3. ROLLD * 3. ROLLD * IAUX JAUX COLAS + + rAUX SWAP - SWAP / NCONT" +   IF NCONT TYPE 0. SAME   THEN NCONT ->STR +   END " 1 - 'NCONT" +
  IF NCONT TYPE 0. SAME   THEN NCONT ->STR +   END "' STO ELSE P2mat IAUX JAUX 2. 
->LIST GET P2mat IAUX COLAL 2. ->LIST GET DUP 2. ^ 1. SWAP - √ P2mat JAUX COLAL 2. 
->LIST GET DUP 3. ROLLD 2. ^ 1. SWAP - √ * 3. ROLLD * 3. ROLL SWAP - SWAP / EVAL
    END
  >>" + OBJ-> EVAL
>>


Me gusta este lenguaje (usr RPL), porque las funciones vienen ya precompiladas, y eso ofrece algunas ventajas y facilidades para usar la calculadora, pero a cambio tiene bastantes limitaciones. Me encanta sortear las limitaciones.

Este programa usa la variable rAUX de forma recursiva para calcular los coeficientes de correlación parcial, y además escribe su propio código para declarar variables dinámicamente. Es decir: si introducimos una matriz o lista de datos con 5 columnas, declarará 5 grupos de variables para poder calcular cada coeficiente. Si introducimos una matriz con 7 columnas, declarará 7 grupos de variables para calcular los coeficientes. Sencillamente sublime!

Se puede descargar el programa y las variables que genera en este enlace. Es un archivo binario, y no hay mas que ponerlo en la carpeta HOME de la calculadora. Yo uso HPConnect para conectar mi calculadora a mi iMac.

Si empleamos el programa que he escrito para calcular los coeficientes de correlación de la matriz de datos anterior, resulta:

r2)12 = 0,8275
r2)13 = 0,7484
r2)23 = 0,9364

(son los elementos (i,j) de la matriz P2mat que crea el programa)

Lo que significa que efectivamente, hay una alta correlación entre altura (1) y coeficiente intelectual (2): La varianza de la altura explica el 82,75% de la varianza en el coeficiente intelectual! - r2)12

También hay cierta correlación entre altura (1) y edad (3), pero sin embargo, la mayor correlación de todas la encontramos entre el coeficiente intelectual (2) y la edad (3).

Antes de aventurarnos a sacar conclusiones, veamos que pasa con los coeficientes de correlación parcial:

r2)12-3 = 0,2963
r2)13-2 = -0,018
r2)23-1 = 0,7251

(son los elementos (i,j) de la matriz r2mat que crea el programa) 

r2)12-3: Al eliminar el efecto de la edad (3), vemos que la correlación entre la estatura (1) y el coeficiente intelectual (2) se queda en un 0,2963, mucho menos significativa de lo que parecía antes - r2)12 -.

r2)23-1 : Aún eliminando la influencia de la estatura, se mantiene una correlación significativa entre la edad y el coeficiente intelectual.

La conclusión que yo saco se corresponde con la intuición habitual: es la edad y no la altura la que guarda cierta relación con el coeficiente intelectual, puesto que maduramos al envejecer.

Repito el enlace por si alguien está interesado en usar mi programa en su hp 50G para calcular los coeficientes de correlación: https://dl.dropbox.com/u/11312619/share/ESTADISTICA.hphttps://dl.dropbox.com/u/11312619/share/ESTADISTICA.hp

No dejéis de comentar sobre el ejercicio de ejemplo o sobre el programa!

domingo, 5 de agosto de 2012

Estadísticas de los "rasca y gana" de Caprabo

Desde principios de verano Caprabo regala unos vales "rasca y gana" cuando compras en sus establecimientos. Te dan mas vales cuanto mayor sea el importe de la compra. En los vales, hay seis casillas para rascar. Si rascas tres de ellas, y descubres tres símbolos iguales, te llevas el premio de ese símbolo.

Hay seis tipos de premios (y de símbolos!):

  • 1€ de descuento en tu próxima compra
  • 2€
  • 3€
  • 100€
  • Un bono de viajes de 1500€ (Viajes Eroski)
  • 6000€

Como da la casualidad de que estoy repasando los apuntes de estadística ahora mismo, me he animado a recopilar datos de todos los vales que consigo, con la esperanza de encontrar algún patrón.

Por cada vale, leo los dos códigos de barras con el teléfono móvil (http://neoreader.com/how-neoreader-works), y rasco todas las casillas, para luego registrarlo en una hoja de Excel que se puede descargar de aquí:



Esta una foto de uno de los boletos. En rojo están señaladas los elementos a los que hacen referencia los elementos de la hoja de Excel:

Parte frontal de un vale cualquiera.

Parte trasera del mismo vale.

Para hacer el archivo mas manejable, he usado el siguiente código:
  • x = 1€
  • y = 2€
  • z = 3€
  • u = 100€
  • v = viaje
  • w = 6000€
Es interesante ver la distribución de los premios y el análisis de frecuencias. Estas son algunas conclusiones:
  • Efectivamente, todos los vales tienen premio.
  • Aún sabiendo que al menos, va a haber tres casillas iguales, la posibilidad de conseguir el premio rascando tres círculos al azar es de sólo un 5% (1 de cada 20).
  • El "nº Eroski" y el "nº Caprabo" son siempre los mismos, e identifican el premio que contiene el vale.
  • El número que identifica la distribución de los premios de cada vale es el que he denominado "secreto". Como es obvio, no nos permiten rascar el nº secreto si queremos reclamar el premio. Por ejemplo, Una papeleta con el nº secreto M4G0 siempre tiene:
    • en el azul:         100€      (u)
    • en el negro:       6000€   (w)
    • en el rojo:          1€         (x)
    • en el amarillo:   1€         (x)
    • en el fucsia:       avión     (v) 
    • en el verde:       1€         (x)
  • Los premios menores (1€, 2€ y 3€) parecen estar aleatoriamente distribuidos entre los diferentes colores, sin embargo, parece haber una mayoría de premios avión en la casilla fucsia y de premios de 6000€ en la negra.
  • Hasta ahora, no me ha tocado ningún vale que tuviera tres casillas iguales de 100€, avión, o 6000€. Me han tocado un 76% de los vales con premio de 1€, un 16% con 2€ y un 8% con 3€.
  • Por otra parte, es perfectamente normal, puesto que según la parte de atrás del vale, hay millones de vales, y sólo 1000 tarjetas regalo de 100€, 200 regalos de avión y 50 de 6000€, y mi muestra (por ahora) es de sólo 25 papeletas.
Desde aquí, animo a todos los visitantes del blog a bajarse la hoja de Excel que mantengo, y a completarla con los vales que consigan. Podéis dejarme vuestras valoraciones y aportaciones en los comentarios :-)