Contenidos de «Complementos de Estadística»

• Descargar el último R desde https://cran.r-project.org/bin/windows/base
• Interfases opcionales:
  • el paquete Rcmdr (commander); instalarlo desde el menú Paquetes del RGui o escribiendo en la R Console la siguiente orden:

    install.packages ("Rcmdr")
    

  • RStudio
  • Emacs Speaks Statistics desde Emacs para Windows

• apt install r-cran-rcmdr rkward ess

• Descargar el instalador de https://www.python.org/downloads/release
• Una vez instalado, ejecutar en el terminal cmd.exe (quizá haya que pulsar intro alguna vez):

  py -m pip install --user --upgrade pip
  py -m pip install --user numpy scipy pandas statsmodels matplotlib sklearn 
  

apt install python3-pandas python3-dill python3-sklearn # módulos necesarios para tratar datos
apt install python-numpy-doc python-scipi-doc python-pandas-doc python-sklearn-doc # documentación

El editor usado en el aula P3:

apt install idle-python3.5

• Menú Datos, submenú Conjunto de datos en paquetes.

• Menú Datos, opción Cargar conjunto de datos.
• Habitualmente, dichos ficheros tienen extensión .rda ó .RData.
• Pueden contener
  • un solo conjunto de datos (cuyo nombre puede coincidir o no con el nombre del fichero)
  • múltiples conjuntos de datos, además de variables, funciones, etc. (entorno de trabajo de R)

• Menú Datos, submenú Importar datos.
• De especial interés es la importación de datos de texto (CSV)
  • ¿campos delimitados o de anchura fija?
  • codificación: utf-8 (unicode) o latin-1 (iso-8859-1, -15)
• Ejercicio: importar los siguientes conjuntos de datos
  • vacas.csv
  • vacas.txt
  • los participantes del campus virtual
  • vacas.blancos1
  • vacas.blancos2

Para encontrar ayuda sobre un comando de R, por ejemplo rep, puedes usar

?rep

Para algunos comandos con nombres especiales, hay que usar comillas; por ejemplo,

?"["
?"if"

Si sólo conocemos parte del nombre del comando (p.ej. norm), usamos

apropos("norm")

Si queremos encontrar información relacionada con un concepto, podemos usar:

??variance
help.search("variance")

• https://cran.r-project.org/doc/contrib/R-intro-1.1.0-espanol.1.pdf
• «Spanish» en https://cran.r-project.org/other-docs.html

help("variance")
help("modules variance")

• Dive into Python 3
• Python Data Science Handbook

• condicionales

  Generemos un vector de temperaturas:

  t <- runif (15, 10, 40)
  

  Para combinar comparaciones vectorialmente, usamos ifelse y un operador «lanzado» («&» = «y»; «|» = «ó»):

  ifelse (t>20 & t<30, "agradable", "desagradable")
  

  Para combinar comparaciones de forma que esperamos un único booleano, usamos operadores «cortocircuito» («&&» = «y»; «||» = «ó»):

  n=10; if (is.numeric(n) && length(n)==1 && n>0) sqrt(n) else "no puedo"
  

• bloques de instrucciones mediante llaves «{}»; se usan en contextos donde tiene que ir más de una instrucción:

  if (n > 5) {print("demasiado grande"); n<-5} else {print("demasiado pequeño"); n<-n+1}
  

  El valor devuelto por un bloque es el valor devuelto por su última instrucción.

• bucles

  for (i in 1:10) cat ("El cuadrado de", i, "es", i^2, "\n")
  

  Detrás de «in» puede ir cualquier vector o lista; por ejemplo

  for (i in letters) {cat (i); cat (toupper (i)); cat ("-")}; cat ("\n")
  for (i in iris) print (summary (i))
  

• bucles con acumulador

  Se define una variable (acumulador) antes del bucle para que recoja los resultados:

  suma <- 0; for (i in 1:10) suma <- suma+i; suma
  producto <- 1; for (i in 1:10) producto <- producto*i; producto
  bector <- c(); for (i in 1:10) bector <- c(bector,i^2); bector
  lista <- list(); for (i in 1:10) lista <- c(lista,list(1:i)); lista
  

• condicionales

  Para combinar comparaciones de forma que esperamos un único booleano, usamos operadores «cortocircuito» («and» = «y»; «or» = «ó»):

  import math
  n=10
  if ((isinstance(n,int) or isinstance(n,float)) and n>0):
     math.sqrt(n) 
  else:
     "no puedo"
  

  Generemos un vector de temperaturas:

  import scipy.stats as st
  t = st.uniform.rvs (10, 40, 15)
  

  Para combinar comparaciones vectorialmente, usamos if así:

  ["agradable" if ti>20 and ti<30 else "desagradable" for ti in t]
  

• nótense los bloques de instrucciones tras dos puntos y con sangría:

  if n > 5:
      print("demasiado grande")
      n = 5
  else:
      print("demasiado pequeño")
      n++
  

• bucles

  for i in range(10):
    print ("El cuadrado de", i, "es", i**2)
  
  mtcars = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/mtcars.csv")
  for v in mtcars:
    mtcars[v].describe()
  

• bucles con acumulador

  Se define una variable (acumulador) antes del bucle para que recoja los resultados:

  • numéricos

  suma = 0
  for i in range(1,11):
    suma += i
  

  producto = 1
  for i in range(1,11):
    producto *= i
  

  • recolección

  vector = []
  for i in range(1,11):
     vector.append (i**2)
  

  lista = []
  for i in range(10):
     lista.append (list (range (i)))
  

• con nombre

  f <- function(x) x^2
  f(15)
  

• anónimas (útiles para definiciones al vuelo, como se ve en la sección siguiente)

  (function(x) x^2) (15)
  

• pueden tener varios argumentos

  var.ponderada <- function (x, pesos)     # varianza ponderada
  {  mp <- weighted.mean (x, pesos)        # media ponderada
     weighted.mean ((x - mp) ^ 2, pesos) }
  

• cada argumento puede tener un valor por omisión

  var.ponderada <- function (x, pesos = rep(1,length(x)))
  {  mp <- weighted.mean (x, pesos)
     weighted.mean ((x - mp) ^ 2, pesos) }
  

