Examen Final Estadística Computacional

Luis Manuel Román García

Inferencia Gráfica

Primera parte

El precio de una casa se ve directamente afectado por la demanda, si la demanda es baja el precio tenderá a ser bajo. Por ende se esperaría que casas con climas desfavorables tengan menor demanda y en consecuencia menor precio. Como hipótesis nula supondremos una relación lineal entre ambas variables.

De tres individuos entrevistados, dos detectaron la gráfica con los valores reales.

> decrypt("ytIA Hiqi Un 0pTUqUpn ba")
[1] "True data in position 16"

Bajo la hipótesis nula, cada ensayo es independiente y se distribuye Bernoulli con . Con los tres ensayos tenemos una Binomial(3, .05). Luego entonces la probabilidad de observar nuestra muestra es de:

\theta = .05

\theta = .05

P(x = 2| H_0) =P(x = 2| n = 3, \theta = .05) = \frac{3!}{2!1!}.05^2.95^1=0.007125 = p-value

P(x = 2| H_0) =P(x = 2| n = 3, \theta = .05) = \frac{3!}{2!1!}.05^2.95^1=0.007125 = p-value

Inferencia Gráfica

Segunda Parte

Bootstrap Paramétrico

Bootstrap e Inferencia Bayesiana

x_1, x_2, ..., x_n.\quad iid\quad N(0,\sigma^2)

x_1, x_2, ..., x_n.\quad iid\quad N(0,\sigma^2)

Supongamos que tenemos

y queremos encontrar un estimador de máxima verosimilitud para . Para lograr esto calculamos la log-verosimilitud e igualamos a cero para obtener: que en este caso se reduce a:

\sigma^2

\sigma^2

log(f_x(x)) = log(\displaystyle\prod_{i = 1}^nf_{x_i}(x_i)) = -\frac{n}{2}log(2\pi\sigma^2)-\displaystyle\sum_{i = 1}^n\frac{(x-\mu)^2}{2\sigma^2}

log(f_x(x)) = log(\displaystyle\prod_{i = 1}^nf_{x_i}(x_i)) = -\frac{n}{2}log(2\pi\sigma^2)-\displaystyle\sum_{i = 1}^n\frac{(x-\mu)^2}{2\sigma^2}

\hat{\sigma^2}=\frac{1}{n}\displaystyle\sum_{i=1}^n(x-\mu)^2

\hat{\sigma^2}=\frac{1}{n}\displaystyle\sum_{i=1}^n(x-\mu)^2

\hat{\sigma^2} = \frac{1}{n}\displaystyle\sum_{i = 1}^nx^2

\hat{\sigma^2} = \frac{1}{n}\displaystyle\sum_{i = 1}^nx^2

Análisis Bayesiano

Bootstrap e Inferencia Bayesiana

Dado que no tenemos mucha información sobre la muestra ni el fenómeno que estas variables describen, lo más razonable es suponer una distribución a priori no muy informativa, o en otras palabras que no sea muy picuda.

\alpha = 10,\quad \beta = 100

\alpha = 10,\quad \beta = 100

Parece ser una opción razonable pues no acumula demasiada probabilidad en ningún intervalo.

Análisis Bayesiano

Bootstrap e Inferencia Bayesiana

Para calcular anlíticamente la distribución posterior, necesitamos de la distribución a priori y de la verosimilitud.

P(\theta|x)=\frac{p(x|\theta)p(\theta)}{p(x)}

P(\theta|x)=\frac{p(x|\theta)p(\theta)}{p(x)}

p(x|\theta)p(\theta) =

p(x|\theta)p(\theta) =

(\frac{1}{(2\pi\sigma^2)^{-\frac{n}{2}}}e^{\displaystyle\sum_{i=1}^n\frac{(x_i-\mu)^2}{2\sigma^2}})(\frac{\beta^{\alpha}}{\Gamma(\alpha)}{\sigma^2}^{(-1-\alpha)}e^{-\frac{\beta}{\sigma^2}})

(\frac{1}{(2\pi\sigma^2)^{-\frac{n}{2}}}e^{\displaystyle\sum_{i=1}^n\frac{(x_i-\mu)^2}{2\sigma^2}})(\frac{\beta^{\alpha}}{\Gamma(\alpha)}{\sigma^2}^{(-1-\alpha)}e^{-\frac{\beta}{\sigma^2}})

{\sigma^2}^{(-(\alpha + \frac{n}{2})-1)}e^{\displaystyle\sum_{i=1}^n\frac{(x_i-\mu)^2}{2\sigma^2}-\frac{\beta}{\sigma^2}}

