Análisis Estadístico con R - Parte I

Distribución Binomial

La distribución binomial es una distribución de probabilidad discreta que cuenta el número de éxitos en una secuencia de n ensayos de Bernoulli independientes entre sí, con una probabilidad fija p de ocurrencia del éxito entre los ensayos. Para \(n = 1\), la binomial se convierte, de hecho, en una distribución de Bernoulli.

\[p(x)={n \choose x}p^x(1-p)^{n-x}\quad x = \{0, 1, 2, \dots , n\}\]

# Semilla aleatoria
set.seed(999)

# Muestras binomiales de tamaño 10,50,100,1000
bin1<-list("N=10"=rbinom(10,10,0.1),"N=50"=rbinom(50,10,0.1),
        "N=100"=rbinom(100,10,0.1),"N=1000"=rbinom(1000,10,0.1))

# Clases del histogramas
b<-(-0.5:10.5)

# Múltiples gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(bin1[i]),freq=F,breaks=b,ylim=c(0,0.5),col='lightblue',xlab="X",
       ylab="Densidad",main=paste0("Binomial(10,0.1)\n",names(bin1)[i]))
  points(b+0.5,dbinom(b+0.5,10,0.1),col=2,pch=19)
  points(b+0.5,dbinom(b+0.5,10,0.1),col=2,type='h',lwd=2)
}

# Semilla aleatoria
set.seed(999)
# Muestras binomiales de tamaño 10,50,100,1000
bin1<-list("N=10"=rbinom(10,10,0.5),"N=50"=rbinom(50,10,0.5),
        "N=100"=rbinom(100,10,0.5),"N=1000"=rbinom(1000,10,0.5))
# Clases del histogramas
b<-(-0.5:10.5)
# Múltiples gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(bin1[i]),freq=F,breaks=b,ylim=c(0,0.3),col='lightblue',xlab="X",
       ylab="Densidad",main=paste0("Binomial(10,0.5)\n",names(bin1)[i]))
  points(b+0.5,dbinom(b+0.5,10,0.5),col=2,pch=19)
  points(b+0.5,dbinom(b+0.5,10,0.5),col=2,type='h',lwd=2)
}

# Semilla aleatoria
set.seed(999)
# Muestras binomiales de tamaño 10,50,100,1000
bin1<-list("N=10"=rbinom(10,10,0.9),"N=50"=rbinom(50,10,0.9),
        "N=100"=rbinom(100,10,0.9),"N=1000"=rbinom(1000,10,0.9))
# Clases del histogramas
b<-(-0.5:10.5)
# Múltiples gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(bin1[i]),freq=F,breaks=b,ylim=c(0,0.5),col='lightblue',xlab="X",
       ylab="Densidad",main=paste0("Binomial(10,0.9)\n",names(bin1)[i]))
  points(b+0.5,dbinom(b+0.5,10,0.9),col=2,pch=19)
  points(b+0.5,dbinom(b+0.5,10,0.9),col=2,type='h',lwd=2)
}

Distribución Geométrica

La variable aleatoria geométrica representa el número de fallas que ocurren antes de que se presente el primer éxito. La función de masa de probabilidad está dada por:

\[p(x) = (1 - p)^{x-1}p\]

# Semilla Aleatoria
set.seed(999)

# Muestras Geométricas de tamaño 10, 50, 100 y 1000
geo<-list("N=10"=rgeom(10,0.3),"N=50"=rgeom(50,0.3),
          "N=100"=rgeom(100,0.3),"N=1000"=rgeom(1000,0.3))

# Clases del histograma
b<-seq(-0.5,max(geo[[4]])+0.5)

# Múltiples Gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(geo[i]),freq=F,breaks=b,ylim=c(0,0.3),col='lightgreen',xlab="X",
       ylab="Densidad",main=paste0("Geométrica(0.3)\n",names(geo)[i]))
  points(b+0.5,dgeom(b+0.5,0.3),col=2,pch=19,cex=0.8)
  points(b+0.5,dgeom(b+0.5,0.3),col=2,type='h',lwd=2)
}

# Semilla Aleatoria
set.seed(999)

# Muestras Geométricas de tamaño 10, 50, 100 y 1000
geo<-list("N=10"=rgeom(10,0.5),"N=50"=rgeom(50,0.5),
        "N=100"=rgeom(100,0.5),"N=1000"=rgeom(1000,0.5))

# Clases del histograma
b<-seq(-0.5,max(geo[[4]])+0.5)

# Múltiples Gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(geo[i]),freq=F,breaks=b,ylim=c(0,0.6),col='lightgreen',xlab="X",
       ylab="Densidad",main=paste0("Geométrica(0.5)\n",names(geo)[i]))
  points(b+0.5,dgeom(b+0.5,0.5),col=2,pch=19)
  points(b+0.5,dgeom(b+0.5,0.5),col=2,type='h',lwd=2)
}

# Semilla Aleatoria
set.seed(999)

# Muestras Geométricas de tamaño 10, 50, 100 y 1000
geo<-list("N=10"=rgeom(10,0.7),"N=50"=rgeom(50,0.7),
        "N=100"=rgeom(100,0.7),"N=1000"=rgeom(1000,0.7))

# Clases del histograma
b<-seq(-0.5,max(geo[[4]])+0.5)

# Múltiples Gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(geo[i]),freq=F,breaks=b,ylim=c(0,0.8),col='lightgreen',xlab="X",
       ylab="Densidad",main=paste0("Geométrica(0.7)\n",names(geo)[i]))
  points(b+0.5,dgeom(b+0.5,0.7),col=2,pch=19)
  points(b+0.5,dgeom(b+0.5,0.7),col=2,type='h',lwd=2)
}

Distribución Binomial Negativa

Esta distribución puede considerarse como una extensión o ampliación de la distribución geométrica. La distribución binomial negativa es un modelo adecuado para tratar aquellos procesos en los que se repite un determinado ensayo o prueba hasta conseguir un número determinado de resultados favorables (por vez primera) .Es por tanto de gran utilidad para aquellos muestreos que procedan de esta manera. Si el número de resultados favorables buscados fuera 1 estaríamos en el caso de la distribución geométrica.

\[p(x)= {x-1 \choose x-k}p^k(1-p)^{x-k} = {x-1 \choose k-1}p^k(1-p)^{x-k}\]

# Semilla Aleatoria
set.seed(999)

# Muestras Geométricas de tamaño 10, 50, 100 y 1000
nbin<-list("N=10"=rnbinom(10,3,0.3),"N=50"=rnbinom(50,3,0.3),
        "N=100"=rnbinom(100,3,0.3),"N=1000"=rnbinom(1000,3,0.3))

# Clases del histograma
b<-seq(-0.5,max(nbin[[4]])+0.5)

# Múltiples Gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(nbin[i]),freq=F,breaks=b,ylim=c(0,0.15),col='moccasin',xlab="X",
       ylab="Densidad",main=paste0("Nbinomial(3,0.3)\n",names(nbin)[i]))
  points(b+0.5,dnbinom(b+0.5,3,0.3),col=2,pch=19,cex=0.8)
  points(b+0.5,dnbinom(b+0.5,3,0.3),col=2,type='h',lwd=1)
}