• pueden pasar sus argumentos a otras funciones mediante «…»

  var.ponderada <- function (x, pesos, ...)
  {  mp <- weighted.mean (x, pesos, ...)
     weighted.mean ((x - mp) ^ 2, pesos, ...) }
  

  Por ejemplo:

  muestra <- c (1, 20, NA, 5); pesos <- c (1,2,3,4)
  weighted.mean (muestra, pesos); 
  weighted.mean (muestra, pesos, na.rm=TRUE)
  var.ponderada (muestra, pesos)                  # produce NA
  var.ponderada (muestra, pesos, na.rm=TRUE)      # pasa «na.rm=TRUE» a weighted.mean
  

• con nombre

  def f (x):
    return x**2
  f(15)
  

• anónimas (útiles para definiciones al vuelo, como se ve en la sección siguiente)

  (lambda x: x**2) (15)
  

  sólo pueden contener una expresión;

• pueden tener varios argumentos

  def media_ponderada (x, pesos):
     suma = 0
     for i in range (len (x)):
        suma += x[i] * pesos[i]
     return suma / sum(pesos)
  

• cada argumento puede tener un valor por omisión

  def media_ponderada (x, pesos = None):
     suma = 0
     for i in range (len (x)):
        suma += x[i] * (pesos[i] if pesos else 1)
     return suma / (sum(pesos) if pesos else len(x))
  

• pueden pasar sus argumentos a otras funciones mediante un argumento con «*»

  def var_ponderada (x, *otros):
    mp        = media_ponderada (x, *otros)
    cuadrados = [(xi-mp)**2 for xi in x]
    return media_ponderada (cuadrados, *otros)
  

  Por ejemplo:

  muestra = [1, 20, 100, 5]; pesos = [1,2,3,4]
  media_ponderada (muestra)
  media_ponderada (muestra, pesos)
  var_ponderada (muestra)                  # pasa *otros como vacío
  var_ponderada (muestra, pesos)           # pasa *otros como «pesos»
  

• sobre secuencias (vectores o listas)
  • sapply (con s de simple) se usa para generar un vector (o una matriz, si el resultado de la función aplicada es vectorial)

    sapply (1:10, sqrt)                       # aplicar una función con nombre
    sapply (1:10, function(x) x**2)           # aplicar una función anónima
    sapply (1:10, function(x) c(x^2,x^3))     # el resultado es una matriz
    

  • lapply (con l de lista) devuelve una lista

    lapply (1:10, function(x) list (x, list(letters[i],month.name[i])))
    

• apply se usa para aplicar una función a filas o columnas de una matriz o dataframe

  m <- matrix (1:12, 3)
  apply (m, 1, sum)                 # aplicar por filas
  apply (m, 2, sum)                 # aplicar por columnas
  

• para aplicar por grupos:

  prod <- 1:20; sexo <- rep(c("M","H"),each=10)
  tapply (prod, sexo, mean)           # por niveles de "sexo"
  by (prod, sexo, mean)               # ídem
  

• la más general (mapeo)

  mapply (rep, 1:4, 4:1)              # aplica funciones con varios argumentos
  

• un comando muy habitual en simulaciones:

  replicate (10000, {m <- rnorm(5); mean(m)})
  

  Funciona parecido a

  sapply (1:10000, function (i) {m <- rnorm(5); mean(m)})
  

  donde la variable i es irrisoria; su valor no se usa.

• trasformación por descripción (mapeo)

  from math import sqrt
  [sqrt(i) for i in range(10)]
  list (map (sqrt, range(10)))               # lo mismo
  list (map (lambda x: x**2, range(10)))     # con función anónima
  

  otra versión de la media ponderada anteriormente definida:

  def media_ponderada (x, pesos):
     return  sum (map (lambda xi,pi: xi*pi, x, pesos)) / sum(pesos)
  

  otra versión con argumento opcional:

  def media_ponderada (x, pesos = None):
     return  sum (map (lambda xi,pi: xi*pi, x, [1]*len(x) if pesos is None else pesos)) / (sum(pesos) if pesos else len(x))
  

• aplicar una función a filas o columnas de un dataframe

  from sklearn import datasets
  iris = datasets.load_iris()        # datos de ejemplo
  import pandas as pd
  I = pd.DataFrame(iris.data, columns=iris.feature_names)
  I.apply(sum)                       # aplicar «sum» por columnas; equivale a I.sum()
  I.apply(sum,axis=0)                # ídem
  I.apply(sum,axis=1)                # aplicar «sum» por filas; equivale a I.sum(axis=1)
  

• para aplicar por grupos:

  I["especie"] = pd.Series(pd.Categorical.from_codes(codes=iris.target,categories=iris.target_names))
  I.groupby("especie").apply(lambda x: x.iloc[:,0:4].median()) # 0:4 para evitar "especie"
  

(Lo mismo en Rcmdr, Python.)

• Guardar todo el entorno de trabajo (funciones, dataframes…)

  save.image ()                         # lo guarda en ./.RData (el fichero cargado al arrancar R)
  save.image (file = "~/trabajo/mi-sesión.rda") # en fichero arbitrario
  

• Guardar uno o varios objetos en formato nativo de R

  save (mi.dataframe, mi.función, file = "~/trabajo/mis-objetos.rda")
  

• Guardar un dataframe como una tabla de texto (CSV)

  write.table (mtcars, "/tmp/mtcars1.csv") # sep=" " dec="." sin campo de rownames en cabecera
  write.csv   (mtcars, "/tmp/mtcars2.csv") # sep="," dec="."
  write.csv2  (mtcars, "/tmp/mtcars3.csv") # sep=";" dec=","
  

(Lo mismo en R, Rcmdr.)

• Guardar entorno de trabajo

  import dill
  dill.dump_session ("/tmp/miTrabajo.pkl") # recuperar con dill.load_session
  

• Guardar un objeto en formato nativo de Python

  misDatos.to_pickle ("/tmp/misDatos.pkl")
  

• Guardar un dataframe como una tabla de texto (CSV)

  misDatos.to_csv ("/tmp/misDatos.csv")
  

(Lo mismo en Rcmdr, Python.)

D    <- mtcars                # dataframe de ejemplo
D$am <- factor (mtcars$am)    # una cualitativa; también: D$am <- factor(D$am, labels=c("auto","manual"))
summary (D)                   # resumen del dataframe
T <- table (D$cyl)            # frecuencias absolutas
prop.table (T)                # frecuencias relativas
summary (D$hp)                # media y cinco cuartiles
mean (D$hp)                   # media
var (D$hp)                    # varianza
sd (D$hp)                     # desviación típica
median (D$hp)                 # mediana
quantile (D$hp, 0.95)         # percentil 95
quantile (D$hp, 0:4/4)        # cinco cuartiles
IRQ (D$hp)                    # recorrido intercuartílico