{\sigma^2}^{(-(\alpha + \frac{n}{2})-1)}e^{\displaystyle\sum_{i=1}^n\frac{(x_i-\mu)^2}{2\sigma^2}-\frac{\beta}{\sigma^2}}

= \frac{\beta^{\alpha}{\sigma^2}^{(-1-\alpha)}}{\Gamma(\alpha)(2\pi\sigma^2)^{-\frac{n}{2}}}e^{\displaystyle\sum_{i=1}^n\frac{(x_i-\mu)^2}{2\sigma^2}-\frac{\beta}{\sigma^2}}

= \frac{\beta^{\alpha}{\sigma^2}^{(-1-\alpha)}}{\Gamma(\alpha)(2\pi\sigma^2)^{-\frac{n}{2}}}e^{\displaystyle\sum_{i=1}^n\frac{(x_i-\mu)^2}{2\sigma^2}-\frac{\beta}{\sigma^2}}

Análisis Bayesiano

Bootstrap e Inferencia Bayesiana

Supongamos ahora que el parámetro que deseamos estimar es . Como los estimadores de máxima verosimilitud son equivariantes bajo transformaciones 1 a 1 tenemos que:

\tau = log(\sigma)

\tau = log(\sigma)

\hat{\tau}_{MV} = log(\sqrt{\hat{\sigma}^2_{MV}})

\hat{\tau}_{MV} = log(\sqrt{\hat{\sigma}^2_{MV}})

Bootstrap Paramétrico

Bootstrap e Inferencia Bayesiana

Estimando con bootstrap paramétrico, obtenemos.

# Función máx verosimilitud
tau.hat <- function( x ){
	log( sqrt( sum( ( x - mean( x ) )^2 ) / length( x ) ) )
}

# Realización bootstrap del estimador de mv de 
# sigma^2
boot.var <- function( x ){
	x <- sample( x, length( x ), TRUE )
	tau.hat( x )
}

# N realizaciones bootstrap del estimador de 
# mv de sigma^2
n.boot.var <- function( x, n = 1000 ){
	replicate( n, boot.var( x ) )
}

set.seed( 123454321 )
data <- data.frame( x = n.boot.var( x ) )
ggplot(data = data,
	   aes(x = x ) ) +
geom_histogram( alpha = .5, fill = "#1A237E") + 
theme( text = element_text( color = "#E91E63",
							face  = "bold" ),
	   panel.background = element_blank( ) )

s.e <- sqrt(sum( ( data$x - 
				   mean( data$x ) )^2 ) /
	   nrow( data ) )

CI_{95\%} = [2.323498,2.541974]

CI_{95\%} = [2.323498,2.541974]

Análisis Bayesiano

Bootstrap e Inferencia Bayesiana

Estimando con una distribución a priori

# Bayesiana
set.seed( 123454321 )
log.sigma <- c()
for(i in 1:1000){
u   <- runif( 1000, .1, 300 ) 
log.sigma[i] <- log( sqrt( u )  )
}
ggplot(data = data.frame(x = log.sigma),
	   aes(x = x ) ) +
geom_histogram( alpha = .5, fill = "#1A237E") + 
theme( text = element_text( color = "#E91E63",
							face  = "bold" ),
	   panel.background = element_blank( ) )

U(a = .1, b = 300)

U(a = .1, b = 300)

Caminata aleatoria

MCMC

En esta primera sección, llevaremos a cabo 6,000 pasos de un modelo MCMC con brincos normales. Llevaremos a cabo el experimento con tres valores distintos de

N(0, 5),\quad N(0,.2),\quad N(0, 20)

N(0, 5),\quad N(0,.2),\quad N(0, 20)

\sigma^2

\sigma^2

Caminata aleatoria

MCMC

Para la tasa de rechazo simplemente ejecutamos el siguiente código:

> rechazo1 / pasos
[1] 0.5705
> rechazo2 / pasos
[1] 0.05533333
> rechazo3 / pasos
[1] 0.8751667
>

\%M_{N(0,.2)} = .05

\%M_{N(0,.2)} = .05

\% M_{N(0,5)} = .57

\% M_{N(0,5)} = .57

\% M_{N(0,20)} = .87

\% M_{N(0,20)} = .87

Caminata aleatoria

MCMC

Histogramas

M_{N(0,5)}

M_{N(0,5)}

M_{N(0,.2)}