# Semilla Aleatoria
set.seed(999)
# Muestras Geométricas de tamaño 10, 50, 100 y 1000
nbin<-list("N=10"=rnbinom(10,3,0.5),"N=50"=rnbinom(50,3,0.5),
        "N=100"=rnbinom(100,3,0.5),"N=1000"=rnbinom(1000,3,0.5))
# Clases del histograma
b<-seq(-0.5,max(nbin[[4]])+0.5)
# Múltiples Gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(nbin[i]),freq=F,breaks=b,ylim=c(0,0.3),col='moccasin',xlab="X",
       ylab="Densidad",main=paste0("Nbinomial(3,0.5)\n",names(nbin)[i]))
  points(b+0.5,dnbinom(b+0.5,3,0.5),col=2,pch=19)
  points(b+0.5,dnbinom(b+0.5,3,0.5),col=2,type='h',lwd=2)
}

# Semilla Aleatoria
set.seed(999)
# Muestras Geométricas de tamaño 10, 50, 100 y 1000
nbin<-list("N=10"=rnbinom(10,3,0.7),"N=50"=rnbinom(50,3,0.7),
        "N=100"=rnbinom(100,3,0.7),"N=1000"=rnbinom(1000,3,0.7))
# Clases del histograma
b<-seq(-0.5,max(nbin[[4]])+0.5)
# Múltiples Gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(nbin[i]),freq=F,breaks=b,ylim=c(0,0.5),col='moccasin',xlab="X",
       ylab="Densidad",main=paste0("Nbinomial(3,0.7)\n",names(nbin)[i]))
  points(b+0.5,dnbinom(b+0.5,3,0.7),col=2,pch=19)
  points(b+0.5,dnbinom(b+0.5,3,0.7),col=2,type='h',lwd=2)
}

# Semilla Aleatoria
set.seed(999)

# Muestras Geométricas de tamaño 1000
nbin<-list("N=1000"=rnbinom(1000,30,0.3),"N=1000"=rnbinom(1000,30,0.5),
        "N=1000"=rnbinom(1000,30,0.7))

# Clases del histograma
b<-list(seq(-0.5,max(nbin[[1]])+0.5),seq(-0.5,max(nbin[[2]])+0.5),
     seq(-0.5,max(nbin[[3]])+0.5))
p<-c(0.3,0.5,0.7)

# Múltiples Gráficos
par(mfrow=c(1,3))
for(i in 1:3)
{
  hist(unlist(nbin[i]),freq=F,breaks=b[[i]],col='moccasin',xlab="X",
       ylab="Densidad",main=paste0("Nbinomial(30,",p[i],")\n",names(nbin)[i]))
  points(b[[i]]+0.5,dnbinom(b[[i]]+0.5,30,p[i]),col=2,pch=19,cex=0.7)
  points(b[[i]]+0.5,dnbinom(b[[i]]+0.5,30,p[i]),col=2,type='h',lwd=1)
}

Distribución Hipergeométrica

La distribución hipergeométrica es una distribución discreta que modela el número de eventos en una muestra de tamaño fijo cuando se conoce el número total de elementos en la población de la cual proviene la muestra. Las muestras no tienen reemplazo, por lo que cada elemento de la muestra es diferente. Cuando se elige un elemento de la población, no se puede volver a elegir. Por lo tanto, la probabilidad de que un elemento en particular sea seleccionado aumenta con cada ensayo, suponiendo que aún no ha sido seleccionado.

\[P(X=x)=\frac{{k \choose x}{{N-k \choose n-x}}}{{N \choose n}} \quad x=0,1,...,k\]

Donde \(N\) es el tamaño de población, \(n\) es el tamaño de la muestra extraída, \(k\) es el número de elementos en la población original que pertenecen a la categoría deseada y \(x\) es el número de elementos en la muestra que pertenecen a dicha categoría.

# Semilla aleatoria
set.seed(999)

# Muestras Hipergeométricas de tamaño 30,100,300,1000
N<-80; k<-10; n<-10
hyp<-list("Tamaño=30"=rhyper(30,n,N-n,k),"Tamaño=100"=rhyper(100,n,N-n,k),
          "Tamaño=300"=rhyper(300,n,N-n,k),"Tamaño=1000"=rhyper(1000,n,N-n,k))

# Clases del histogramas
b<-seq(-0.5,max(hyp[[4]])+0.5)

# Múltiples gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
   hist(unlist(hyp[i]),freq=F,breaks=b,ylim=c(0,0.45),col='lightblue',xlab="X",
        main=paste0("Hipergeométrica k = ",k,' N = ',N,' n = ',n,'\n',names(hyp)[i]),
        cex.main=0.8,ylab="Densidad")
   points(b+0.5,dhyper(b+0.5,n,N-n,k),col=2,pch=19,cex=0.8)
   points(b+0.5,dhyper(b+0.5,n,N-n,k),col=2,type='h',lwd=2)
}

# Semilla aleatoria
set.seed(999)
# Muestras Hipergeométricas de tamaño 30,100,300,1000
N<-60; k<-10; n<-25
hyp<-list("Tamaño=30"=rhyper(30,n,N-n,k),"Tamaño=100"=rhyper(100,n,N-n,k),
          "Tamaño=300"=rhyper(300,n,N-n,k),"Tamaño=1000"=rhyper(1000,n,N-n,k))
# Clases del histogramas
b<-seq(-0.5,max(hyp[[4]])+0.5)
# Múltiples gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
   hist(unlist(hyp[i]),freq=F,breaks=b,ylim=c(0,0.35),col='lightblue',xlab="X",
        main=paste0("Hipergeométrica k = ",k,' N = ',N,' n = ',n,'\n',names(hyp)[i]),
        cex.main=.8,ylab="Densidad")
   points(b+0.5,dhyper(b+0.5,n,N-n,k),col=2,pch=19,cex=0.8)
   points(b+0.5,dhyper(b+0.5,n,N-n,k),col=2,type='h',lwd=2)
}

# Semilla aleatoria
set.seed(999)
# Muestras Hipergeométricas de tamaño 30,100,300,1000
N<-60; k<-35; n<-25
hyp<-list("Tamaño=30"=rhyper(30,n,N-n,k),"Tamaño=100"=rhyper(100,n,N-n,k),
          "Tamaño=300"=rhyper(300,n,N-n,k),"Tamaño=1000"=rhyper(1000,n,N-n,k))
# Clases del histogramas
b<-seq(-0.5,max(hyp[[4]])+0.5)
# Múltiples gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
   hist(unlist(hyp[i]),freq=F,breaks=b,ylim=c(0,0.25),col='lightblue',xlab="X",
        main=paste0("Hipergeométrica k = ",k,' N = ',N,' n = ',n,'\n',names(hyp)[i]),
        cex.main=.8,ylab="Densidad")
   points(b+0.5,dhyper(b+0.5,n,N-n,k),col=2,pch=19,cex=0.5)
   points(b+0.5,dhyper(b+0.5,n,N-n,k),col=2,type='h',lwd=1)
}