En muchas funciones se puede incluir la opción na.rm=TRUE para omitir los ausentes:

mean (c (1:10, NA))
mean (c (1:10, NA), na.rm=TRUE)

Gráficas:

barplot (table (D$cyl))   # diagrama de barras
pie (table (D$am))        # diagramo de sectores
hist (D$hp)               # histograma
boxplot (hp ~ am, D)      # diagrama de cajas
plot (mpg ~ hp, D)        # diagrama de dispersión (nube de puntos)
pdf ("fichero.pdf")       # guarda gráficas subsiguientes en 1 fichero PDF (vectorial)
for (i in 1:ncol(mtcars)) {
  if (is.numeric (mtcars[[i]]))
    hist (mtcars[[i]])
  else 
    barplot (table (mtcars[[i]])) }
dev.off ()                # cierra el PDF
png ("fichero%d.png")     # guarda gráficas subsiguientes en ficheros PNG (pixelados)
for (i in 1:3)
  boxplot (mtcars[[i]] ~ mtcars$am)
dev.off ()                # cierra el PNG

Ejercicio:

• Escribir un programa en R que realice gráficas adecuadas para las variables del conjunto de datos «Chile» del paquete «car».

  data (Chile, package="car")
  Chile$education <- factor (Chile$education, levels=c("P","S","PS"))
  

(Lo mismo en R, Rcmdr.)

Con los datos clásicos de flores iris:

import pandas as pd
i = pd.read_csv ("datos/iris.csv")
i.describe()                           # numérico
i["Species"].describe()                # moda
iSn = i['Species'].value_counts()      # recuentos
iSn / iSn.sum()                        # frecuencias relativas
i["Petal.Length"].mean()               # media
i["Petal.Length"].var()                # varianza
i["Petal.Length"].std()                # desviación típica
i["Petal.Length"].median()             # mediana
i["Petal.Length"].quantile(0.5)        # ídem
i["Petal.Length"].quantile([.05,.95])  # percentiles 5 y 95
import numpy as np                     # para «nan»
pd.Series([1,2,np.nan,4]).mean()       # excluye los ausentes

Gráficas:

import matplotlib.pyplot as plt
i["Petal.Length"].plot(kind="hist")                # histograma
plt.show()                                         # mostrar en pantalla
i["Petal.Length"].value_counts().plot(kind="bar")  # diagrama de barras
plt.savefig("barras.png")                          # guardar en fichero
i["Petal.Length"].value_counts().plot(kind="pie")  # diagrama de sectores
i["Petal.Length"].plot(kind="box")                 # diagrama de caja
i.boxplot("Petal.Length","Species")                # cajas por grupo

(Lo mismo en Python, Rcmdr.)

• Crear nueva variable de un dataframe:

  mtcars$consumo <- (3.8 / mtcars$mpg * 1.6) * 100
  

• Recodificar:

  mtcars$am <- factor (mtcars$am, labels=c("auto","manual"))
  mtcars$consumo3 <- cut (mtcars$consumo, c(-Inf, 10, 20, 30, Inf), c("muy poco", "poco", "mucho", "demasiado"))
  

• Reordenar niveles de factor:

  data (Chile, package="car")
  Chile$educación <- factor (Chile$education, levels = c("P","S","PS")) # P=primaria, S=secundaria, PS=postsecundaria
  

• Filtrar

  edades.mujeres <- Chile$age [Chile$sex == "F"]
  mtcars4v <- mtcars [mtcars$cyl==4 & mtcars$vs==0, ]
  ?==                           # ayuda sobre < > == != 
  ?&                            # ayuda sobre ! & |
  

(Lo mismo en R, Rcmdr.)

import pandas as pd
Chile = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/Chile.csv")

• Crear nueva variable en un dataframe:

  Chile.euro = 1.2 * Chile.income
  

• Recodificar:

  Chile.sexo = Chile.sex.astype("category")
  Chile.sexo.values.categories = ["mujer", "varón"]
  Chile.clase = pd.cut (Chile.income, [-np.inf,30000,100000,np.inf], labels=["baja","media","alta"])
  

• Reordenar niveles de categórica:

  Chile["educación"] = Chile["education"].astype("category",categories=["P","S","PS"])
  

• Filtrar:

  edades_mujeres = Chile.age[Chile.sex=="F"]
  varones_jovenes = Chile[(Chile.sex=="M") & (Chile.age<35)]
  

Más en la documentación de Pandas.

(Lo mismo en Python, Rcmdr.)

• letra inicial
  • r para generar aleatorios
  • p para la función de distribución \(F_X\,(t) = \Pr\,[X\leq t]\)
  • d para la función de densidad de variables continuas o la función de probabilidad (o de cuantía, o de masa) de variables discretas \(f_X\,(t) = \Pr\,[X=t]\)
  • q para la función cuantil
• abreviatura del nombre de la familia de la distribución
  • unif uniforme
  • norm gausiana
  • exp exponencial
  • binom binomial
  • pois puasón
  • t de Estiúdent (Student)
  • chisq \(\chi^2\) (ji-cuadrado)

runif (10, 0, 1)                   # genera 10 aleatorios con distribución uniforme entre 0 y 1
pnorm (150, 170, 15)               # probabilidad de que una gausiana N(170;15) sea menor que 150
plot (x <- seq(0,5,.1), dexp(x,1)) # curva de densidad de una exponencial con media 1
qpois (0.5, 1.7)                   # mediana de una distribución de Puasón con media 1,7

(Lo mismo en R, Rcmdr.)

• métodos
  • rvs para generar aleatorios
  • cdf función de distribución (cumulative density function)
  • pdf función de densidad (probability density function)
  • pmf función masa de probabilidad (probability mass function)
  • mean, var, std, median - media, varianza, desviación típica, mediana
  • ppf función cuantil (percent point function)
• distribuciones
  • uniform uniforme
  • norm gausiana
  • expon exponencial
  • binom binomial
  • pois puasón
  • t de Estiúdent (Student)
  • chi2 ji-cuadrado (\(\chi^2\))

import scipy.stats as st
st.uniform.rvs (2, 0.1, 30)     # 30 valores uniformes entre 2 y 2+0,1
st.norm.cdf (200, 170, 10)      # función de distribución en 200 de una N(170;10)
import matplotlib.pyplot as plt # función de densidad (de la docstring de expon)
fig, ax = plt.subplots(1, 1)
x = np.linspace(expon.ppf(0.01),
                expon.ppf(0.99), 100)
