Distribución de Probabilidad
PHOBOS
2 de agosto de 2015
Funciones Basicas en R
Graficas y subgrupos
tab = read.csv("msleep_ggplot2.csv")
class(tab)## [1] "data.frame"dim(tab)## [1] 83 11c(tab$sleep_total, 1000)## [1] 12.1 17.0 14.4 14.9 4.0 14.4 8.7 7.0 10.1 3.0
## [11] 5.3 9.4 10.0 12.5 10.3 8.3 9.1 17.4 5.3 18.0
## [21] 3.9 19.7 2.9 3.1 10.1 10.9 14.9 12.5 9.8 1.9
## [31] 2.7 6.2 6.3 8.0 9.5 3.3 19.4 10.1 14.2 14.3
## [41] 12.8 12.5 19.9 14.6 11.0 7.7 14.5 8.4 3.8 9.7
## [51] 15.8 10.4 13.5 9.4 10.3 11.0 11.5 13.7 3.5 5.6
## [61] 11.1 18.1 5.4 13.0 8.7 9.6 8.4 11.3 10.6 16.6
## [71] 13.8 15.9 12.8 9.1 8.6 15.8 4.4 15.6 8.9 5.2
## [81] 6.3 12.5 9.8 1000.0plot(tab$brainwt, tab$sleep_total)
plot(tab$brainwt, tab$sleep_total, log="x")
tab[ c(1,2), ]## name genus vore order conservation sleep_total sleep_rem
## 1 Cheetah Acinonyx carni Carnivora lc 12.1 NA
## 2 Owl monkey Aotus omni Primates <NA> 17.0 1.8
## sleep_cycle awake brainwt bodywt
## 1 NA 11.9 NA 50.00
## 2 NA 7.0 0.0155 0.48tab[ tab$sleep_total > 18, ]## name genus vore order conservation
## 22 Big brown bat Eptesicus insecti Chiroptera lc
## 37 Thick-tailed opposum Lutreolina carni Didelphimorphia lc
## 43 Little brown bat Myotis insecti Chiroptera <NA>
## 62 Giant armadillo Priodontes insecti Cingulata en
## sleep_total sleep_rem sleep_cycle awake brainwt bodywt
## 22 19.7 3.9 0.1166667 4.3 0.00030 0.023
## 37 19.4 6.6 NA 4.6 NA 0.370
## 43 19.9 2.0 0.2000000 4.1 0.00025 0.010
## 62 18.1 6.1 NA 5.9 0.08100 60.000tab$sleep_total[ c(1,2) ]## [1] 12.1 17.0mean(tab$sleep_total[tab$sleep_total>18])## [1] 19.275Funciones Logicas, de Ordenamiento y de Coincidencias
which(tab$sleep_total > 18)## [1] 22 37 43 62tab$sleep_total[ which(tab$sleep_total > 18)[1] ]## [1] 19.7which(tab$sleep_total > 18 & tab$sleep_rem < 3)## [1] 43sort(tab$sleep_total)## [1] 1.9 2.7 2.9 3.0 3.1 3.3 3.5 3.8 3.9 4.0 4.4 5.2 5.3 5.3
## [15] 5.4 5.6 6.2 6.3 6.3 7.0 7.7 8.0 8.3 8.4 8.4 8.6 8.7 8.7
## [29] 8.9 9.1 9.1 9.4 9.4 9.5 9.6 9.7 9.8 9.8 10.0 10.1 10.1 10.1
## [43] 10.3 10.3 10.4 10.6 10.9 11.0 11.0 11.1 11.3 11.5 12.1 12.5 12.5 12.5
## [57] 12.5 12.8 12.8 13.0 13.5 13.7 13.8 14.2 14.3 14.4 14.4 14.5 14.6 14.9
## [71] 14.9 15.6 15.8 15.8 15.9 16.6 17.0 17.4 18.0 18.1 19.4 19.7 19.9order(tab$sleep_total)## [1] 30 31 23 10 24 36 59 49 21 5 77 80 11 19 63 60 32 33 81 8 46 34 16
## [24] 48 67 75 7 65 79 17 74 12 54 35 66 50 29 83 13 9 25 38 15 55 52 69
## [47] 26 45 56 61 68 57 1 14 28 42 82 41 73 64 53 58 71 39 40 3 6 47 44
## [70] 4 27 78 51 76 72 70 2 18 20 62 37 22 43tab$sleep_total[30]## [1] 1.9tab$sleep_total[ order(tab$sleep_total) ]## [1] 1.9 2.7 2.9 3.0 3.1 3.3 3.5 3.8 3.9 4.0 4.4 5.2 5.3 5.3
## [15] 5.4 5.6 6.2 6.3 6.3 7.0 7.7 8.0 8.3 8.4 8.4 8.6 8.7 8.7
## [29] 8.9 9.1 9.1 9.4 9.4 9.5 9.6 9.7 9.8 9.8 10.0 10.1 10.1 10.1
## [43] 10.3 10.3 10.4 10.6 10.9 11.0 11.0 11.1 11.3 11.5 12.1 12.5 12.5 12.5
## [57] 12.5 12.8 12.8 13.0 13.5 13.7 13.8 14.2 14.3 14.4 14.4 14.5 14.6 14.9
## [71] 14.9 15.6 15.8 15.8 15.9 16.6 17.0 17.4 18.0 18.1 19.4 19.7 19.9rank(tab$sleep_total)## [1] 53.0 77.0 66.5 70.5 10.0 66.5 27.5 20.0 41.0 4.0 13.5 32.5 39.0 55.5
## [15] 43.5 23.0 30.5 78.0 13.5 79.0 9.0 82.0 3.0 5.0 41.0 47.0 70.5 55.5
## [29] 37.5 1.0 2.0 17.0 18.5 22.0 34.0 6.0 81.0 41.0 64.0 65.0 58.5 55.5
## [43] 83.0 69.0 48.5 21.0 68.0 24.5 8.0 36.0 73.5 45.0 61.0 32.5 43.5 48.5
## [57] 52.0 62.0 7.0 16.0 50.0 80.0 15.0 60.0 27.5 35.0 24.5 51.0 46.0 76.0
## [71] 63.0 75.0 58.5 30.5 26.0 73.5 11.0 72.0 29.0 12.0 18.5 55.5 37.5match(c("Cow","Owl monkey","Cheetah"), tab$name)## [1] 5 2 1idx = match(c("Cotton rat","Owl monkey","Cheetah"), tab$name)
tab[idx,]## name genus vore order conservation sleep_total sleep_rem
## 68 Cotton rat Sigmodon herbi Rodentia <NA> 11.3 1.1
## 2 Owl monkey Aotus omni Primates <NA> 17.0 1.8
## 1 Cheetah Acinonyx carni Carnivora lc 12.1 NA
## sleep_cycle awake brainwt bodywt
## 68 0.15 12.7 0.00118 0.148
## 2 NA 7.0 0.01550 0.480
## 1 NA 11.9 NA 50.000Funciones aplicados a listas o vectores
vec = c("red","blue","red","green","green","yellow","orange")
fac = factor(vec)
fac## [1] red blue red green green yellow orange
## Levels: blue green orange red yellowlevels(fac)## [1] "blue" "green" "orange" "red" "yellow"vec == "blue"## [1] FALSE TRUE FALSE FALSE FALSE FALSE FALSEfac2 = factor(vec, levels=c("blue","green","yellow","orange","red"))
fac2## [1] red blue red green green yellow orange
## Levels: blue green yellow orange redlevels(fac2)## [1] "blue" "green" "yellow" "orange" "red"table(tab$order)##
## Afrosoricida Artiodactyla Carnivora Cetacea
## 1 6 12 3
## Chiroptera Cingulata Didelphimorphia Diprotodontia
## 2 2 2 2
## Erinaceomorpha Hyracoidea Lagomorpha Monotremata
## 2 3 1 1
## Perissodactyla Pilosa Primates Proboscidea
## 3 1 12 2
## Rodentia Scandentia Soricomorpha
## 22 1 5s = split(tab$sleep_total, tab$order)
s[["Rodentia"]]## [1] 14.4 7.0 9.4 12.5 8.3 14.9 14.2 14.3 12.8 12.5 14.6 7.7 14.5 11.5
## [15] 13.0 8.7 11.3 10.6 16.6 13.8 15.9 15.8mean(s[[17]])## [1] 12.46818# lapply(s, mean)
sapply(s, mean)## Afrosoricida Artiodactyla Carnivora Cetacea
## 15.600000 4.516667 10.116667 4.500000
## Chiroptera Cingulata Didelphimorphia Diprotodontia
## 19.800000 17.750000 18.700000 12.400000
## Erinaceomorpha Hyracoidea Lagomorpha Monotremata
## 10.200000 5.666667 8.400000 8.600000
## Perissodactyla Pilosa Primates Proboscidea
## 3.466667 14.400000 10.500000 3.600000
## Rodentia Scandentia Soricomorpha
## 12.468182 8.900000 11.100000tapply(tab$sleep_total, tab$order, mean)## Afrosoricida Artiodactyla Carnivora Cetacea
## 15.600000 4.516667 10.116667 4.500000
## Chiroptera Cingulata Didelphimorphia Diprotodontia
## 19.800000 17.750000 18.700000 12.400000
## Erinaceomorpha Hyracoidea Lagomorpha Monotremata
## 10.200000 5.666667 8.400000 8.600000
## Perissodactyla Pilosa Primates Proboscidea
## 3.466667 14.400000 10.500000 3.600000
## Rodentia Scandentia Soricomorpha
## 12.468182 8.900000 11.100000tapply(tab$sleep_total, tab$order, sd)## Afrosoricida Artiodactyla Carnivora Cetacea
## NA 2.5119050 3.5021638 1.5716234
## Chiroptera Cingulata Didelphimorphia Diprotodontia
## 0.1414214 0.4949747 0.9899495 1.8384776
## Erinaceomorpha Hyracoidea Lagomorpha Monotremata
## 0.1414214 0.5507571 NA NA
## Perissodactyla Pilosa Primates Proboscidea
## 0.8144528 NA 2.2098951 0.4242641
## Rodentia Scandentia Soricomorpha
## 2.8132994 NA 2.7046257Introducción a Variables Aleatorias
- Los valores obtenidos de la dieta menos los valores del control, para observar si el cambio es relativamente grande.
data<-read.csv("femaleMiceWeights.csv")
mean(data[13:24,2])-mean(data[1:12,2])## [1] 3.020833Obsvamos un cambio de 3 gramos, esto el 12.69%, la cual es un incremento significativo biologicamente. La pregunta ahora es si esta diferencia es obtenido aleatoriamente?, en otro grupo de ratones al azar, vamos a ver la diferencia de nuevo? o si esto es algo que se ve por casualidad.
Decripción de datos
- Graficamos los pesos de los ratones en el control y de la dieta con alta grasa.
s = split(data[,2], data[,1])
stripchart(s, vertical=TRUE, col=1:2)
abline(h=sapply(s, mean), col=1:2)
3 de los ratones de la dieta con mucha grasa son menores de la media del control y 3 de los ratones del control son mayores de la media con la dieta de mucha grasa.
Ejemplo de Muestras Aleatorias
- La función
samplegenerar muestras aleatorias de una población.
highfat = s[["hf"]]
sample(highfat, 6)## [1] 26.37 25.71 28.14 34.02 25.34 29.58sample(highfat, 6, replace=TRUE)## [1] 29.58 26.37 21.90 29.58 26.37 29.58- Para calcular la proporción, se pueden utilizar vectores lógicos
highfat > 30## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE TRUE
## [12] FALSEas.numeric(highfat > 30)## [1] 0 0 0 0 0 0 0 1 1 0 1 0sum(highfat > 30)## [1] 3mean(highfat > 30)## [1] 0.25Ejemplo Distribución Nula y p-valor
- La hipotesis nula, tiene la premisa que los cambios
(mean(dat[13:24,2]))que vimos al comparar los ratones controles y con la dieta son por cambios aleatorios. Teniendo toda la población
poblacion<-read.csv("femaleControlsPopulation.csv")
dim(poblacion)## [1] 225 1- Ahora si tomamos de la población 12 muestras aleatorias como control, calculamos la media. Observamos que los valores cambian porque cada vez seleccionamos muestras diferentes esto denota una variabilidad natural y aleatoriedad
control<-sample(poblacion[,1],12)
mean(control)## [1] 24.185control<-sample(poblacion[,1],12)
mean(control)## [1] 23.11333- Realizamos una Distribución nula, esto es, generamos controles y tratamienos de muestras aleatorias del origen, y realizamos una diferencia.
n<-1000
nula<-vector("numeric",n)
for(i in 1:n){
control<-sample(poblacion[,1],12)
tratamiento<-sample(poblacion[,1],12)
nula[i]<-mean(tratamiento)-mean(control)
}- Observamos la diferencia entre los datos originales y los datos creados aleatoriamente.
diffOrigen<-mean(data[13:24,2])-mean(data[1:12,2])
mean(nula>diffOrigen)## [1] 0.009Obteniendo la proporción de las veces que esto sucede que es 0.9% veces. Es el calculo del p-valor esto es la probabilidad de ver una diferencia tan grande como vimos en nuestro experimento, cuando la Ho es verdadera. La creación de la Ho es comparando una y otra ves los grupos (control y dieta), creando una distribución nula.
- Ilustración de la Distribucón Nula
n <- 100
plot(0,0,xlim=c(-5,5),ylim=c(1,30),type="n")
totals <- vector("numeric",11)
for(i in 1:n){
control <- sample(poblacion[,1],12)
treatment <- sample(poblacion[,1],12)
nulldiff <- mean(treatment) - mean(control)
j <- pmax(pmin(round(nulldiff)+6,11),1)
totals[j]<-totals[j]+1
text(j-6,totals[j],pch=15,round(nulldiff,1))
if(i<15)scan()
}
Visualialización de Distribución Nula
null = replicate(10000, mean(sample(poblacion[,1], 12)) - mean(sample(poblacion[,1], 12)))
hist(nula)
abline(v=diffOrigen, col="red")
abline(v=-diffOrigen, col="red")
Cual es la probabilidad de distribución nula para una cola mean(null > abs(diffOrigen)) = 0.0129, y para dos colas es mean(abs(null) > abs(diffOrigen)) = 0.0257
Distribución de Probabilidad
- Distribución es como una descripción compacta de muchos números. Para definir una distribución calculamos, para todos los valores posibles de \(\mathbf{a}\) la proporción de números en nuestra lista que se encuentra por debajo de \(\mathbf{a}\). Usamos la siguiente notación \[ F(a) \equiv \mbox{Pr}(x \leq a) \] Esto se llama la función de distribución acumulativa empírica.
Usando esto podemos decir cuántos individuos en cualquier intervalo \[ \mbox{Pr}(a \leq x \leq b) = F(b) - F(a) \] Si describimos la distribución, describimos la lista entera. - Función de Distribución Empirica
values <- seq(min(null),max(null),len=300)
myecdf <- ecdf(null)
plot(values,myecdf(values),type="l")
Ejemplo
- Este conjunto de datos contiene la esperanza de vida, el PIB per cápita y la población por países, cada cinco años, de 1952 a 2007.
Se trata de un extracto de un conjunto más grande y más completa de los datos disponibles sobre Gapminder.org
library(gapminder)
data(gapminder)
head(gapminder)## country continent year lifeExp pop gdpPercap
## 1 Afghanistan Asia 1952 28.801 8425333 779.4453
## 2 Afghanistan Asia 1957 30.332 9240934 820.8530
## 3 Afghanistan Asia 1962 31.997 10267083 853.1007
## 4 Afghanistan Asia 1967 34.020 11537966 836.1971
## 5 Afghanistan Asia 1972 36.088 13079460 739.9811
## 6 Afghanistan Asia 1977 38.438 14880372 786.1134data1952<-gapminder[gapminder$year == 1952,]La proporción de paises en 1952 con expectativa de vida menor o igual a 40 años e igual a 0.29. La proporción de paises en 1952 con expectativa de vida entre 40 y 60 años es igual a 0.4647887
Distribución Normal
Si en lugar de los números totales nos informan de las proporciones, entonces el histograma es una distribución de probabilidad. La distribución de probabilidad se aproxima a uno que es muy común en una naturaleza:la curva de campana o la distribución normal o distribución de Gauss. Cuando el histograma de una lista de números se aproxima a la distribución normal se puede utilizar una fórmula matemática conveniente para aproximar la proporción de individuos en cualquier intervalo dado: \[\mbox{Pr}(a < x < b) = \int_a^b \frac{1}{\sqrt{2\pi\sigma^2}} \exp{\left( \frac{-(x-\mu)^2}{2 \sigma^2} \right)} \, dx\] Donde \(\mu\) y \(\sigma\) son la media y la desviación estandar. Estos dos valores decribe la proporción de cualquier intervalo \[\mu = \frac{1}{M} \sum_{i=1}^{M} x_{i}\] \[\sigma^{2} = \frac{1}{M} \sum_{i=1}^{M} (x_{i} - \mu )^{2}\]
Unidades estandarizadas: Si sabemos que los datos tienen una aproximación normal, entonces puede convertirlo en lo que se llaman unidades estándar restando la media, y dividiéndolo por la desviación estándar para cada punto de la distribución. Vamos a tener una distribución normal cuando la
media = 0ySD = 1
\[Z_{k}=\frac{ x_{i}-\overline{X}}{ s_{x} }\]
Ejemplo de Distribución Normal
library(gapminder)
data(gapminder)
data1952<-gapminder[gapminder$year == 1952,]
par(mfrow=c(1,2))
hist(data1952$pop)
hist(log10(data1952$pop))
La media 6.6016108 y desviación estandar de los datos en log10 0.7070292
Distribución Estandarizada
pop52<-log10(data1952$pop)
qqnorm(pop52)
z<-(pop52-mean(pop52))/sd(pop52)
qqnorm(z);abline(0,1)
max(z)## [1] 3.03194