Distribución Poisson

La distribución de Poisson es una distribución de probabilidad discreta que expresa, a partir de una frecuencia de ocurrencia media, la probabilidad de que ocurra un determinado número de eventos durante cierto período de tiempo. Concretamente, se especializa en la probabilidad de ocurrencia de sucesos con probabilidades muy pequeñas, o sucesos “raros”.

\[p(x;\lambda)=\frac{e^{-\lambda} \lambda^x}{x!}\]

# Semilla aleatoria
set.seed(999)

# Muestras poisson de tamaño 30,100,300,1000
pois<-list("N=30"=rpois(30,2),"N=100"=rpois(100,2),
        "N=300"=rpois(300,2),"N=1000"=rpois(1000,2))

# Clases del histogramas
b<-seq(-0.5,max(pois[[4]])+0.5)

# Múltiples gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(pois[i]),freq=F,breaks=b,ylim=c(0,0.3),col='lemonchiffon',xlab="X",
       ylab="Densidad",main=paste0("Poisson (2)\n",names(pois)[i]))
  points(b+0.5,dpois(b+0.5,2),col=2,pch=19)
  points(b+0.5,dpois(b+0.5,2),col=2,type='h',lwd=2)
}

# Semilla aleatoria
set.seed(999)
# Muestras poisson de tamaño 30,100,300,1000
pois<-list("N=30"=rpois(30,5),"N=100"=rpois(100,5),
        "N=300"=rpois(300,5),"N=1000"=rpois(1000,5))
# Clases del histogramas
b<-seq(-0.5,max(pois[[4]])+0.5)
# Múltiples gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(pois[i]),freq=F,breaks=b,ylim=c(0,0.20),col='lemonchiffon',xlab="X",
       ylab="Densidad",main=paste0("Poisson (5)\n",names(pois)[i]))
  points(b+0.5,dpois(b+0.5,5),col=2,pch=19)
  points(b+0.5,dpois(b+0.5,5),col=2,type='h',lwd=2)
}

# Semilla aleatoria
set.seed(999)
# Muestras poisson de tamaño 30,100,300,1000
pois<-list("N=30"=rpois(30,10),"N=100"=rpois(100,10),
        "N=300"=rpois(300,10),"N=1000"=rpois(1000,10))
# Clases del histogramas
b<-seq(-0.5,max(pois[[4]])+0.5)
# Múltiples gráficos
par(mfrow=c(2,2))
for(i in 1:4)
{
  hist(unlist(pois[i]),freq=F,breaks=b,ylim=c(0,0.15),col='lemonchiffon',xlab="X",
       ylab="Densidad",main=paste0("Poisson (10)\n",names(pois)[i]))
  points(b+0.5,dpois(b+0.5,10),col=2,pch=19,cex=0.5)
  points(b+0.5,dpois(b+0.5,10),col=2,type='h',lwd=1)
}

Distribución Multinomial

require(mvtnorm)
require(plot3D)
require(scatterplot3d)

La distribución multinomial es una generalización de la distribución binomial, la cual es la probabilidad de un número de éxitos en \(N\) sucesos de Bernoulli independientes, con la misma probabilidad de éxito en cada suceso. En una distribución multinomial, el análogo a la distribución de Bernoulli es la distribución categórica, donde cada suceso concluye en únicamente un resultado de un número finito \(K\) de los posibles, con probabilidades \(p_1, \dots , p_k\) (tal que \(p_i \geq 0\) para \(i\) entre \(1\) y \(K\) y \(\sum_{i=1}^k p_i = 1\)) y con \(n\) sucesos independientes.

Entonces sea la variable aleatoria \(x_i\), que indica el número de veces que se ha dado el resultado \(i\) sobre los n sucesos. El vector \(x=(x_1,..., x_k)\) sigue una distribución multinomial con parámetros \(n\) y \(p\), donde \(p=(p_1,...,p_k)\) la distribución multinomial es una generalización de la distribución binomial.

Para el caso de la distribución trinomial, la función de masa de probabilidad es:

\[ p(x,y) = { \displaystyle {n! \over x!y!(n-x-y)!}p_1^{x}p_2^{y}(1-p_1-p_2)^{n-x-y}}, \quad \mbox{cuando } x+y \leq n \]

# Semilla aleatoria
set.seed(999)

# Vector de probablidades
p<-c(0.1,0.8);p[3]<-1-sum(p)

# Número de ensayos
n<-10

# Muestra aleatoria multinomial de tamaño 1000
x<-t(rmultinom(10000,n,p))

# Matriz de densidades
d<-table(x[,1],x[,2])/100

# Histograma 3D
par(mar=c(2,2,2,2))
hist3D(z=d,col='aquamarine',border=1,phi = 30,theta=-25,d=2,bty='g',
       shade = .35, zlab = "Densidad",space = 0.1,cex.axis = 1e-9)

# Semilla aleatoria
set.seed(999)

# Vector de probablidades
p<-c(0.45,0.45);p[3]<-1-sum(p)

# Número de ensayos
n<-10

# Muestra aleatoria multinomial de tamaño 1000
x<-t(rmultinom(10000,n,p))

# Matriz de densidades
d<-table(x[,1],x[,2])/100

# Histograma 3D
par(mar=c(2,2,2,2))
hist3D(z=d,col='aquamarine',border=1,phi = 30,theta=-25,d=2,bty='g',
       shade = .35, zlab = "Densidad",space = 0.1,cex.axis = 1e-9)

# Semilla aleatoria
set.seed(999)

# Vector de probablidades
p<-c(0.8,0.1);p[3]<-1-sum(p)

# Número de ensayos
n<-10

# Muestra aleatoria multinomial de tamaño 1000
x<-t(rmultinom(10000,n,p))

# Matriz de densidades
d<-table(x[,1],x[,2])/100

# Histograma 3D
par(mar=c(2,2,2,2))
hist3D(z=d,col='aquamarine',border=1,phi = 30,theta=-25,d=2,bty='g',
       shade = .35, zlab = "Densidad",space = 0.1,cex.axis = 1e-9)

Distribución Uniforme

Sea \(X\) una v.a continua definida en el intervalo \((a,b)\) tal que \(\forall_I \subset (a,b)\) la \(P(X \in I)\) es proporcional a la longitud del intervalo.

Entonces se dice que \(X\) tiene una distribución uniforme en el intervalo \((a,b)\). Para hallar la p.d.f. de X, como la probabilidad de que X pertenezca a cualquier subintervalo de \((a,b)\) es proporcional a su longitud, necesariamente la p.d.f. de X es una función constante;

\[ f(x) = \begin{cases} k & ; \quad a<x<b\\ 0 & ; \quad \text{en otro caso}\\ \end{cases} \]

\[P(X \in (a,b))=k(b-a)=1 \Rightarrow k=\cfrac{1}{b-a}\]

De esta manera, la p.d.f. de la variable aleatoria X es de la forma:

\[ f(x) = \begin{cases} \cfrac{1}{b-a} & ; \quad a<x<b\\ \\ 0 & ; \quad \text{en otro caso}\\ \end{cases} \]