ax.plot(x, expon.pdf(x),
        'r-', lw=5, alpha=0.6, label='exponencial')
## también se pueden fijar previamente los parámetros:
v = st.binom (10, .3)           # binomial n=10, p=0,3
v.mean(), v.var(), v.std()      # media, varianza, desviación típica
v = st.poisson (100)            # puasón con media 100
v.ppf([.05,.95])                # percentiles 5 y 95

(Lo mismo en Python, Rcmdr.)

?cars                            # velocidad y distancia de frenada
library (car)                    # para datos Chile
## contraste de gausianidad
shapiro.test (cars$speed)        # comprobar mediante stem o hist
shapiro.test (Chile$statusquo)
shapiro.test (Chile$age)
## información del objeto "contraste"
shapiro.test (cars$speed) -> contraste
names (contraste)                # "statistic" "p.value"   "method"    "data.name"
contraste$p.value                # nivel crítico o p-valor
## simulación
n <- 100                                                        # tamaño muestral
m <- rnorm(n); stem(m); shapiro.test(m)$p.value                 # relación entre forma y p-valor
mean (replicate (10000, shapiro.test(rnorm(n))$p.value) < 0.05) # proporción de rechazos
## contraste sobre la media
t.test (cars$speed, mu=15)                        # bilateral
t.test (cars$speed, mu=15, alternative="less")    # H0: mu=15; H1: mu<15
t.test (cars$speed, mu=15, alternative="greater") # H0: mu=15; H1: mu>15
## contraste sobre proporción
summary (Chile$sex)
?prop.test
prop.test (sum(Chile$sex=="F"), nrow(Chile), 0.5) # H0: p = 1:2; H1: p =/= 1:2
## contraste sobre la mediana
wilcox.test (cars$speed, mu=15)                        # bilateral
wilcox.test (cars$speed, mu=15, alternative="less")    # H0: mu=15; H1: mu<15
wilcox.test (cars$speed, mu=15, alternative="greater") # H0: mu=15; H1: mu>15
M  <- rexp (25)                                        # exponencial con media=1
Me <- qexp (.5)                                        # mediana teórica = 0,693.147.2
wilcox.test (M, mu=Me)                                 # poco recomendable porque exige simetría
mean (replicate (100000, {M <- rexp(25); wilcox.test(M,mu=Me)$p.value<.05})) # 11% en vez de 5%
prop.test (sum(M>Me), 25, .5)                          # prueba de signos: poco potente, pero cumple
mean (replicate (100000, {M <- rexp(25); prop.test(sum(M>Me),25,.5)$p.value<.05})) # <5%
mean (replicate (100000, {M <- rexp(25); t.test(M,mu=1)$p.value<.05})) # 7,5% (pero con media, no mediana)
## contraste sobre dos proporciones
## ¿es significativamente mayor la proporción de mujeres que votan Y?
## p = proporción que votan Y; H0: pF=pM; H1: pF>pM
mujeres <- Chile [Chile$sex=="F" & !is.na(Chile$vote),]
varones <- Chile [Chile$sex=="M" & !is.na(Chile$vote),]
prop.test (c (sum (mujeres$vote=="Y"), sum (varones$vote=="Y")),
           c (nrow (mujeres),          nrow (varones)),
           alternative = "greater")
## contraste de igualdad de varianzas (auxiliar para el de igualdad de medias)
? iris                  # comparar varianza de Petal.Length en «versicolor» y «virginica»
with (iris, plot (Petal.Length, Petal.Width, col=Species))
boxplot (Petal.Length ~ Species, iris)
with (iris, by (Petal.Length, Petal.Species, shapiro.test))
PLversicolor <- iris [iris$Species=="versicolor", "Petal.Length"]
PLvirginica  <- iris [iris$Species=="virginica",  "Petal.Length"]
var.test (PLversicolor, PLvirginica)      # contraste de igualdad de varianzas
var.test      (Petal.Length ~ Species, iris, subset = Species %in% c("versicolor","virginica")) # lo mismo
bartlett.test (Petal.Length ~ Species, iris, subset = Species %in% c("versicolor","virginica")) # otro similar
## contraste de igualdad de medias (muestras independientes)
t.test (Petal.Length ~ Species, iris, subset = Species %in% c("versicolor","virginica")) # estimando dos varianzas distintas
t.test (Petal.Length ~ Species, iris, subset = Species %in% c("versicolor","virginica"),
        var.equal = TRUE)                                                                # casi no cambia el p-valor
## contraste de igualdad de medias (muestras pareadas)
? ChickWeight # ¿Hay incremento significativo de peso a partir del día 20?
boxplot (weight ~ Time, ChickWeight)
cw <- ChickWeight [ChickWeight $ Time %in% c (20, 21), ]
t.test (weight ~ Time, cw, alternative="less") # como si fueran independientes
cw20       <- cw [cw$Time == 20,]; rownames(cw20) <- cw20$Chick
cw21       <- cw [cw$Time == 21,]; rownames(cw21) <- cw21$Chick
ambos      <- intersect (rownames(cw20), rownames(cw21))
incremento <- cw21[ambos,"weight"] - cw20[ambos,"weight"]
shapiro.test (incremento) # innecesario por el tamaño de muestra (45)
t.test (cw20[ambos,"weight"], cw21[ambos,"weight"], paired=TRUE, alternative="less")
t.test (incremento, mu=0, alternative="greater") # equivale al anterior
## contraste sobre la mediana de la diferencia (muestras independientes)
## comparar Petal.Width en «setosa» y «virginica» dentro del conjunto de datos «iris»
?iris
shapiro.test (iris [iris$Species == "setosa",    "Petal.Width"])
shapiro.test (iris [iris$Species == "versicolor","Petal.Width"])
t.test      (Petal.Width ~ Species, iris, subset = Species %in% c("setosa","versicolor")) # suficiente tamaño muestral
wilcox.test (Petal.Width ~ Species, iris, subset = Species %in% c("setosa","versicolor")) # resultado parecido
## contraste sobre la mediana de la diferencia (muestras pareadas)
?sleep            # ¿hay diferencias significativas en sueño extra según somnífero?
boxplot (extra ~ group, sleep)
extra1 <- sleep$extra [sleep$group == '1']
extra2 <- sleep$extra [sleep$group == '2']
dif    <- extra1 - extra2
shapiro.test (dif)
wilcox.test (extra1, extra2, paired=TRUE)       # asegurarse de que el orden se corresponde en ambos vectores
wilcox.test (sleep ~ group, sleep, paired=TRUE) # otra forma de expresarlo
wilcox.test (dif)                               # es equivalente a un contraste de Wilcoxon de 1 muestra
t.test (extra1, extra2, paired=TRUE)            # no muy distinto

