Multi taller 2

---
title: "R Notebook"
output:
  pdf_document: default
  html_notebook: default
---

# Problemas prácticos de datos, distribuciones y correlación con respuestas parciales

### 1) Si
$$X \sim N_{2}\left(\begin{bmatrix}15\\20 \end{bmatrix},\begin{bmatrix}1&0.5\\0.5&1.25 \end{bmatrix}\right)$$ responda lo siguiente

### a) Encuentre $f(x)$ multiplicando las matrices necesarias (no uses dmvnorm()), tal que $X = \mu$

$$f(X|\mu,\Sigma)=\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}exp\left\{-\frac{1}{2}(x-\mu)'\Sigma^{-1}(x-\mu)\right\}$$

simplificando cuando $x=\mu$

$$f(X|\mu,\Sigma)=\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}$$

```{r}
library(expm)
library(mvtnorm)
library(MASS)

mu<-c(15,20)
sigma<-matrix(data = c(1, 0.5, 0.5,1.25),byrow = TRUE, ncol = 2)

#X <- mvrnorm(10000, mu, sigma)
dmvnorm(x=c(15,20), mean = mu, sigma = sigma)
#A%ˆ%20

#det(sigma)
f<- 1/(2*pi)*sqrt(det(sigma))
f
```

### Estimación raíz cuadrada de una matriz

```{r}
sigma%^%(1/2)

sigma%^%(1/2)

vect.val.sigma <- eigen(sigma)
vect.sigma <-vect.val.sigma$vectors

#vect.sigma<-matrix(vect.sigma)
val.sigma <- vect.val.sigma$values
val.sigma <- sqrt(diag(val.sigma))
(sigma.raiz <- vect.sigma %*% val.sigma%*%solve(vect.sigma))

```

### b) Halle $f(x)$ en $x =[14, 19]$ multiplicando realmente las matrices necesarias y utilizando dmvnorm().

```{r}
x<-c(14,19)
mu <- c(15, 20)
sigma <- matrix(data = c(1, 0.5, 0.5, 1.25), nrow = 2, ncol = 2, byrow = TRUE)
dmvnorm(x = x, mean = mu, sigma = sigma)
```

### c)  Represente gráficamente la distribución y estudie qué ocurre con ella cuando se realizan cambios en el vector de la media y en la matriz de covarianza. Asegúrese de examinar también los valores propios.

```{r}
mu<-c(15,20)
sigma<-matrix(data = c(1, 0.5, 0.5,1.25),byrow = TRUE, ncol = 2)

#X <- mvrnorm(10000, mu, sigma)
dmvnorm(x=c(15,20), mean = mu, sigma = sigma)
```

Al variar el vector de medias, el centro de las elipses se desplazan y al variar las covarianzas las elipses se
inclinan hacia la derecha o izquierda dependiendo del signo de la covarianza. Cuando la covarianza es cero
las elipses son circunfenrecias.

```{r}
# Mean vector and covariance matrix
mu1<-c(5,11)
mu2<-c(10,15)
mu3<-c(15,20)
mu4<-c(20,25)
mu5<-c(25,30)

sigma1<-matrix(data = c(1, 0, 0,1.25),byrow = TRUE, ncol = 2)
sigma2<-matrix(data = c(1, 0.25, 0.25,1.25),byrow = TRUE, ncol = 2)
2
sigma3<-matrix(data = c(1, 0.5, 0.5,1.25),byrow = TRUE, ncol = 2)
sigma4<-matrix(data = c(1, 0.75, 0.75,1.25),byrow = TRUE, ncol = 2)
sigma5<-matrix(data = c(1, 0.95, 0.95,1.25),byrow = TRUE, ncol = 2)

mu<-list(mu1,mu2,mu3,mu4,mu5)
sigma <-list(sigma1,sigma2,sigma3,sigma4,sigma5)

# To plot the ellipses
plot(x=NULL, y=NULL, xlab='Peso (kg)', ylab='Estatura (cm)',
xlim=c(2.5, 7.5), ylim=c(5, 15), las=1)
car::ellipse(center=mu1, shape=sigma1, col='tomato', fill=TRUE,
radius=sqrt(qchisq(p=0.90, df=2)))
car::ellipse(center=mu1, shape=sigma1, col='blue', fill=TRUE,
radius=sqrt(qchisq(p=0.50, df=2)))
legend('topleft', col=c('blue', 'tomato'), lwd=2, bty='n',
legend=c('50%', '90%'))

```

```{r}
plot(x=NULL, y=NULL, xlab='Peso (kg)', ylab='Estatura (cm)',
xlim=c(8, 12.5), ylim=c(10, 20), las=1)
car::ellipse(center=mu2, shape=sigma2, col='tomato', fill=TRUE,
radius=sqrt(qchisq(p=0.90, df=2)))
car::ellipse(center=mu2, shape=sigma2, col='blue', fill=TRUE,
radius=sqrt(qchisq(p=0.50, df=2)))
legend('topleft', col=c('blue', 'tomato'), lwd=2, bty='n',legend=c('50%', '90%'))

```

```{r}
plot(x=NULL, y=NULL, xlab='Peso (kg)', ylab='Estatura (cm)',
xlim=c(13, 18), ylim=c(15, 23), las=1)
car::ellipse(center=mu3, shape=sigma3, col='tomato', fill=TRUE,
radius=sqrt(qchisq(p=0.90, df=2)))
car::ellipse(center=mu3, shape=sigma3, col='blue', fill=TRUE,
radius=sqrt(qchisq(p=0.50, df=2)))
legend('topleft', col=c('blue', 'tomato'), lwd=2, bty='n',
legend=c('50%', '90%'))

```