\[ F(x) = \begin{cases} 0 & ; \quad x<a \\\\ \cfrac{x-a}{b-a} & ; \quad a \leq x \leq b \\\\ 1 & ; \quad x>b\\ \end{cases} \]

Entonces \(X \sim U(a,b)\), además:

\[M_X(t)= \begin{cases} \cfrac{\large e^{tb}-\large e^{ta}}{t(a-b)} \quad t \neq 0 \\ \\ \quad \quad 1 \quad \quad \quad t=0 \\ \end{cases}\]

Así, \(E[X]=\cfrac{a+b}{2}\) y \(\text{Var}[X]=\cfrac{(b-a)^2}{12}\)

Distribución Normal

La Distribución Normal o Gaussiana es indudablemente la más importante y la de mayor uso de todas las distribuciones de probabilidad continua. Es la piedra angular en la aplicación de la inferencia estadística en el analisis de datos, puesto que las distribuciones de muchas estadísticas muestrales tienden hacia la distribución normal conforme crece el tamaño de muestra.

Algunos fenómenos presentan este tipo de distribución:

• Precipitación pluvial
• Temperaturas
• Mediciones en organismos vivos
• Calificaciones de pruebas de actitud
• Mediciones físicas en manufactura
• Errores de instrumentación

Debe tenerse mucho cuidado al suponer para una situación dada un modelo de probabilidad normal sin previa comprobación. si bien es cierto que esta distribución es la de mayor uso, también es de la que más se abusa. Quizás por la mala interpretación de la palabra “normal”, especialmente si se aplica su significado literal de “patrón o estándar aceptado”. Es posible que una distribución normal proporcione de manera razonable una aproximación alrededor de una media de una variable aleatoria; sin embargo, puede resultar que para valores extremos se encuentren en cualquier dirección.

La apariencia gráfica de la distribución normal es una curva simétrica con forma de campana, que se extiende sin límite tanto en dirección positiva como negativa.

Propiedades

Sea \(X\) una variable aleatoria tal que que \(X \sim N(\mu,\sigma^2)\). Las siguientes afirmaciones se cumplen para la variable aleatoria \(X\).

• El área abajo la curva de una distribución normal comprendida entre \((- \infty ; \infty)\) es igual a 1.
• La normal es simétrica respecto a \(\mu\). Como consecuencia de esto, el área bajo de la curva en el intervalo \([\mu, \infty)\) es igual a 0.5.
• La distribución normal tiene forma de campana.
• La normal queda completamente caracterizada con el conocimiento de sus parámetros de localización y de escala.
• En la práctica la normalidad se alcanza de manera aproximada.
• \(E[X]=\mu\) y \(Var[X]=\sigma^2\).
• Cuanto más grande es el parámetro \(\sigma^2\) el gráfico de la función es mas “achatada” y de colas más largas.

Normal Estándar

Cuando \(\mu=0\) y \(\sigma^2=1\), se obtiene una distribución normal especial, conocida como Normal Estándar y usualmente es denotada cola letra \(Z\) y se escribe \(Z \sim N(0,1)\).

\[P(x_1 < X < x_2)= \int\limits_{x_1}^{x_2} \cfrac{1}{\sqrt{2 \pi} \sigma} \large e^{- \frac{1}{2} \frac{(x-\mu)^2}{\sigma^2}} \small dx\]

Esta probabilidad no puede ser resuelta de manera explícita, además cada vez que se modifique \(\mu\) o \(\sigma^2\), se debe calcular de nuevo la integral.

Normalización

Para evitar este problema, se utiliza el cambio de variable \(Z=\cfrac{X-\mu}{\sigma}\), del cual se obtiene que \(dz=\cfrac{1}{\sigma}dx\).

\[ \begin{aligned} P(x_1 < X < x_2) &= \int\limits_{z_1}^{z_2} \cfrac{1}{\sqrt{2 \pi}} \mathit{e}^{- \cfrac{z^2}{2}} dz \\ &= P(z_1 < Z < z_2)\\ \end{aligned} \]

donde \(z_1=\cfrac{x_1-\mu}{\sigma}\) y \(z_2=\cfrac{x_2-\mu}{\sigma}\). Esto indica que cualquier cálculo de probabilidades para una variable aleatoria \(N(\mu, \sigma^2)\) puede reducirse al cálculo de probabilidades con una variable aleatoria \(N(0,1)\). El proceso de transformar una v.a normal en una normal estándar se conoce como Estandarización o Normalización.

Para el cálculo de probabilidades con una variable \(N(0, 1)\), la c.d.f. estará dada por:

\[\Phi(z)=P(Z \leq z)=\int\limits_{-\infty}^{z} \cfrac{1}{\sqrt{2 \pi}} \large e^{- \frac{z^2}{2}} \small dz\]

De esta manera se tiene que:

\[P(x_1 < X < x_2)=P(z_1 < Z < z_2)=\Phi(z_2)-\Phi(z_1)\]

Los valores de la función \(\Phi(z)\) están tabulados en casi todos los textos de estadística.

Tabla Distribución Normal

Distribución Exponencial

La distribución exponencial tiene una gran utilidad práctica ya que podemos considerarla como un modelo adecuado para la distribución de probabilidad del tiempo de espera entre dos hechos que sigan un proceso de Poisson. De hecho la distribución exponencial puede derivarse de un proceso experimental de Poisson con las mismas características que las que enunciábamos al estudiar la distribución de Poisson, pero tomando como variable aleatoria, en este caso, el tiempo que tarda en producirse un hecho.

Sea \(X\) una v.a continua. \(X \sim Exp(\lambda)\).

\[f(x) = \begin{cases} \lambda \mathit{e}^{-\lambda x} &;\qquad x > 0\\ 0 &;\qquad \text{en otro caso}\\ \end{cases}\]

La función generadora de momentos está dada por:

\[M_X(t)=E\left[ e^{tx}\right]=\cfrac{\lambda}{\lambda-t}\]

Así,\(E[X]=\frac{1}{\lambda}\) y \(\text{Var}[X]=\frac{1}{\lambda^2}\).

La función acumulada está dada por:

\[F(x) = \begin{cases} 1 - \mathit{e}^{-\lambda x} &;\qquad x \geq 0\\ 0 &;\qquad x < 0\\ \end{cases}\]

Por lo general en la distribución exponencial se da el tiempo promedio hasta la ocurrencia de un evento Poisson, es útil igualar el promedio a \(E[X]\) para hallar el parámetro \(\lambda\) despejando de la ecuación.

Distribución Log-Normal

Una variable aleatoria \(X\), no negativa, tiene una p.d.f Lognormal si \(Y = ln(X)\) es una v.a con p.d.f normal. Haciendo \(\mu = E[Y]\) y \(\sigma^2 = Var[Y]\), la p.d.f de \(X\) es de la forma

\[f(x)=\cfrac{1}{x\sigma \sqrt{2 \pi}} \large e^{- \frac{1}{2} \frac{(\text{ln}(x)-\mu)^2}{\sigma^2}}; \quad \small \forall_x >0\]