• Otro ejemplo:
  • Comparar las dos secuencias «x» y «y» de seudoaleatorios de cierto generador:

    ?randu                                             # comparar «x» con «y»
    dif <- randu$x - randu$y
    shapiro.test (dif)                                 # irrelevante por el tamaño muestral, en principio...
    wilcox.test (randu$x, randu$y, paired=TRUE)
    t.test (randu$x, randu$y, paired=TRUE)             # ...pero allende el umbral del 5%
    wilcox.test (runif(400), runif(400), paired=TRUE)  # comparar con seudoaleatorios decentes
    

(Lo mismo en R, Rcmdr.)

import numpy  as np                # para nan, isnan
import pandas as pd                # para dataframes
from scipy import stats            # para contrastes
import matplotlib.pyplot as plt    # para histograma
cars  = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/cars.csv") 
Chile = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/Chile.csv")
## contraste de gausianidad
stats.shapiro (cars.speed)
stats.shapiro (cars.speed) [1]                     # p-valor
stats.shapiro (Chile.age [~ np.isnan (Chile.age)]) # hay NA
## simulación
n = 100                                                     # tamaño muestral
## relación entre forma y p-valor
m = stats.norm.rvs(0,1,n); m.plot(kind="hist"); plt.show(); stats.shapiro(m)[1]
(pd.Series([stats.shapiro(stats.norm.rvs(0,1,n))[1] for i in range(10)]) < 0.05) . mean() # proporción de rechazos
## contraste sobre la media
pvalor2 = stats.ttest_1samp (cars.speed, 15) . pvalue                       # bilateral
## ya que mediamuestral = 15,4 > 15 entonces
pvalor1 = pvalor2 / 2                                                       # unilateral H0: mu=15; H1: mu>15
1 - pvalor1                                                                 # unilateral H0: mu=15; H1: mu<15
## contraste sobre proporción H0: p = 1:2; H1: p =/= 1:2
Chile.sex.value_counts()
n = Chile.sex.count()
p = Chile.sex.value_counts(True)['F']
from math import sqrt
2 * (1 - stats.norm.cdf (abs ((p - 0.5) / sqrt (0.25/n))))  # p-valor bilateral
stats.chisquare ([1379,1321]) . pvalue                      # otra posibilidad: usar contraste ji-2 equivalente
## contraste sobre la mediana
M  = stats.expon.rvs (0, 1, 25)                        # exponencial con media=1
Me = stats.expon.ppf (.5)                              # mediana teórica = 0,693.147.2
stats.wilcoxon (M - Me)                                # poco recomendable porque exige simetría
nM = sum (M-Me > 0); stats.chisquare ([nM, len(M)-nM]) # prueba de signos: poco potente
## contraste sobre dos proporciones
## ¿es significativamente mayor la proporción de mujeres que votan Y?
## p = proporción que votan Y; H0: pF=pM; H1: pF>pM
tabla = Chile.groupby(['sex','vote']).size().unstack()                       # tabla de contingencia
stats.chi2_contingency(tabla)[1]                                             # P-valor de contraste ji-2 de homogeneidad
mujeres = Chile [(Chile.sex=="F") & ~pd.isnull(Chile.vote)]
varones = Chile [(Chile.sex=="M") & ~pd.isnull(Chile.vote)]
pvalor2 = stats.chi2_contingency([[sum(varones.vote=="Y"),sum(varones.vote!="Y")],
                                  [sum(mujeres.vote=="Y"),sum(mujeres.vote!="Y")]]) [1]
pF = sum(mujeres.vote=="Y")/len(mujeres) # 37%
pM = sum(varones.vote=="Y")/len(varones) # 32%
pvalor1 = pvalor2 / 2                    # porque pF > pM
## contraste de igualdad de varianzas (auxiliar para el de igualdad de medias)
## ¿hay diferencias significativas en sueño extra entre los dos grupos?
sleep = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/sleep.csv")
stats.shapiro (sleep.extra [sleep.group == 1])
stats.shapiro (sleep.extra [sleep.group == 2])
stats.bartlett (sleep.extra [sleep.group == 1], sleep.extra [sleep.group == 2])
## contraste de igualdad de medias (muestras independientes)
stats.ttest_ind (sleep.extra [sleep.group == 1], sleep.extra [sleep.group == 2]) # suponiendo varianzas iguales
stats.ttest_ind (sleep.extra [sleep.group == 1], sleep.extra [sleep.group == 2], equal_var=False)
## contraste de igualdad de medias (muestras pareadas)
## ¿hay incremento significativo de peso a partir del día 20?
CW = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/ChickWeight.csv")
CW.boxplot ("weight", "Time")
plt.show()
CW.weight[CW.Time==20].mean()                                              # menor
CW.weight[CW.Time==21].mean()                                              # mayor
stats.ttest_ind (CW.weight[CW.Time==20], CW.weight[CW.Time==21]).pvalue/2  # como si fueran independientes
cw = CW[CW.Time.isin([20,21])]
cw20 = cw [cw.Time == 20]; cw20.index = cw20.Chick
cw21 = cw [cw.Time == 21]; cw21.index = cw21.Chick
ambos = cw20.index.intersection(cw21.index)
stats.shapiro (cw20.loc[ambos,"weight"] - cw21.loc[ambos,"weight"]) # innecesario por el tamaño de muestra (45)
stats.ttest_rel (cw20.loc[ambos,"weight"], cw21.loc[ambos,"weight"]).pvalue / 2
## contraste sobre la mediana de la diferencia (muestras independientes)
## comparar Petal.Width en «setosa» y «virginica» dentro del conjunto de datos «iris»
iris = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/iris.csv")
stats.shapiro (iris.loc[iris.Species == "setosa",    "Petal.Width"])
stats.shapiro (iris.loc[iris.Species == "virginica", "Petal.Width"])
stats.ttest_ind (iris.loc[iris.Species=="setosa","Petal.Width"], 
                 iris.loc[iris.Species=="virginica","Petal.Width"], equal_var=False)
stats.mannwhitneyu (iris.loc[iris.Species=="setosa","Petal.Width"], 
                    iris.loc[iris.Species=="virginica","Petal.Width"], alternative='two-sided') # parecido
