practica4

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

n <- 100
set.seed(123)
x1 <- round(runif(n, 0,10),2)
x2 <- round(runif(n, 0,10),2)

mediciones <- data.frame(x1=x1, x2=x2)
head(mediciones)

##     x1   x2
## 1 2.88 6.00
## 2 7.88 3.33
## 3 4.09 4.89
## 4 8.83 9.54
## 5 9.40 4.83
## 6 0.46 8.90

b0 <- 5
b1 <- 2
b2 <- 3

alpha <- 0.9
percentil_derecho<- qt(alpha/2, df = n-3)
percentil_izquierdo<- qt(1-(alpha/2), df = n-3)
run_and_fit <- function(rep) {
  epsilons <- rnorm(n)
  y <- b0 + b1*x1 + b2*x2 + epsilons
  replicacion <- data.frame(x1=x1, x2=x2, y=y)
  ajustado <- lm(y ~ x1+x2, data=replicacion, x=TRUE)
  error_standard <- summary(ajustado)$sigma * sqrt(solve(t(ajustado$x)%*%ajustado$x)[2,2])
  
  to_return <- as.list(coefficients(ajustado))
  to_return$error_standard <- error_standard
  # Calculo a mano el intervalo de confianza
  to_return$left<-  to_return$x1 - percentil_izquierdo * error_standard
  to_return$right <-  to_return$x1 - percentil_derecho * error_standard
  to_return
  
}

n_rep <- 1000
repeticiones<-purrr::map_dfr(1:n_rep, run_and_fit)
head(repeticiones)

## # A tibble: 6 × 6
##   `(Intercept)`    x1    x2 error_standard  left right
##           <dbl> <dbl> <dbl>          <dbl> <dbl> <dbl>
## 1          4.86  2.01  3.00         0.0346  2.01  2.01
## 2          5.02  2.06  2.96         0.0332  2.06  2.06
## 3          5.20  1.96  2.99         0.0370  1.96  1.97
## 4          5.15  2.03  2.96         0.0350  2.03  2.04
## 5          5.11  2.00  2.97         0.0335  1.99  2.00
## 6          5.01  2.00  2.97         0.0366  2.00  2.01

Por teorema, sabemos que el estimador B1 tiene distribucion normal con media B1 y varianza \(\sigma^2 (X^tX)^{-1}\).

Notese como los estimadores no son idependientes dado que tienen todos la misma varianza.

epsilons <- rnorm(n)
y <- b0 + b1*x1 + b2*x2 + epsilons
replicacion <- data.frame(x1=x1, x2=x2, y=y)
X2 <- lm(y ~ x1+x2, data=mediciones, x=TRUE)$x
X <- as.matrix(mediciones %>% select(x1,x2))

varianza_teorica <- solve(t(X) %*% X)[1,1] # sigma = 1
varianza_teorica_2 <- solve(t(X2) %*% X2)[2,2]# sigma = 1. Se agrega la columna de 1s
varianza_teorica_2

## [1] 0.001253859

list(media_muetral = mean(repeticiones$x1),
     media_teorica = b1,
     varianza_muestral = var(repeticiones$x1),
     varianza_teorica = varianza_teorica,
     varianza_teorica_2 = varianza_teorica)

## $media_muetral
## [1] 2.000783
## 
## $media_teorica
## [1] 2
## 
## $varianza_muestral
## [1] 0.001284684
## 
## $varianza_teorica
## [1] 0.0007061633
## 
## $varianza_teorica_2
## [1] 0.0007061633

solve(t(X) %*% X)

##               x1            x2
## x1  0.0007061633 -0.0005294574
## x2 -0.0005294574  0.0006970093

sum((repeticiones$x1 - mean(repeticiones$x1))**2) / (n_rep-1)

## [1] 0.001284684

repeticiones %>% ggplot() + geom_histogram(aes(x=x1))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

La distribucion del estadistico es T student con n-p grados de libertad, en este caso 3.

t = (repeticiones$x1-b1) / sqrt(var(repeticiones$x1))
plot(density(t))
lines(density(rt(10000, n-3)), col="red")
lines(density(rnorm(10000)), col="blue")

# Intervalos de confianza: veamos el cubrimiento empirico

repeticiones$cubre <- repeticiones$x1 > repeticiones$left &  repeticiones$x1 < repeticiones$right
sum(repeticiones$cubre)

## [1] 1000