\(\mu\) y \(\sigma^2\) NO son la media y varianza de \(X\). Son la media y varianza de \(\text{ln}(X)\).

\[E[X]=\mathit{e}^{\mu+\frac{\sigma^2}{2}} \quad\quad Var[X]=\mathit{e}^{2\mu+\sigma^2}\Big(\mathit{e}^{\sigma^2}-1\Big)\]

El cálculo de probabilidades con una p.d.f Lognormal es algo complicado. Pero debido al hecho de que el logaritmo natural de un v.a Lognormal es una v.a normal, podemos usar las tablas para una normal estándar para calcular dichas probabilidades. Como \(\text{ln}(X)\) es una función estrictamente creciente, entonces:

\[ \begin{aligned} P(X \leq a) &= P\left(\text{ln}[X] \leq \text{ln}[a]\right)=P\left(\cfrac{\text{ln}[X]-\mu}{\sigma} \leq \cfrac{\text{ln}[a]-\mu}{\sigma}\right) \\\\ &= P\left(Z \leq \cfrac{\text{ln}[a]-\mu}{\sigma}\right)= \Phi\left(\cfrac{\text{ln}[X]-\mu}{\sigma}\right) \end{aligned} \]

Así, la f.d.a de \(X\) es de la forma:

\[F(X)=P(X \leq x)=\Phi\left(\cfrac{\text{ln}[X]-\mu}{\sigma}\right); \quad \forall_x >0\]

Distribución Gamma

En distintos ámbitos interesa modela fenomenos que involucren conteos y tiempos al mismo tiempo. Hasta ahora hemos visto una aproximación con las distribuciones Exponencial y Poisson, siempre y cuando los parámetros sean proporcionales. ¿Qué pasaría si esto no se cumple?

Conocer el comportamiento de muchos fenómenos que relacionen el tiempo supervivencia (vida util) entre un número de ocurrencias de un evento, se de gran interés. En este caso es muy importante conocer una distribución que nos permita conocer las probabilidades, un valor esperado y una variabilidad de la ocrurrencia de el evento de interés y para ello vamos a considerar la siguiente distribución.

Esta distribución es usualmente usada en el análisis de superivencia, en confiabilidad y demás áreas que tienen un enfoque relacionado en observar el tiempo de falla o de vida de una persona, duración de productos industriales, tiempo que transcurre un número de personas en una fila.

Recuerde que la función Gamma está dada por:

\[\Gamma(\alpha)=\int_{0}^{\infty}x^{\alpha-1}e^{-x}dx \qquad x>0,\quad a>0 \]

Sea \(X\) una variable aletoria continua tal que \(X\sim Gamma(\alpha,\beta)\). Así la f.d.p está dada por:

\[f_{(\alpha,\beta)}(x)= \frac{1}{\Gamma(\alpha)\beta^{\alpha}}x^{\alpha-1}e^{-x/\beta} \quad \alpha,\beta,x>0\]

• El parámetro \(\alpha\) se le conoce como de escala y se puede interpretar como el número necesario de veces que se debe realizar un experimento hasta observar el evento de interés hasta un tiempo \(t\), por ejemplo se requieren 2 ciclos de un proceso para que una pieza falle.

• Por otro lado \(\beta\) se le conoce como un parámetro de forma y se interpreta como el número de ocurrencias en una unidad de tiempo o intervalo específico, por ejemplo 3 piezas falla cada 100 horas.

La probabilidad de llegar hasta un valor de x \(P(X<x)\) es decir, la f.d.a es:

\[P(X<x)=\int_{0}^{x}\frac{1}{\Gamma(\alpha)\beta^{\alpha}}t^{\alpha-1}e^{-t/\beta}dt \quad \alpha,\beta,t>0\]

\[M_{x}(t)=(1-\beta t)^{-\alpha}, \quad t< \frac{1}{\beta}\]

Entonces \(E[x]=\alpha \beta\) y \(\text{Var}[X]=\alpha\beta^{2}\).

par(mfrow=c(2,3))
curve(dgamma(x,shape=2,scale=1/4),xlim=c(0,5),lwd=2,ylab='Densidad',
     main=expression("Distribución Gamma "*alpha==2))
curve(dgamma(x,shape = 2,scale = 1/2),lty=2,lwd=2,add=T)
curve(dgamma(x,shape = 2,scale = 1),lty=4,lwd=2,add=T)
legend('topright',col=1,lty=c(1,2,4),lwd=2,legend=c(expression(beta==1/4),
                                                    expression(beta==1/2),
                                                    expression(beta==1)))

curve(dgamma(x,shape=3,scale=1/4),xlim=c(0,5),col=2,lwd=2,ylab='Densidad',
     main=expression("Distribución Gamma "*alpha==3))
curve(dgamma(x,shape = 3,scale = 1/2),lty=2,col=2,lwd=2,add=T)
curve(dgamma(x,shape = 3,scale = 1),lty=4,col=2,lwd=2,add=T)
legend('topright',col=2,lty=c(1,2,4),lwd=2,legend=c(expression(beta==1/4),
                                                    expression(beta==1/2),
                                                    expression(beta==1)))

curve(dgamma(x,shape = 4,scale = 1/4),xlim=c(0,6),col=4,lwd=2,ylab='Densidad',
     main=expression("Distribución Gamma "*alpha==4))
curve(dgamma(x,shape = 4,scale = 1/2),lty=2,col=4,lwd=2,add=T)
curve(dgamma(x,shape = 4,scale = 1),lty=4,col=4,lwd=2,add=T)
legend('topright',col=4,lty=c(1,2,4),lwd=2,legend=c(expression(beta==1/4),
                                                    expression(beta==1/2),
                                                    expression(beta==1)))

curve(dgamma(x,shape = 2,scale = 1/4),xlim=c(0,4),lwd=2,ylab='Densidad',
     main=expression("Distribución Gamma "*beta==1/4))
curve(dgamma(x,shape = 3,scale = 1/4),col=2,lwd=2,add=T)
curve(dgamma(x,shape = 4,scale = 1/4),col=4,lwd=2,add=T)
legend('topright',col=c(1,2,4),lty=1,lwd=2,legend=c(expression(alpha==2),
                                                    expression(alpha==3),
                                                    expression(alpha==4)))

curve(dgamma(x,shape = 2,scale = 1/2),xlim=c(0,5),lwd=2,lty=2,ylab='Densidad',
     main=expression("Distribución Gamma "*beta==1/2))
curve(dgamma(x,shape = 3,scale = 1/2),col=2,lwd=2,lty=2,add=T)
curve(dgamma(x,shape = 4,scale = 1/2),col=4,lwd=2,lty=2,add=T)
legend('topright',col=c(1,2,4),lty=2,lwd=2,legend=c(expression(alpha==2),
                                                    expression(alpha==3),
                                                    expression(alpha==4)))

curve(dgamma(x,shape = 2,scale = 1),xlim=c(0,7),lwd=2,lty=4,ylab='Densidad',
     main=expression("Distribución Gamma "*beta==1))