## contraste sobre la mediana de la diferencia (muestras pareadas)
randu = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/randu.csv") # comparar «x» con «y»
dif = randu.x - randu.y
stats.shapiro (dif)                   # irrelevante por el tamaño muestral
stats.wilcoxon (dif)                 
stats.wilcoxon  (randu.x, randu.y)    # otra forma de hacer lo mismo
stats.ttest_rel (randu.x, randu.y)    # parecido
stats.ttest_1samp (dif, 0)            # lo mismo

(Lo mismo en Python, Rcmdr.)

library(car) # para Chile

### contraste ji-2

## ¿hay relación entre región y sexo?
table (Chile$region, Chile$sex)
prop.table(table (Chile$region, Chile$sex), 1)
chisq.test (table (Chile$region, Chile$sex))
## ¿hay relación entre educación y sexo?
## ¿hay relación entre sexo y voto?
## ¿hay relación entre educación y voto?

### contraste de correlación lineal

## ¿hay relación entre población, edad, ingresos y statusquo?
variables <- c("population", "age", "income", "statusquo") # names (which (sapply (Chile, is.numeric)))
cor    (Chile[variables], use="pair")
sapply (Chile[variables], shapiro.test)                    # no gausianas => mejor Spearman
plot (Chile$population, Chile$age)
cor.test (Chile$population, Chile$age)                     # por omisión, Pearson (exige gausianidad)
cor.test (Chile$population, Chile$age, method="spearman")
plot (Chile$population, Chile$income)
plot (jitter(Chile$population), jitter(Chile$income))      # jitter ayuda para representar discretas con pocos niveles
## etc.

(Lo mismo en R, Rcmdr.)

import numpy as np                 # para isnan
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt    # para scatter
Chile = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/Chile.csv")

### contraste ji-2

## ¿hay relación entre región y sexo?
tabla = Chile.groupby(["region","sex"]).size().unstack()
tabla = pd.crosstab (Chile.region, Chile.sex)               # otra forma
pd.crosstab (Chile.region, Chile.sex, normalize="index")    # index=filas columns=columnas
stats.chi2_contingency (tabla)
## ¿hay relación entre educación y sexo?
## ¿hay relación entre sexo y voto?
## ¿hay relación entre educación y voto?

### contraste de correlación lineal

## ¿hay relación entre población, edad, ingresos y statusquo?
Chile.corr()
variables = ["population", "age", "income", "statusquo"]
Chile[variables].apply(stats.shapiro)                                  # pero hay «nan»es
Chile[variables].apply(lambda x:stats.shapiro(x[~np.isnan(x)]))        # no gausianas => mejor Spearman
plt.scatter (Chile.population, Chile.age)
def cor_test (x, y):
    xsinnan = ~ np.isnan (x)
    ysinnan = ~ np.isnan (y)
    ambos   = xsinnan & ysinnan
    return stats.pearsonr (x[ambos], y[ambos])
cor_test (Chile.population, Chile.age)                                 # pearson
stats.spearmanr (Chile.population, Chile.age)                          # spearman
## etc.

(Lo mismo en Python, Rcmdr.)

cor (mtcars)                       # buscamos predictor para «mpg»
sort (abs (cor(mtcars)[,"mpg"]))   # el peso parece el mejor predictor lineal
plot (mpg ~ wt, mtcars)
modelo <- lm (mpg ~ wt, mtcars)
abline (modelo, col=2, lwd=3)      # añade recta; lwd = linewidth = anchura
summary (modelo)                   # coeficientes, R2, p-valor, estima de sigma
## p-valor del modelo = p-valor de la pendiente
predict (modelo, data.frame(wt=3)) # para un coche de 3.000 libras
coef(modelo) %*% c(1,3)            # lo mismo; 1 multiplica al intercepto

(Lo mismo en R, Rcmdr.)

import pandas as pd
mtcars = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/mtcars.csv")
mtcars.corr()
abs(mtcars.corr().mpg).sort_values() # «wt» parece el mejor predictor lineal
import matplotlib.pyplot as plt
x = mtcars.wt;              y = mtcars.mpg
x = x.values.reshape(-1,1); y = y.values.reshape(-1,1)
plt.scatter (x, y)
plt.show()
from sklearn import linear_model
regre = linear_model.LinearRegression()
regre.fit (x, y)
plt.scatter (x, y)
plt.plot (x, regre.predict (x), color='red', linewidth=3)
plt.show()
regre.predict(3)                     # predicción para 3.000 libras
regre.intercept_ + regre.coef_ * 3   # lo mismo
regre.score (x, y)                   # R2

Con otro módulo (quizá requiera «apt install python-statsmodels»):

import pandas as pd
mtcars = pd.read_csv ("http://bellman.ciencias.uniovi.es/~carleos/master/iot/ces/datos/mtcars.csv")
import statsmodels.formula.api as smf
regre = smf.ols ("mpg ~ wt", data = mtcars)
resul = regre.fit()
print (resul.summary ())
print (resul.summary2 ()) # incluye estimación de varianza residual (scale)

• conectar a carleos2.epv.uniovi.es
• equivalentemente, a 156.35.173.1
• la conexión SSH usa la puerta 22

• existen éstas:

  abel  daniel  guillermo  ivana  javiers  jorge  jorgel
  

• para cambiar la contraseña:

  passwd
  

• para evitar problemas por desconexión usar GNU screen

  screen -DUR
  

  ejecutar esa orden antes de empezar a trabajar y, en caso de desconexión, volver a ejecutarla para recuperar la sesión

• si quieres cambiar el idioma del sistema:
  • entorno en español

    export LC_ALL=es_ES.UTF-8
    

    • es = idioma castellano
    • ES = país España
    • UTF-8 = codificación de Unicode
  • entorno básico «internacional»

    export LC_ALL=C
    

    • mensajes en inglés
    • codificación ASCII

apropos compiler            # buscar sobre un término
man gcc                     # página de manual
info                        # ayuda de GNU

pwd                         # conocer directorio actual
ls                          # listar contenido
ls -l                       #   versión larga
ls -lrt                     #   deja al final los últimos modificados
ls -h                       #   tamaños legibles para humanos
mkdir ejemplo1              # crear directorio «ejemplo1»
cd ejemplo1                 # entrar en «ejemplo1»
cd ..                       # salir del directorio
rmdir                       # borrar un directorio vacío (seguro)
rm -rf                      # borrar un directorio y su contenido (atrevido)