M_{N(0,.2)}

M_{N(0,20)}

M_{N(0,20)}

# histogramas
ggplot(data = data.frame( x = 1:4001,
						   camino1 = camino1[ 2e3:6e3 ],
						   camino2 = camino2[ 2e3:6e3 ],
						   camino3 = camino3[ 2e3:6e3 ] ),
		aes(x = camino1 )) +
  geom_histogram(aes(y = ..density..), fill ="#1A237E", alpha = .7) +
   theme( text = element_text( color = "#E91E63",
							face  = "bold" ),
	   panel.background = element_blank( ) )
  
ggplot(data = data.frame( x = 1:4001,
						   camino1 = camino1[ 2e3:6e3 ],
						   camino2 = camino2[ 2e3:6e3 ],
						   camino3 = camino3[ 2e3:6e3 ] ),
		aes(x = camino2 )) +
  geom_histogram(aes(y = ..density..), fill ="#1565C0", alpha = .7) +
   theme( text = element_text( color = "#E91E63",
							face  = "bold" ),
	   panel.background = element_blank( ) )

ggplot(data = data.frame( x = 1:4001,
						   camino1 = camino1[ 2e3:6e3 ],
						   camino2 = camino2[ 2e3:6e3 ],
						   camino3 = camino3[ 2e3:6e3 ] ),
		aes(x = camino3 )) +
  geom_histogram(aes(y = ..density..), fill ="#D500F9", alpha = .7) +
   theme( text = element_text( color = "#E91E63",
							face  = "bold" ),
	   panel.background = element_blank( ) )

Caminata aleatoria

MCMC

Histogramas

Es interesante notar que el modelo que usa no alcanza a muestrear adecuadamente la distribución (esto se debe a que la distribución propuesta es muy estrecha). Mientras que el modelo con muestra ciertas irregularidades (esto se debe a la alta tasa de rechazo). Por otro lado, la elección de parece balancear bien ambos efectos pues nos da larepresentación más cercana a la verdadera distribución posterior.

\sigma = .02

\sigma = .02

\sigma = 20

\sigma = 20

\sigma = 5

\sigma = 5

Caminata aleatoria

MCMC

Predictiva Posterio

p(y)=\int p(y|\theta)p(\theta|x)d\theta

p(y)=\int p(y|\theta)p(\theta|x)d\theta

Recordando que podemos generar un histograma y un intervalo de 95% de confianza para la distribución de predicción. Cómo vimos, el modelo que mejor reproduce la distribución posterior es el que usa por lo que es el que usaremos.

\sigma = 5

\sigma = 5

#Distribución de predicción
camino1 <- camino1[ 2e3:6e3 ]
y_sims <- rnorm(1:length(camino1), mean(camino1), sd(camino1))

ggplot(data = data.frame(x = y_sims), aes(x = x)) + 
  geom_histogram(fill = "#D500F9") +
  geom_vline(aes(xintercept = mean(y_sims)), color = "#1565C0") +
   theme( text = element_text( color = "#E91E63",
							face  = "bold" ),
	   panel.background = element_blank( ) )
# banda de confianza.
   set.seed( 123454321 )
   mean.y <- c()
   sd.y   <- c()
   for(i in 1:1000){
   		samp.y      <- sample( y_sims, length( y_sims ), TRUE )
   		mean.y[ i ] <- mean( samp.y )
   		sd.y[ i ]   <- sd( samp.y )
	}
	samp.y <- sample(y_sims,1000, TRUE)
 ggplot(data = data.frame( x   = 1:1000, 
 						   y   = samp.y, 
 						   up  = mean.y + 2*sd.y,
 						   low = mean.y - 2*sd.y ),
 aes(x = x, y = y) ) + geom_point( col = "#E91E63" ) +
 geom_segment( aes(x = x, xend = x, y = up, yend = low),
 alpha = .3, col =  "#1565C0" ) + ylim(120, 140)+
   theme( text = element_text( color = "#E91E63",
							face  = "bold" ),
	   panel.background = element_blank( ) )

134.2846

126.4218

Estimación de Coeficientes

Convergencia

En esta sección buscamos estimar los valores de los coeficientes de un modelo de regresión lineal. En particular estamos suponiendo la siguientes distribuciones iniciales:

y_i \quad N(\alpha + \beta x_i, \sigma^2)

y_i \quad N(\alpha + \beta x_i, \sigma^2)

con

\beta\quad N(0, 1000)

\beta\quad N(0, 1000)