```{r}
plot(x=NULL, y=NULL, xlab='Peso (kg)', ylab='Estatura (cm)',
xlim=c(18, 23), ylim=c(22, 28), las=1)
car::ellipse(center=mu4, shape=sigma4, col='tomato', fill=TRUE,
radius=sqrt(qchisq(p=0.90, df=2)))
car::ellipse(center=mu4, shape=sigma4, col='blue', fill=TRUE,
radius=sqrt(qchisq(p=0.50, df=2)))
legend('topleft', col=c('blue', 'tomato'), lwd=2, bty='n',
legend=c('50%', '90%'))

```

```{r}
plot(x=NULL, y=NULL, xlab='Peso (kg)', ylab='Estatura (cm)',
xlim=c(22.5, 28), ylim=c(27, 33), las=1)
car::ellipse(center=mu5, shape=sigma5, col='tomato', fill=TRUE,
radius=sqrt(qchisq(p=0.90, df=2)))
car::ellipse(center=mu5, shape=sigma5, col='blue', fill=TRUE,
radius=sqrt(qchisq(p=0.50, df=2)))
legend('topleft', col=c('blue', 'tomato'), lwd=2, bty='n',
legend=c('50%', '90%'))

```

```{r}
eigen(sigma1)
eigen(sigma2)
eigen(sigma3)
eigen(sigma4)
eigen(sigma5)

```

A medida que se las covarianzas pasa de cero a valores mayores de cero, el pirmer valor comienza a incrementar
y el segundo valor propio disminuye. Un análisis similar ocurre con los vectores propios, es decir, para el
primer vector propio su primera componente comienza a incrementarse y su segunda componente en cambio
disminuye de valor

### 2) Simule datos de una distribución normal multivariante de su elección (no es necesario que sea bivariante). Compare los estadísticos resumen multivariantes $\mu$ y $\Sigma$. Examine qué ocurre con los estadísticos resumen relativos $\mu$ y $\Sigma$a medida que N aumenta o disminuye.

```{r}
sigma_matrix=matrix(c(5,1,6,1,4,0,6,0,8), byrow = TRUE, nrow = 3)
mu = c(1, 2,4)

n=10
data2=mvrnorm(n,mu,sigma_matrix)
colMeans(data2)
cov(data2)
```

```{r}
n=30
data2=mvrnorm(n,mu,sigma_matrix)
colMeans(data2)
cov(data2)
```

```{r}
n=100
data2=mvrnorm(n,mu,sigma_matrix)
colMeans(data2)
cov(data2)
```

```{r}
n=500
data2=mvrnorm(n,mu,sigma_matrix)
colMeans(data2)
cov(data2)
```
A medida que N aumenta, los estadísticos van acercándose más al valor de los parámetros.

### 3) La matriz de covarianzas y correlaciones es la misma para variables aleatorias estandarizadas. ¿Por qué?

Si Z1 y Z2 son variables aleatorias estandarizadas, entonces ambas tienen una media de 0 y una varianza de
1. La correlación entre las dos variables es

$$\rho_{12}=\frac{Cov(Z_1,Z_2)}{\sqrt{Var(Z_1)Var(Z_2)}}=\frac{Cov(Z_1,Z_2)}{\sqrt{1\times1}}=Cov(Z_1,Z_2)$$

### 4)  Con respecto a los datos de goblet, complete lo siguiente:

### a) ¿Cuál es la unidad experimental?

Un goblet

### b)  Los investigadores del tema están interesados en agrupar copas que tienen la misma forma aunque puedan tener tamaños diferentes. Una forma sugerida por Manly y Alberto (2016) es ajustar los datos dividiendo cada medida por $X_3$ (altura) y luego completar análisis adicionales. A continuación se muestra cómo se realiza el ajuste de datos:

```{r}
#install.packages("readr")
suppressWarnings(library(readr))

#goblet <- read.csv("C:/Users/USUARIO/Downloads/goblet.csv")
goblet <- read.csv("D:/Descargas/DATABASES/goblet.csv")

head(goblet)

goblet2 <- data.frame(ID = goblet$goblet, w1 = goblet$x1/goblet$x3,
  w2 = goblet$x2/goblet$x3,
  w4 = goblet$x4/goblet$x3,
  w5 = goblet$x5/goblet$x3,
  w6 = goblet$x6/goblet$x3)

head(goblet2)

```

### c) Exponer los datos ajustados en una matriz.

```{r}
as.matrix(goblet2)

```



### d) Halla los estadísticos de resumen multivariantes para los datos ajustados.

```{r}
covar<-cov(goblet2[,2:6])
covar
```

### e)  Calcule la covarianza de las variables $X_1$ y $X_2$ ajustadas sin utilizar cov() o cualquier función que encuentre inmediatamente la covarianza. En otras palabras, programe en R la fórmula real que calcula la covarianza de dos variables y aplíquela a los datos de estas variables.

Lo que buscamos es:

$$Cov(X_1,X_2)=\frac{1}{n-1}\sum_{i=1}^n(x_{i1}-\bar x_1)(x_{i2}-\bar x_2)$$
```{r}
n=length(goblet$x1)
cov_12=(1/(n-1))*sum((goblet$x1-mean(goblet$x1))*(goblet$x2-mean(goblet$x2)))
cov_12
```