• El directorio actual se representa por « . »
• El directorio superior se representa por « .. »
• El directorio personal se representa por « ~ »
• El directorio público está en « ~/public_html »

  Los archivos que el usuario mengano meta en su directorio público serán accesibles desde un navegador de hipertexto a través de la dirección

  http://carleos2.epv.uniovi.es/~mengano
  

  (salvo que, por ejemplo, se cree un fichero index.html en ese directorio).

touch fichero.txt           # crear fichero vacío
cp fichero.txt fichero1.txt # copiar
mv fichero.txt fichero0.txt # renombrar
mv fichero0.txt /tmp        # mover
rm fichero1.txt             # borrar fichero

nano fichero.c              # editor básico
jed fichero.c               # más completo, amigable
emacs fichero.c             # casi un sistema operativo

gcc fichero.c               # compilar
./a.out                     # ejecutar

gcc -o ej1.exe ejemplo1.c   # compilar
./ej1.exe                   # ejecutar

Si la biblioteca GSL se ha instalado con el sistema de paquetería oficial del sistema operativo (por ejemplo, APT en Debian y Ubuntu):

gcc ejemplo.c -lgsl -lgslcblas -lm

Si la biblioteca GSL se ha instalado desde código fuente (por ejemplo, descargándola desde la página de GNU para tener la última versión publicada):

LD_LIBRARY_PATH=/usr/local/lib gcc ejemplo.c -lgsl -lgslcblas -lm

Otra opción en este caso es usar compilación estática…

gcc -static ejemplo.c -lgsl -lgslcblas -lm

…pero compárense los tamaños de los ejecutables:

gcc -o ejemplo-dinamico.out ejemplo.c
gcc -o ejemplo-estatico.out -static ejemplo.c
ls -lh ejemplo-*.out

stdlib no añade mucho al tamaño; gsl, un poco más.

Según https://www.gnu.org/software/gsl/doc/html/statistics.html#examples

#include <stdio.h>
#include <gsl/gsl_sort.h>       /* para la mediana */
#include <gsl/gsl_statistics.h>

int main ()
{
  double datos[5] = {17.2, 18.1, 16.5, 18.3, 12.6};
  double media, mediana;

  printf ("Los datos son %g, %g, %g, %g, %g\n",
          datos[0], datos[1], datos[2], datos[3], datos[4]);
  media = gsl_stats_mean   (datos, 1, 5);
  gsl_sort (datos, 1, 5);
  printf ("Datos ordenados: %g, %g, %g, %g, %g\n",
          datos[0], datos[1], datos[2], datos[3], datos[4]);
  mediana = gsl_stats_median_from_sorted_data (datos, 1, 5);

  printf ("La media muestral es %g\n", media);
  printf ("La mediana es %g\n", mediana);
  return 0;
}

En caso de no querer modificar el vector original por el cálculo de la mediana, podría usarse memcpy de stdlib o bien usar vectores de GLS como en el siguiente listado:

#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_statistics.h>

int main ()
{
  gsl_vector * datos = gsl_vector_alloc (5);
  gsl_vector_set (datos, 0, 17.2);
  gsl_vector_set (datos, 1, 18.1);
  gsl_vector_set (datos, 2, 16.5);
  gsl_vector_set (datos, 3, 18.3);
  gsl_vector_set (datos, 4, 12.6);

  double media, mediana;

  printf ("Los datos son %g, %g, %g, %g, %g\n",
    gsl_vector_get (datos, 0),
    gsl_vector_get (datos, 1),
    gsl_vector_get (datos, 2),
    gsl_vector_get (datos, 3),
    gsl_vector_get (datos, 4));

  media = gsl_stats_mean
    (datos->data, datos->stride, datos->size);

  gsl_vector * copia = gsl_vector_alloc (5);
  gsl_vector_memcpy (copia, datos);
  gsl_sort_vector (copia);
  printf ("Los datos ordenados son %g, %g, %g, %g, %g\n",
    gsl_vector_get (copia, 0),
    gsl_vector_get (copia, 1),
    gsl_vector_get (copia, 2),
    gsl_vector_get (copia, 3),
    gsl_vector_get (copia, 4));

  mediana = gsl_stats_median_from_sorted_data
    (copia->data, copia->stride, copia->size);

  printf ("La media muestral es %g\n", media);
  printf ("La mediana es %g\n", mediana);

  printf ("Los datos originales siguen siendo %g, %g, %g, %g, %g\n",
    gsl_vector_get (datos, 0),
    gsl_vector_get (datos, 1),
    gsl_vector_get (datos, 2),
    gsl_vector_get (datos, 3),
    gsl_vector_get (datos, 4));

  return 0;
}

• nominales: nacionalidad, novela favorita…
• ordinales: calificación (malo, regular, bueno), nivel de estudios, rango…
• circulares: días de la semana, meses del año…

• distribución de frecuencias (apartado 2.4 del capítulo 1 de Estadística descriptiva y probabilidad)

  table (x)              # absolutas
  prop.table (table (x)) # relativas
  

• diagrama de sectores (para atributos nominales o circulares)

  pie (table (x))
  

• diagrama de barras

  barplot (table (x))
  

• discreta: número de hijos, número de averías…
• continua: diámetro del ánima, peso del proyectil…

• centralización

  mean (x)           # media
  median (x)         # mediana
  mean (x, trim=0.1) # media truncada
  

• dispersión

  sd (x)             # desviación típica
  IQR (x)            # recorrido intercuartílico
  

• posición

  quantile (x, 0.95) # percentil 95
  

• histograma

  hist (x)
  hist (x, col=2, breaks=100)
  stem (x)           # gráfico de tallo: a menudo más informativo
  

• afirmación sobre una distribución de probabilidad
• tipos
  • paramétrica
    • la media es 107
    • las varianzas son iguales en ambos grupos
  • no paramétrica
    • la distribución es gausiana
    • la mediana de las diferencias es negativa

• es un procedimiento para elegir entre dos hipótesis contrapuestas
• H0
  • hipótesis nula
  • siempre contiene la «igualdad»
    • la media es igual a 100
    • la proporción de defectos es menor o igual al 1%
    • las varianzas son iguales
    • la distribución es gausiana (densidad = campana de Gaus)
• H1
  • hipótesis alternativa
  • no suele contener la «igualdad» («H1 compuesta»)
    • la media es distinta de 100
    • la proporción de defectos es mayor que el 1%
    • las varianzas son distintas
    • la distribución no es gausiana (densidad ≠ campana de Gaus)

    hay excepciones («H1 simple»)

    • H0: proporción = 10%
    • H1: proporción = 15%
  • tipos según lateralidad:
    • H1: media ≠ 100 → contraste bilateral
    • H1: media < 100 → contraste unilateral izquierdo
    • H1: media > 100 → contraste unilateral derecho