\sigma^2\quad U(0,1000)

\sigma^2\quad U(0,1000)

Caminata aleatoria

Modelos Jerárquicos

Si hacemos un simil entre los conejos del experimento y un lanzamiento de volados podemos ver que:

Un conejo desarrollará un tumor con probabilidad . Esto que puede ser visto como la probabilidad de éxito en un volado. Los conejos representan monedas.
Los experimentos representan volados. Mi incertidumbre sobre que impacto puede tener el tratamiento sobre los conejos se representa con el hiperparámetro .

\theta

\theta

\mu

\mu

Caminata aleatoria

Modelos Jerárquicos

Si hacemos un simil entre los conejos del experimento y un lanzamiento de volados podemos ver que:

Un conejo desarrollará un tumor con probabilidad . Esto que puede ser visto como la probabilidad de éxito en un volado. Los conejos representan monedas.
Los experimentos representan volados. Mi incertidumbre sobre que impacto puede tener el tratamiento sobre los conejos se representa con el hiperparámetro .

\theta

\theta

\mu

\mu

Caminata aleatoria

Modelos Jerárquicos

Estimaremos la distribución posterior de los parámetros con un muestreador de Gibbs.

\theta

\theta

\mu

\mu

K

Beta(a,b),\quad a = \mu k, \quad b = (1-\mu)k

Beta(a,b),\quad a = \mu k, \quad b = (1-\mu)k

Beta(A_\mu,B_\mu)

Beta(A_\mu,B_\mu)

Gamma(S_k, R_k)

Gamma(S_k, R_k)

#Modelo Jerarquico
load("data/rabbits.RData")
modelo_rabbit.txt <-
'
model{
  for(t in 1:N){
    x[t] ~ dbern(theta[coin[t]])
  }
  for(j in 1:nCoins){
    theta[j] ~ dbeta(a, b)
  }
  a <- mu * kappa
  b <- (1 - mu) * kappa
  # A_mu = 2, B_mu = 2
  mu ~ dbeta(1, 1)
  kappa ~ dgamma(1, 0.1)
}
'


cat(modelo_rabbit.txt, file = 'modelo_rabbit.bugs')

x <- rabbits$tumor
coin <- rabbits$experiment

jags.inits <- function(){
  list("mu" = runif(1, 0.1, 0.9),
    "kappa" = runif(1, 5, 20))
}

jags_fit <- jags(
  model.file = "modelo_rabbit.bugs",    # modelo de JAGS
  inits = jags.inits,   # valores iniciales
  data = list(x = x, coin = coin, nCoins = length(coin),  N = length(x)),    # lista con los datos
  parameters.to.save = c("mu", "kappa", "theta"),  # parámetros por guardar
  n.chains = 3,   # número de cadenas
  n.iter = 10000,    # número de pasos
  n.burnin = 1000   # calentamiento de la cadena
  )
sims_df <- data.frame(n_sim = 1:jags_fit$BUGSoutput$n.sims,
                      jags_fit$BUGSoutput$sims.matrix) %>% 
  dplyr::select(-deviance) %>%
  gather(parametro, value, -n_sim)

Caminata aleatoria

Modelos Jerárquicos

Gráficas

\mu

\mu

K

ggplot(filter(sims_df, parametro == "kappa"), aes(x = value)) +
  geom_histogram(alpha = 0.8, fill = "#1565C0") +
     theme( text = element_text( color = "#E91E63",
							face  = "bold" ),
	   panel.background = element_blank( ) )

ggplot(filter(sims_df, parametro == "mu"), aes(x = value)) +
  geom_histogram(alpha = 0.8, fill = "#1A237E") +
     theme( text = element_text( color = "#E91E63",
							face  = "bold" ),
	   panel.background = element_blank( ) )

adquiere la forma de una distribución normal mientras que K se inclina de manera importante por valores cercanos a cero.

\mu

\mu

Caminata aleatoria

Modelos Jerárquicos

Gráficas

	ggplot(filter(sims_df, ! parametro %in% c("mu", "kappa" ) ),
	aes(x = parametro, y = value ) ) + geom_boxplot( color = "#1A237E", 
	fill = "#1565C0") +
     theme( text = element_text( color = "#E91E63",
							face  = "bold" ),
	   panel.background = element_blank( ) )

Vemos que para prácticamente todas las estimaciones de se mantiene aproximadamente la misma distribución, sin embargo el número de outliers es importante.

\theta

\theta

Examen Final Estadística Computacional Luis Manuel Román García