curve(dgamma(x,shape = 3,scale = 1),col=2,lwd=2,lty=4,add=T)
curve(dgamma(x,shape = 4,scale = 1),col=4,lwd=2,lty=4,add=T)
legend('topright',col=c(1,2,4),lty=4,lwd=2,legend=c(expression(alpha==2),
                                                    expression(alpha==3),
                                                    expression(alpha==4)))

Distribución Beta

En muchos casos no es posble encontrar proporciones constantes en un fenómeno como el caso de una distribución uniforme, sin embargo hay una distribución que nos permite darle forma y controlar ciertos parámetros para obtener una f.d.p que es “conocida” y tiene una forma matematicamente cerrada, a esta distribución se le conoce como la distribución Beta.

Hoy en día la distribución beta ha tomado mucha importancia en el campo de la estadística bayesiana tomándose como distribución a priori para modelar proporciones por ejemplo de una distribución binomial, realizar estimaciones de una prevalencia o para modelar índices de riesgos.

Sea \(X\) una v.a continua tal que \(X\sim Beta(\alpha,\beta)\) con dominio solo para \(x\in[0,1]\)

Cuya f.d.p viene dada por:

\[f_{(\alpha,\beta)}(x)= \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\beta^{\alpha}}x^{\alpha-1}(1-x)^{\beta-1} \quad \alpha,\beta>0; \quad 0<x<1\]

La probabilidad de estar entre \(a\) y \(b\) con \(0<a<b<1\) será:

\(P(a<X<b)= \int_{a}^{b} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\beta^{\alpha}}x^{\alpha-1}(1-x)^{\beta-1}dt \quad \alpha,\beta>0; \quad 0<x<1\)

La probailidad de llegar hasta un valor x \(P(X<x)\) la f.d.a se calcula como:

\(P(X<x)= \int_{0}^{x} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\beta^{\alpha}}t^{\alpha-1}(1-t)^{\beta-1}dt \quad \alpha,\beta>0; \quad 0<t<1\)

\[M_{x}(t)=1+\sum_{k=1}^{\infty}\left(\prod_{r=0}^{k-1}\frac{\alpha+r}{\alpha+\beta+r}\right)\frac{t^{k}}{k!} \]

Entonces \(E[X]=\frac{\alpha}{\alpha+\beta}\) y \(\text{Var}[X]=\frac{\alpha \beta}{(\alpha+\beta)^{2}(\alpha+\beta+1)}\).

par(mfrow=c(1,3))
curve(dbeta(x,1,1),xlim=c(0,1),ylim=c(0,3),ylab='Densidad',lwd=2,
      main=expression('Distribución Beta '*alpha==beta))
curve(dbeta(x,2,2),col=2,lwd=2,add=T)
curve(dbeta(x,4,4),col=3,lwd=2,add=T)
curve(dbeta(x,7,7),col=4,lwd=2,add=T)
curve(dbeta(x,0.5,0.5),col=6,lwd=2,add=T)
legend(x=0.3,y=0.6,lty=1,col=c(1:4,6),lwd=2,
       legend=c(expression(alpha*'='*beta==1),expression(alpha*'='*beta==2),
                expression(alpha*'='*beta==4),expression(alpha*'='*beta==7),
                expression(alpha*'='*beta==0.5)))

curve(dbeta(x,2,2),xlim=c(0,1),ylim=c(0,4),ylab='Densidad',lwd=2,
      main=expression('Distribución Beta '*alpha<=beta))
curve(dbeta(x,2,4),col=2,lwd=2,add=T)
curve(dbeta(x,2,7),col=3,lwd=2,add=T)
curve(dbeta(x,2,9),col=4,lwd=2,add=T)
legend('toprigh',lty=1,col=1:4,lwd=2,title=expression(alpha==2),
       legend=c(expression(beta==2),expression(beta==4),
                expression(beta==7),expression(beta==9)))

curve(dbeta(x,2,2),xlim=c(0,1),ylim=c(0,4),ylab='Densidad',lwd=2,
      main=expression('Distribución Beta '*alpha>=beta))
curve(dbeta(x,4,2),col=2,lwd=2,add=T)
curve(dbeta(x,7,2),col=3,lwd=2,add=T)
curve(dbeta(x,9,2),col=4,lwd=2,add=T)
legend('topleft',lty=1,col=1:4,lwd=2,title=expression(beta==2),
       legend=c(expression(alpha==2),expression(alpha==4),
                expression(alpha==7),expression(alpha==9)))

Distribución \(\chi^2\)

La distribución de Pearson, llamada también ji cuadrada(o) o chi cuadrado(a) (\(\chi^2\)), es una distribución de probabilidad continua con un parámetro \(r\) que representa los grados de libertad de la variable aleatoria.

\[ f(x)= \begin{cases}\displaystyle \frac{1}{2^{r/2}\Gamma(r/2)}\,x^{(r/2) - 1} e^{-x/2}&\text{para }x>0,\\ 0&\text{en otro caso} \end{cases} \]

curve(dchisq(x,2),lwd=2,ylab='Densidad',xlim=c(0,25),ylim=c(0,0.3),col='salmon',
     main = expression('Distribución '*chi[(r)]^2))
curve(dchisq(x,3),col=1,lwd=2,add=T)
curve(dchisq(x,5),col=2,lwd=2,add=T)
curve(dchisq(x,7),col=3,lwd=2,add=T)
curve(dchisq(x,9),col=4,lwd=2,add=T)
curve(dchisq(x,11),col=6,lwd=2,add=T)
legend('topright',legend=c('r = 2','r = 3','r = 5','r = 7','r = 9','r = 11'),
       col=c('salmon',1:4,6),lty=1,lwd=2,title='Grados de Libertad')

Distribución t-student

La distribución t-de Student es una distribución de probabilidad que surge del problema de estimar la media de una población normalmente distribuida cuando el tamaño de la muestra es pequeño. Su f.d.p está dada por:

\[f(x) = \frac{\Gamma\left(\dfrac{\nu + 1}{2}\right)}{\Gamma\left(\dfrac{\nu}{2}\right)\sqrt{\pi\nu}\sigma} \left[1+\frac{1}{\nu}\left(\frac{x-\mu}{\sigma}\right)^2\right]^{-\tfrac{\nu+1}{2}}\]

curve(dnorm(x),xlim=c(-5,5),main='Distribución t-Student',ylab='Densidad',lwd=2)
curve(dt(x,df = 1),col=2,add=T,lwd=2)
curve(dt(x,df = 2),col=3,add=T,lwd=2)
curve(dt(x,df = 5),col=4,add=T,lwd=2)

legend('toprigh',lty=1,col=1:4,lwd=2,cex=0.9,
       legend=c('Normal (0,1)','Cauchy','t(v=2)','t(v=5)'))

En esta gráfica podemos apreciar la realción que existe entre las distribuciones Cauchy, t-de Student y la Normal Estandarizada. Se observa que la distribución t-de Student con \(v = 1\) grado de libertad es una distribución Cauchy la cual es la más achatada y presenta las colas mas pesadas en comparación las demás distribuciones presentes en la gráfica. Además se puede evidenciar que a medida que aumentan los grados de libertad en la distribución t, ésta tiende a la Normal.