• son equivalentes las siguientes formas de expresar contrastes unilaterales: \[ \left.\begin{array}{l} H_0:\quad \mu = 100 \\ H_1:\quad \mu < 100 \end{array}\right\} \equiv \left\{\begin{array}{l} H_0:\quad \mu \geq 100 \\ H_1:\quad \mu < 100 \end{array}\right. \] \[ \left.\begin{array}{l} H_0:\quad \mu = 100 \\ H_1:\quad \mu > 100 \end{array}\right\} \equiv \left\{\begin{array}{l} H_0:\quad \mu \leq 100 \\ H_1:\quad \mu > 100 \end{array}\right. \]

• errores

  realidad \ decisión	elijo H0	elijo H1
  H0 es cierta	👍	error 1
  H1 es cierta	error 2	👍

• presunción de inocencia
  • paradigma jurídico
    • H0: inocente
    • H1: culpable
  • paradigma de control de calidad
    • H0: proceso bajo control
    • H1: proceso fuera de control
  • nemotecnia: en caso de duda, hipótesis nula
• umbral de elección = nivel de significación = α
  • es la probabilidad del error de tipo 1 (probabilidad de condenar a un inocente)
  • el más habitual: 5%
  • depende mucho del campo de trabajo (control de calidad 3σ ⇒ α = 0,3%)

• cálculo del nivel crítico o P-valor
  • brevemente: es la probabilidad bajo H0 de obtener una muestra más rara que la mía;
  • con algo más de rigor: es la probabilidad, suponiendo que H0 es cierta, de obtener una muestra (del mismo tamaño que la que ya tengo) al menos tan rara como la mía;
  • por tanto, 0 ≤ P-valor ≤ 1
  • P-valor cercano a 0 cuando mi muestra sea muy rara bajo H0, luego no me creo H0;
  • P-valor grande cuando mi muestra no es rara bajo H0, luego puedo creerme H0;
• comparar P-valor y nivel de significación
  • P-valor ≥ α ⇒ elijo H0
  • P-valor < α ⇒ elijo H1
  • nemotecnia
    • P-valor alto, hipótesis alta (H0)
    • P-valor bajo, hipótesis de abajo (H1)

• (bajo H0) tamaño de un contraste
  • es la proporción de veces que se elige H1 siendo cierta H0;
  • estimable fácilmente mediante simulación:

    alfa <- 0.05; n <- 25; mean (replicate (1e5, shapiro.test (rnorm (n)) $ p.value < alfa))
    

  • debería parecerse a α
• bajo H1
  • potencia es la probabilidad de elegir H1 siendo cierta H1, es decir, 1 - Pr (error 2);
  • suele calcularse en función de un valor concreto dentro de H1

    ## p = proporción de defectuosos
    ## H0: p = 10%
    ## H1: p > 10%
    ## tamaño del contraste considerando muestras de tamaño 100:
    mean (replicate (1e5, binom.test (rbinom(1,100,0.1), 100, 0.1, alternative="greater") $ p.value < 0.05))
    ## potencia del contraste suponiendo que en realidad p = 20%
    mean (replicate (1e5, binom.test (rbinom(1,100,0.2), 100, 0.1, alternative="greater") $ p.value < 0.05))
    

  • puede calcularse la función de potencia, en función de los valores que toman los parámetros bajo H1:

    potencia1 <- function (p) 
      mean (replicate (1e5, binom.test (rbinom(1,100,p), 100, 0.1, alternative="greater") $ p.value < 0.05))
    potencia <- Vectorize (potencia1)
    plot (Pd <- seq (0.05, 0.95, 0.05), potencia (Pd), type="l")
    

• en poblaciones gausianas, un contraste de Wilcoxon es casi tan potente como un contraste t:

  > mean (replicate (1e5, t.test (rnorm(25), mu=0) $ p.value < 0.05))
  [1] 0.05036
  > mean (replicate (1e5, wilcox.test (rnorm(25), mu=0) $ p.value < 0.05))
  [1] 0.0496
  > mean (replicate (1e5, t.test (rnorm(25), mu=0.5) $ p.value < 0.05))
  [1] 0.66888
  > mean (replicate (1e5, wilcox.test (rnorm(25), mu=0.5) $ p.value < 0.05))
  [1] 0.64707
  

• en poblaciones no gausianas y con tamaño de muestra pequeño, el contraste t es conservador: su tamaño es menor que el nivel de significación;

  > mean (replicate (100000, t.test (rexp(15), rexp(15)) $ p.value < 0.05)) ; gc()
  [1] 0.04267
  > mean (replicate (100000, wilcox.test (rexp(15), rexp(15)) $ p.value < 0.05))
  [1] 0.0455
  

  el de Wilcoxon se acerca algo más al tamaño óptimo (5%)

• constraste de independencia
  • rachas
    • secuencia de valores iguales
    • calculables en R mediante «rle»
    • para variables continuas, tomar la mediana como referencia: «rle (m > median(m))»
  • número de rachas
    • implementación sencilla

      numrachas <- function (m,
                             FUN = function (m)
                                   length (rle (m > median(m)) $ lengths),
                             B = 9999)
      {
          Pder <- 1/(B+1) + mean (replicate (B, FUN (sample (m))) >= FUN (m))
          Pizq <- 1/(B+1) + mean (replicate (B, FUN (sample (m))) <= FUN (m))
          2 * min (Pder, Pizq)
      }
      ## ejemplo
      > numrachas (mtcars$mpg, B=99999)
      [1] 0.04646046
      

    • «runs.test» en «snpar» o «tseries»

      library (snpar)
      > runs.test(mtcars$mpg)
      
      Approximate runs rest
      
      data:  mtcars$mpg
      Runs = 11, p-value = 0.03215
      alternative hypothesis: two.sided
      
      > ?runs.test
      > runs.test(mtcars$mpg,exact=TRUE)
      
      Exact runs test
      
      data:  mtcars$mpg
      Runs = 11, p-value = 0.04741
      alternative hypothesis: two.sided
      

  • racha más larga

    rachamaslarga <- function (...)
    numrachas (FUN = function (m)
                     max (rle (m > median(m)) $ lengths),
               ...)
    ## ejemplo (es menos potente que el número de rachas)
    > rachamaslarga(mtcars$mpg,B=99999)
    [1] 0.09080091
    

• contraste de gausianidad
  • «shapiro.test»