Distribuciones de probabilidad
¿Por qué?
Distribuciones de probabilidad
Distribuciones de probabilidad
Continuas
Normal o Gaussiana
Beta
Gamma
Uniforme
Discretas
Poisson
Bernoulli
Binomial
Distribución Tipo
\[\begin{equation}
X \sim Nombre(parametro_1, parametro_2 )
\end{equation}\]
Dominio : \(X \in (-\infty,\infty)\) ; \(parametro_1 \in (-\infty,\infty)\) ; \(parametro_2 > 0\) (real)
En R/NIMBLE : dnombre(parametro1, parametro2)
Probability Distribution App
Distribución Normal o Gaussiana
Distribución Normal o Gaussiana
\[\begin{equation}
X \sim Normal(\mu,\sigma )
\end{equation}\]
Dominio : \(X \in (-\infty,\infty)\) ; \(\mu \in (-\infty,\infty)\) ; \(\sigma > 0\) (real)
En R/NIMBLE : dnorm(mean, sd)
Distribución Normal o Gaussiana
\[\begin{equation}
X \sim Normal(\mu,\tau )
\end{equation}\]
\[\begin{equation}
\tau = 1/\sigma
\end{equation}\]
Dominio : \(X \in (-\infty,\infty)\) ; \(\mu \in (-\infty,\infty)\) ; \(\sigma > 0\) (real)
En R/NIMBLE : dnorm(mean, sd)
Distribución Beta
Distribución Beta
\[\begin{equation}
X \sim Beta(\alpha,\beta )
\end{equation} \]
Dominio : \(X \in (0,1)\) ; \(\alpha > 0\) (real); \(\beta > 0\) (real)
En R/NIMBLE : dbeta(alpha, beta)
Distribución Beta
\[\begin{equation}
X \sim Beta(\mu,\sigma )
\end{equation}\]
\[\begin{equation}
\mu = \frac{\alpha}{\alpha + \beta} \text{ ; }\sigma = \sqrt{\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}}
\end{equation}\]
Dominio : \(X \in (0,1)\) ; \(\alpha > 0\) (real); \(\beta > 0\) (real)
En R/NIMBLE : dbeta(alpha, beta)
Distribución Gamma
Distribución Gamma
\[\begin{equation}
X \sim Gamma(\alpha,\beta)
\end{equation} \]
Dominio : \(X \in (0,\infty)\) ; \(\alpha > 0\) (real); \(\beta > 0\) (real)
En R/NIMBLE : dgamma(shape, rate)
Distribución Gamma
\[\begin{equation}
X \sim Gamma(\mu,\sigma)
\end{equation}\]
\[\begin{equation}
\mu = \frac{\alpha}{\beta} \text{ ; }\sigma = \sqrt(\frac{\alpha}{\beta^2})
\end{equation}\]
Dominio : \(X \in (0,\infty)\) ; \(\alpha > 0\) (real); \(\beta > 0\) (real)
En R/NIMBLE : dgamma(shape, rate)
Distribución de Poisson
Distribución de Poisson
\[\begin{equation}
X \sim Poisson(\lambda)
\end{equation}\]
Dominio : \(X \in (0,\infty)\) (natural); \(\lambda \in (0,\infty)\) (real)
En R/NIMBLE : dpois(lambda)
Distribución de Bernoulli
Distribución de Bernoulli
\[\begin{equation}
X \sim Bernoulli(p)
\end{equation}\]
Dominio : \(X = \{0,1\}\) ; \(p \in (0,1)\)
En R/NIMBLE : dbinom(1, prob) / dbern(prob)
Distribución Binomial
Distribución Binomial
\[\begin{equation}
X \sim Binomial(n, p)
\end{equation}\]
Dominio : \(X = \{0,1\}\) ;\(n \in (0,\infty)\) (natural) \(p \in (0,1)\)
En R/NIMBLE : dbinom(size, prob) / dbinom(prob, size)
Take home messages
Es necesario reflexionar sobre las variables que queremos modelizar
La “nomenclatura” puede ser una pesadilla: mejor entender lo que se está haciendo
Atentos a las diferentes parametrizaciones posibles
No somos estadísticos, pero Google es nuestro amigo!
Modelos Lineales Generalizados
\[\begin{equation}
Y \sim Normal(\mu,\sigma )
\end{equation}\]
\[\begin{equation}
Y \sim Normal(1000,300 )
\end{equation}\]
Modelo Lineal General
\[\begin{equation}
Y_i \sim Normal(\mu_i,\sigma )
\end{equation}\]
\[\begin{equation}
\mu_i = \beta_0 + \beta_1X_i
\end{equation}\]
Modelo Lineal General
\[\begin{equation}
Y_i \sim Normal(\mu_i,\sigma )
\end{equation}\]
\[\begin{equation}
\mu_i = \beta_0 + \beta_1X_i
\end{equation}\]
Modelo Lineal General
\[\begin{equation}
Y_i \sim Normal(\mu_i,\sigma )
\end{equation}\]
\[\begin{equation}
\mu_i = \beta_0 + \beta_1X_i
\end{equation}\]
Modelo Lineal General
\[\begin{equation}
Y_i \sim Normal(\mu_i,\sigma )
\end{equation}\]
\[\begin{equation}
\mu_i = \beta_0 + \beta_1X_i
\end{equation}\]
Modelo Lineal General
\[\begin{equation}
Y_i \sim Normal(\mu_i,\sigma )
\end{equation}\]
\[\begin{equation}
\mu_i = \beta_0 + \beta_1X_i
\end{equation}\]
Modelo Lineal Generalizado
¡Podemos usar cualquier otra distribución utilizando la misma maquinaria!
Se necesitan funciones vínculo (links) para unir covariables con parámetros: log , logit , etc.
Una vez entendidas las fórmulas es muy fácil programarlo en NIMBLE \[~\]
Gaussiana(\(\mu,\sigma\) )
Identidad \(\mu = \beta X\) \(\mu = \beta X\)
Poisson(\(\lambda\) )
Log \(log(\lambda) = \beta X\) \(\lambda = e^{(\beta X)}\)
Bernoulli(p )
Logit , probit, cloglog\(log(\frac{p}{1-p}) = \beta X\) \(p = \frac{e^{(\beta X)}}{1 + e^{(\beta X)}}\)
Binomial(n ,p )
Logit , probit, cloglog\(log(\frac{p}{1-p}) = \beta X\) \(p = \frac{e^{(\beta X)}}{1 + e^{(\beta X)}}\)
Hemos muestreado 100 cuadrículas contando excrementos de corzo y anotando la temperatura media: \[~\]
conteo temp
1 6 1.2
2 8 2.5
3 28 4.9
4 220 8.9
5 6 0.4
6 225 8.8
7 294 9.3
8 50 5.9
9 46 5.5
10 1 -1.3
11 2 0.5
12 4 0.1
13 59 6.2
14 7 2.6
15 92 7.2
16 13 4.0
17 86 6.6
18 386 9.9
19 17 2.6
20 104 7.3
21 274 9.2
22 2 0.5
23 53 5.8
24 1 -0.5
25 6 1.2
26 6 2.6
27 2 -1.8
28 11 2.6
29 183 8.4
30 7 2.1
31 18 3.8
32 37 5.2
33 16 3.9
34 5 0.2
35 139 7.9
36 52 6.0
37 121 7.5
38 2 -0.7
39 86 6.7
40 10 2.9
41 156 7.9
42 53 5.8
43 123 7.4
44 22 4.6
45 30 4.4
46 123 7.5
47 3 -1.7
48 17 3.7
49 78 6.8
50 58 6.3
51 21 3.7
52 182 8.3
53 11 3.3
54 2 0.9
55 1 -1.2
56 2 -0.8
57 4 1.8
58 19 4.2
59 42 5.9
60 14 2.9
61 250 9.0
62 10 1.5
63 14 3.5
64 5 2.0
65 46 5.8
66 5 1.1
67 14 3.7
68 109 7.2
69 0 -1.0
70 186 8.5
71 9 2.1
72 162 8.1
73 7 2.2
74 7 2.0
75 15 3.7
76 195 8.7
77 178 8.4
78 15 2.7
79 119 7.3
80 318 9.5
81 8 3.2
82 70 6.6
83 9 2.8
84 9 1.9
85 99 7.1
86 1 0.4
87 75 6.5
88 1 -0.5
89 1 0.9
90 1 -0.3
91 3 0.9
92 0 -1.3
93 43 5.7
94 177 8.5
95 85 7.3
96 117 7.6
97 13 3.5
98 8 2.9
99 151 7.7
100 38 5.3
\[~\] Podemos crear un modelo para relacionar el número de excrementos de corzo con la temperatura en cada cuadrícula. ¿Qué distribución podría utilizar?
\[\begin{equation}
N_i \sim Poisson(\lambda_i)
\end{equation}\]
\[\begin{equation}
N_i \sim Poisson(\lambda_i)
\end{equation}\]
\[\begin{equation}
\lambda_i = \beta_0 + \beta_1Temperatura_i
\end{equation}\]
conteo temp
1 6 1.2
2 8 2.5
3 28 4.9
4 220 8.9
5 6 0.4
6 225 8.8
7 294 9.3
8 50 5.9
9 46 5.5
10 1 -1.3
11 2 0.5
12 4 0.1
13 59 6.2
14 7 2.6
15 92 7.2
16 13 4.0
17 86 6.6
18 386 9.9
19 17 2.6
20 104 7.3
21 274 9.2
22 2 0.5
23 53 5.8
24 1 -0.5
25 6 1.2
26 6 2.6
27 2 -1.8
28 11 2.6
29 183 8.4
30 7 2.1
31 18 3.8
32 37 5.2
33 16 3.9
34 5 0.2
35 139 7.9
36 52 6.0
37 121 7.5
38 2 -0.7
39 86 6.7
40 10 2.9
41 156 7.9
42 53 5.8
43 123 7.4
44 22 4.6
45 30 4.4
46 123 7.5
47 3 -1.7
48 17 3.7
49 78 6.8
50 58 6.3
51 21 3.7
52 182 8.3
53 11 3.3
54 2 0.9
55 1 -1.2
56 2 -0.8
57 4 1.8
58 19 4.2
59 42 5.9
60 14 2.9
61 250 9.0
62 10 1.5
63 14 3.5
64 5 2.0
65 46 5.8
66 5 1.1
67 14 3.7
68 109 7.2
69 0 -1.0
70 186 8.5
71 9 2.1
72 162 8.1
73 7 2.2
74 7 2.0
75 15 3.7
76 195 8.7
77 178 8.4
78 15 2.7
79 119 7.3
80 318 9.5
81 8 3.2
82 70 6.6
83 9 2.8
84 9 1.9
85 99 7.1
86 1 0.4
87 75 6.5
88 1 -0.5
89 1 0.9
90 1 -0.3
91 3 0.9
92 0 -1.3
93 43 5.7
94 177 8.5
95 85 7.3
96 117 7.6
97 13 3.5
98 8 2.9
99 151 7.7
100 38 5.3
Distribución de Poisson
\[\begin{equation}
X \sim Poisson(\lambda)
\end{equation}\]
Dominio : \(X \in (0,\infty)\) (natural); \(\lambda \in (0,\infty)\) (real)
En R/NIMBLE : dpois(lambda)
\[\begin{equation}
N_i \sim Poisson(\lambda_i)
\end{equation}\]
\[\begin{equation}
log(\lambda_i) = \beta_0 + \beta_1Temperatura_i
\end{equation}\]
<- read.csv ("1_poisson_datos.csv" )<- glm (conteo ~ temp, family = poisson (link = "log" ))summary (m1)
Call:
glm(formula = conteo ~ temp, family = poisson(link = "log"))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.930115 0.055375 16.80 <2e-16 ***
temp 0.509376 0.007002 72.75 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 9052.170 on 99 degrees of freedom
Residual deviance: 79.327 on 98 degrees of freedom
AIC: 568.48
Number of Fisher Scoring iterations: 4
\[\begin{equation}
N_i \sim Poisson(\lambda_i)
\end{equation}\]
\[\begin{equation}
log(\lambda_i) = 0.93 + 0.5*Temperatura_i
\end{equation}\]
Take home messages
Un modelo es una simplificación de la realidad que nos ayuda a entenderla y/o predecirla
Es importante ser conscientes de las limitaciones de los modelos
Merece la pena empezar siempre por lo más sencillo e ir incrementando la complejidad
No somos estadísticos, pero Google es nuestro amigo!
Prácticas con R
R nos ofrece la posibilidad de simular números aleatorios obtenidos a partir de una distribución de probabilidad dada. Por ejemplo, podemos pedirle 20 números aleatorios obtenidos a partir de una distribución normal de media=0 y sd=1
rnorm (n= 20 , mean = 0 , sd = 1 )
[1] -1.28630053 -1.64060553 0.45018710 -0.01855983 -0.31806837 -0.92936215
[7] -1.48746031 -1.07519230 1.00002880 -0.62126669 -1.38442685 1.86929062
[13] 0.42510038 -0.23864710 1.05848305 0.88642265 -0.61924305 2.20610246
[19] -0.25502703 -1.42449465
Podemos graficar un histograma de frecuencias de estos números, para comprobar cómo se distribuyen realmente los números simulados
hist (rnorm (n= 20 , mean = 0 , sd = 1 ))
Cuanto mayor sea la muestra de números aleatorios, más se parecerá el histograma a su distribución real
hist (rnorm (n= 100000 , mean = 0 , sd = 1 ))
Cuanto mayor sea la muestra de números aleatorios, más se parecerá el histograma a su distribución real
hist (rnorm (n= 10 , mean = 0 , sd = 1 ))
Podemos generar números aleatorios sacados de una distribución cuyos parametros varían. Queremos 1000 números aleatorios de los cuales 500 provengan de una distribución normal media=-2 y sd=1 y otros 500 provengan de una distribución normal media=2 y sd=1
<- rep (c (- 2 ,2 ), each = 500 )hist (rnorm (n= 1000 , mean = mu, sd = 1 ), breaks = 20 )
Bucles
Son secuencias de instrucciones de código que se ejecuta repetidas veces, hasta que una condición deja de cumplirse. Uno de los bucles más utilizados es el denominado “for”. Podemos pedirle que nos imprima en pantalla las 5 primeras cifras de una serie de 100 números obtenidos a partir de una distribución normal de media=0 y sd=1
<- rnorm (n = 100 , mean = 0 , sd = 1 )for (i in 1 : 5 ){print (paste ("El número" ,i, "es el" , numeros[i]))
[1] "El número 1 es el -0.131220053875131"
[1] "El número 2 es el 0.955595821361478"
[1] "El número 3 es el 0.358458885309875"
[1] "El número 4 es el 1.16006162387297"
[1] "El número 5 es el -0.320703067252054"
Podemos consultar todas las distribuciones disponibles en R base en utilizando help("distribution"). Vamos a simular 1000 números obtenidos a partir de una distribución de Poisson con parámetro lambda=4
hist (rpois (n= 1000 , lambda = 4 ))#hist(rpois(n=1000, lambda = -4))
Vamos a simular 1000 números a partir de una distribución de Poisson cuyo parámetro \(\lambda\) va a variar en función a una variable aleatoria \(x_1 \sim Uniforme(-2,10)\) , con coeficientes \(b_0=1\) (intercepto) y \(b_1=0.5\) (pendiente).
<- round (runif (1000 ,- 2 ,10 ),1 )<- 1 <- 0.5 # log(lam) = b0 + b1*x1 <- exp (b0 + b1 * x1) # Inverso de la función vínculo log() <- rpois (1000 , lam)1 : 10 ]
[1] 38 296 78 2 16 30 3 14 1 96
Podemos ver la relación de nuestra variable predictora con los números aleatorios generados
Ahora podemos aplicar la maquinaria de los modelos lineales generalizados para ver si somos capaces de recuperar los parámetros simulados (\(b_0=1\) y \(b_1=0.5\) )
<- data.frame (conteos = conteos, x1 = x1)head (datos)
conteos x1
1 38 5.1
2 296 9.1
3 78 6.7
4 2 -0.8
5 16 3.0
6 30 5.0
<- glm (conteos ~ x1, family = poisson (link = "log" ), data = datos)summary (modelo)$ coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9839318 0.016960300 58.01382 0
x1 0.5022825 0.002043188 245.83276 0
Las simulaciones permiten recrear condiciones controladas (laboratorio). En palabras de Michael Schaub y Marc Kéry:
En las simulaciones la verdad es conocida , permite estudiar si el modelo produce estimas no sesgadas, investigar la violación de asunciones…
Útiles para hacer análisis de potencia (power analysis )
Se puede verificar la identificabilidad de nuestros parámetros
Es una prueba de que hemos entendido el modelo
<- datos[sample (1 : 1000 ,100 ),]<- datos[sample (1 : 1000 ,50 ),]<- datos[sample (1 : 1000 ,10 ),]<- datos[sample (1 : 1000 ,5 ),]<- glm (conteos ~ x1, family = poisson (link = "log" ), data = datos100)<- glm (conteos ~ x1, family = poisson (link = "log" ), data = datos50)<- glm (conteos ~ x1, family = poisson (link = "log" ), data = datos10)<- glm (conteos ~ x1, family = poisson (link = "log" ), data = datos5)#summary(modelo)$coefficients
0.9839318
0.0169603
0.5022825
0.0020432
1.0419032
0.0494432
0.4977652
0.0058699
0.9076790
0.0881629
0.5106495
0.0109850
0.9719651
0.1439569
0.5145671
0.0180198
0.7805831
0.2522782
0.5274006
0.0291127
Las simulaciones son ideales para entender ciertos fenómenos que siempre hemos escuchado en torno a la modelización:
¿Qué pasa si mis datos no siguen la distribución que estoy utilizando?
¿Qué ocurre si introduzco variables predictoras correlacionadas entre sí?
¿Qué efecto puede tener la heterogeneidad no modelada?
¿Cómo afectan a mis estimas los errores en mi base de datos?