Distribución Normal Bivariada

Una distribución normal multivariante, también llamada distribución gaussiana multivariante, es una generalización de la distribución normal unidimensional a dimensiones superiores.

Para el caso bivariado el vector aleatorio \((X,Y)\), que toma todos lo valores en el plano euclidiano, tiene la siguiente función de distribución de probabilidad:

\[f(x,y)=\frac{1}{2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp \left\{ -\frac{1}{2 (1-\rho^2)} \left[ \left(\frac{x-\mu_x}{\sigma_x} \right)^2 - 2 \rho \frac{(x-\mu_x)(y-\mu_y)}{\sigma_x \sigma_y}+ \left(\frac{y-\mu_y}{\sigma_y} \right)^2 \right] \right\}\]

x<-seq(-6,6,length=40);y<-x

# Función Normal Bivariada
f<-function(x,y,mu=c(0,0),s=c(1,0,1))
{
  mu1<-mu[1];mu2<-mu[2]
  s11<-s[1];s12<-s[2];s22<-s[3];rho<-s12/(s11*s22)
  term1<-1/(2*pi*sqrt(s11*s22*(1-rho^2)))
  term2<--1/(2*(1-rho^2))
  term3<-(x-mu1)^2/s11
  term4<-(y-mu2)^2/s22
  term5<--2*rho*((x-mu1)*(y-mu2))/(sqrt(s11)*sqrt(s22))
  term1*exp(term2*(term3+term4-term5))
} 

z<-outer(x,y,f,s=c(5,0,5)) # matriz de densidades

persp(x, y, z,col="lightsalmon",theta=30, phi=20,r=50,d=0.1,expand=0.5,ltheta=90,
      main="Distribución Normal Bivariada",lphi=180,shade=0.3,ticktype="detailed",
      nticks=5,cex.axis=0.7,zlab = 'Densidad')
mtext(expression(list(mu[x]==0,mu[y]==0,sigma[x]==5,sigma[y]==5,
                      sigma[xy]==0,rho==0)), side=3)

x<-seq(-6,6,length=40);y<-x
z<-outer(x,y,f,s=c(5,12.5,5)) # matriz de densidades

persp(x, y, z,col="lightsalmon",theta=30,phi=20,r=50,d=0.1,expand=0.5,ltheta=90,
      main="Distribución Normal Bivariada",lphi=180,shade=0.3,ticktype="detailed",
      nticks=5,cex.axis=0.7,zlab = 'Densidad')
mtext(expression(list(mu[x]==0,mu[y]==0,sigma[x]==5,sigma[y]==5,
                      sigma[xy]==12.5,rho==0.5)), side=3)

x<-seq(-6,6,length=40);y<-x
z<-outer(x,y,f,s=c(5,-12.5,5)) # matriz de densidades

persp(x, y, z,col="lightsalmon",theta=30, phi=20,r=50,d=0.1,expand=0.5,ltheta=90,
      main="Distribución Normal Bivariada",lphi=180,shade=0.3,ticktype="detailed",
      nticks=5,cex.axis=0.7,zlab = 'Densidad')
mtext(expression(list(mu[x]==0,mu[y]==0,sigma[x]==5,sigma[y]==5,
                      sigma[xy]==-12.5,rho==-0.5)), side=3)

x<-seq(-6,6,length=40);y<-x
z<-outer(x,y,f,s=c(5,24,5)) # matriz de densidades

persp(x, y, z,col="lightsalmon",theta=30, phi=20,r=50,d=0.1,expand=0.5,ltheta=90,
      main="Distribución Normal Bivariada",lphi=180,shade=0.3,ticktype="detailed",
      nticks=5,cex.axis=0.7,zlab = 'Densidad')
mtext(expression(list(mu[x]==0,mu[y]==0,sigma[x]==5,sigma[y]==5,
                      sigma[xy]==24,rho==0.96)), side=3)

Solución usando R

Lo primero es descargar la base de datos en la carpeta que se considere, luego debemos cargar la base de datos en RStudio para realizar los análisis. Antes de comenzar con el análisis necesitamos instalar o cargar las librerías necesarias para el análisis

# Para manipular datos
if(!require(dplyr)) install.packages("dplyr")
# Gráficos con estilo
if(!require(ggplot2)) install.packages("ggplot2")
# Pruebas de bondad de ajuste
if(!require(nortest)) install.packages("nortest")
# Gráficos QQ-plot con IC
if(!require(car)) install.packages("car")
# Importar/Exportar datos de Excel
if(!require(openxlsx)) install.packages("openxlsx")

El siguiente paso es cargar la base de datos, esto lo podemos realizar directamente de RStudio con el siguiente comando:

Licores <- read.xlsx(file.choose(),sheet = 1)

Licores <- read.xlsx("Licores.xlsx",sheet = 1)

Revisemos la estructura de los datos, para verificar si fueron cargados correctamente.

str(Licores)

## 'data.frame':    300 obs. of  5 variables:
##  $ TL : chr  "Ron" "Aguardiente" "Ron" "Aguardiente" ...
##  $ PI : num  34719 26900 27062 29317 32471 ...
##  $ GAE: num  35 29 35 29 38 29 29 38 35 40 ...
##  $ GAQ: num  30.2 29 32.2 29 32.6 ...
##  $ CE : num  221 253 249 254 267 ...

# cambiando caracter por factor
Licores <- Licores %>% mutate_if(is.character,as.factor)

El primer punto nos piden calcular el precio total de incautación, y se pregunta si el promedio de esa nueva variable es superior a 7.5 millones. Entonces tenemos las hipótesis:

\[H_o: \mu_{pt} = 7'500.000 \qquad H_a:\mu_{pt} > 7'500.000\]

Para seleccionar el estadístico de prueba más indicado, debemos probar si la distribución del precio total de incautación es Normal.

\[H_0:PT \sim N(\mu,\sigma^2) \qquad H_a: PT\not\sim N(\mu,\sigma^2)\]

Licores$PT <- Licores$CE*Licores$PI
hist(Licores$PT,col="lightblue",las=1)

# pruebas de hipótesis
shapiro.test(Licores$PT)

## 
##  Shapiro-Wilk normality test
## 
## data:  Licores$PT
## W = 0.99485, p-value = 0.4134

ad.test(Licores$PT)

## 
##  Anderson-Darling normality test
## 
## data:  Licores$PT
## A = 0.25534, p-value = 0.7249

lillie.test(Licores$PT)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Licores$PT
## D = 0.032027, p-value = 0.6393

Entonces, siguiendo el árbol de decisión tenemos que los datos se distribuyen normal y no conocemos los parámetros poblacionales, por lo tanto el estadístico de prueba está basado en la t-student.

t.test(Licores$PT,mu = 7500000,alternative = "greater",conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  Licores$PT
## t = -0.42571, df = 299, p-value = 0.6647
## alternative hypothesis: true mean is greater than 7500000
## 95 percent confidence interval:
##  7348999     Inf
## sample estimates:
## mean of x 
##   7469031

Como \(\text{valor-p} \, > \, \alpha\) entonces no se rechaza \(H_0\), existe evidencia muestral suficiente de que el promedio del precio de incautaciones mesual es de $7’500.000, esto implica anualmente se esperaría que alrededor de 90 millones de pesos deje de circular ilegalmente gracias a lo no pago de impuestos o adulteración de licores.

En la siguiente pregunta se pide calcular el ipoconsumo de acuerdo a una formulación, y luego comparar si el impuesto de los licores de tipo Ron son menores a los de Whisky basado en los promedios. Entonces tenemos la siguiente hipótesis:

\[H_0: \mu_{Ir} = \mu_{Iw} \qquad H_a:\mu_{Ir} < \mu_{Iw}\]

De acuerdo al árbol de decisión tenemos que hacer los siguientes pasos:

• Calcular el ipoconsumo
• Separar la base de datos en Ron y Whisky
• Probar la normalidad de ambos grupos
• Comprobar la igualdad de varianzas de los grupos
• Escoger el mejor estadístico de prueba y realizar el análisis.
• Concluir

Los primeros dos pasos son sencillos:

Licores <- Licores %>%
   mutate(IPO = GAQ*CE*400)
hist(Licores$IPO,col="lightblue",las=1)

SoloRW <- Licores %>% 
   filter(TL %in% c("Ron","Whisky")) %>%
   droplevels
SoloRW %>% group_by(TL) %>% 
   summarise(n=length(IPO),
             media=mean(IPO),
             desv=sd(IPO))

SoloRW %>% ggplot(aes(x=IPO,fill=TL))+
   geom_density()+xlim(c(2e06,5e06))+theme_bw()

Luego debemos probar la normalidad de los datos, esto lo podemos hacer de dos formas. La primera es separar la base de datos y construir dos nuevas con los grupos de interés, la otra es evaluar al mismo tiempo los valores-p de las pruebas usando dplyr.

# Primer método - Ron
SoloRW %>% filter(TL=="Ron") %>%
   dplyr::select(IPO) %>% unlist %>% shapiro.test()

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.96894, p-value = 0.1904

SoloRW %>% filter(TL=="Ron") %>%
   dplyr::select(IPO) %>% unlist %>% ad.test()

## 
##  Anderson-Darling normality test
## 
## data:  .
## A = 0.50493, p-value = 0.1943

SoloRW %>% filter(TL=="Ron") %>%
   dplyr::select(IPO) %>% unlist %>% lillie.test()

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  .
## D = 0.095664, p-value = 0.2758

# Primer método - Whisky
SoloRW %>% filter(TL=="Whisky") %>%
   dplyr::select(IPO) %>% unlist %>% shapiro.test()

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.96678, p-value = 0.5193

SoloRW %>% filter(TL=="Whisky") %>%
   dplyr::select(IPO) %>% unlist %>% ad.test()

## 
##  Anderson-Darling normality test
## 
## data:  .
## A = 0.26387, p-value = 0.671

SoloRW %>% filter(TL=="Whisky") %>%
   dplyr::select(IPO) %>% unlist %>% lillie.test()

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  .
## D = 0.087079, p-value = 0.8639

# Segundo método
SoloRW %>% group_by(TL) %>%
   summarise(valor_p=shapiro.test(IPO)$p.val)

SoloRW %>% group_by(TL) %>%
   summarise(valor_p=ad.test(IPO)$p.val)

SoloRW %>% group_by(TL) %>%
   summarise(valor_p=lillie.test(IPO)$p.val)

par(mfrow=c(1,2))
with(SoloRW %>% filter(TL=="Ron"),
     qqPlot(IPO,pch=19,main = "QQ-Plot Ron",id=F,col.lines = 2))
with(SoloRW %>% filter(TL=="Whisky"),
     qqPlot(IPO,pch=19,main = "QQ-Plot Whisky",id=F,col.lines = 4))

Como ambas muestras distribuyen normal, es necesario verificar si existe igualdad en sus varianzas, esto se comprueba con la función var.test.

with(SoloRW,var.test(IPO~TL))

## 
##  F test to compare two variances
## 
## data:  IPO by TL
## F = 0.73017, num df = 51, denom df = 26, p-value = 0.3333
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.3562148 1.3878542
## sample estimates:
## ratio of variances 
##           0.730172

Entonces, como los datos son normales y tienen varianzas iguales, de acuerdo con el árbol de decisión el mejor estadístico de prueba está basado en la t-student y es necesario estimar la varianza conjunta, no obstante, el software lo hace:

with(SoloRW,
     t.test(IPO[TL=="Ron"],
            IPO[TL=="Whisky"],
            var.equal = T,alternative = "less"))

## 
##  Two Sample t-test
## 
## data:  IPO[TL == "Ron"] and IPO[TL == "Whisky"]
## t = -5.9561, df = 77, p-value = 3.629e-08
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##       -Inf -292527.9
## sample estimates:
## mean of x mean of y 
##   3185673   3591694

En conclusión, el ipoconsumo del ron es menor al del whisky con un nivel de significancia del 5%. Entonces la cantidad de impuestos que deja de percibir el departamento tiene mayor aporte en el Whisky, sin embargo, las cantidades de Ron son casi el doble que las de Whisky lo que sugiere que este licor es preferido para la adulteración. Los controles se deben aumentar para incautar más Ron y establecer una estrategia para que no aumente la adulteración del Whisky que es la que más afecta los recursos para salud.

El tercer punto nos piden calcular una nueva variable, si \(GAQ=GAE\) el licor se clasifica como contrabando en caso contrario se llamará fraudulento.

Licores <- Licores %>%
   mutate(TipoLicor = factor(ifelse(GAQ==GAE,"Contrabando","Fraudulento")))

Luego preguntan si la proporción de licores Fraudulentos es superior al 65%. Entonces las hipótesis son:

\[H_0: p_F = 0.65 \qquad H_a: p_F >0.65\]

with(Licores,table(TipoLicor))

## TipoLicor
## Contrabando Fraudulento 
##          83         217

Como en R los factores se organizan alfabéticamente, entonces debemos cambiar la hipótesis hacia los licores de contrabando, así:

\[H_0: p_C = 0.35 \qquad H_a: p_C < 0.35\]

with(Licores,
     prop.test(table(TipoLicor),
               p = 0.35,
               alternative = "less",
               conf.level = 0.95,
               correct = F))

## 
##  1-sample proportions test without continuity correction
## 
## data:  table(TipoLicor), null probability 0.35
## X-squared = 7.0916, df = 1, p-value = 0.003872
## alternative hypothesis: true p is less than 0.35
## 95 percent confidence interval:
##  0.0000000 0.3210025
## sample estimates:
##         p 
## 0.2766667

Como el valor-p es mucho menor al nivel de significancia, entonces se rechaza \(H_0\) por lo tanto existe evidencia muestral para afirmar que la proporción de licores Fraudulentos es mayor del 65%. Por lo tanto, las rentas ilegales de licores están basadas en mayor medida de la adulteración, el impacto para la salud pública es muy grave, pues los efectos nocivos de los licores que no son correctamente destilados aumenta los riesgos clínicos como la ceguera, el cáncer de órganos digestivos o la muerte.