Colección de ejercicios para resolver de forma tradicional con ayuda de R y RStudio.

Grupo 10. Marcos Flórez Fernández, Julio González Escudero, Marina Martínez González, Lucía Ojer Guerra, Roberto Panero Hoz, Sergi Sos Lluch.

EJERCICIO 1. MODELOS DE REGRESIÓN PARA LOS CASOS ACUMULADOS DE COVID-19 EN REINO UNIDO

Un modelo de regresión es una herramienta matemática que nos ayuda a determinar una relacíon entre una variable dependiente con respecto a otras variables llamadas explicativas o independientes. Estos, a lo largo de los últimos meses, se han puesto muy de moda en los noticiarios y periódicos de todo el mundo, ya que nos permite representar y estudiar la evolución del virus en un país, en otras palabras, hacer un seguimiento y predicciones, además de proporcionar gráficas de “fácil” comprensión para que la gente aprenda como funciona, desde un punto de vista matemático, la propagación del virus.

En nuestro caso estudiaremos la evolución de los casos acumulados de nuevas infecciones que han ocurrido en Reino Unido (UK).

Apartado a): Modelar el número acumulado de personas infectadas mediante un ajuste lineal.
Para ello, recurriremos a los comando que R nos proporciona. Empezaremos, lo primero de todo, por introducir los datos oficales:

days<-read.table("AcummulativeDaysUK")
NewCasesUK<-read.table("NewCasesUK",sep=".")

days<-c(days$V1);days

 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
[45] 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
[89] 89 90 91 92 93

NewCasesUK<-c(NewCasesUK$V1);NewCasesUK

 [1]    2    0    0    0    0    1    0    0    1    4    0    0    1    0    0    0    0    0    0    0    0    0    0    4    0    0
[27]    0    6    4   12    5   11   34   29   46   46   65   50   52   83  134  207  264  330  152  407  676  643  714 1040  665  967
[53] 1430 1450 2130 2890 2560 2500 2670 3250 4570 4520 4670 4000 6200 4140 3890 5870 4680 5710 5230 5290 4340 5250 4600 4620 5600 5530
[79] 5850 4680 4300 4450 4580 5390 4910 4460 4310 4000 4080 6030 6200 4810 4340

Como podemos observar, la primera lista de números, se trata de el número de días que han pasado desde el día 1 correspondiente con 31/1/2020, cuando se detectó el primer infectado. La segunda lista se corresponde con el número de nuevos infectados en cada día desde el día uno hasta hoy.
Para visualizarlo mejor, representaremos en una gráfica los picos de infecciones por día:

plot(days,NewCasesUK,xlim=c(0,100),ylim=c(0,7000),ylab="Number of Cases",xlab="Days since 31-1-20", main="New Cases UK")
lines(spline(days,NewCasesUK),lwd="3",col="blue")

Ahora vamos a realizar el ajuste lineal de la suma acumulada por día de los datos representados en la siguente gráfica:

AllCasesUK<-cumsum(NewCasesUK)
AllCasesUK

 [1]      2      2      2      2      2      3      3      3      4      8      8      8      9      9      9      9      9      9
[19]      9      9      9      9      9     13     13     13     13     19     23     35     40     51     85    114    160    206
[37]    271    321    373    456    590    797   1061   1391   1543   1950   2626   3269   3983   5023   5688   6655   8085   9535
[55]  11665  14555  17115  19615  22285  25535  30105  34625  39295  43295  49495  53635  57525  63395  68075  73785  79015  84305
[73]  88645  93895  98495 103115 108715 114245 120095 124775 129075 133525 138105 143495 148405 152865 157175 161175 165255 171285
[91] 177485 182295 186635

Si nos fijamos, hasta día de hoy el número total de infectados es de 177485.

plot(days,AllCasesUK,ylab="Number of Cases",xlab="Days since 31-1-20", main="New Cases UK")
lines(spline(days,AllCasesUK),col="red",lwd="3")

Lo que buscamos es minimizar el error cuadrático medio de los puntos de esta gráfica y modelarlo mediante la expresión de una recta de regresión: \[y=a+bx\]
Para ello usaremos principalmente el comando “lineal model” o lm() de R y las herramientas que nos proporciona para obtener los parámetros \(a\) y \(b\):

LinealModel<-lm(AllCasesUK~days)
summary(LinealModel)

Call:
lm(formula = AllCasesUK ~ days)

Residuals:
   Min     1Q Median     3Q    Max 
-43570 -28214  -1559  26043  61752 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -46280       6387  -7.246 1.35e-10 ***
days            1840        118  15.598  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 30550 on 91 degrees of freedom
Multiple R-squared:  0.7278,    Adjusted R-squared:  0.7248 
F-statistic: 243.3 on 1 and 91 DF,  p-value: < 2.2e-16

Hemos usado el comando summary() para mostrar algunos datos aparte de los parámetros \(a\) y \(b\) que también son de interés. Más concretamente nos referimos a “Multiple R-squared”, lo cual nos da el coeficiente de determinación (número comprendido entre 0 y 1 que de una forma general nos dice como de bueno es el ajuste), que en este caso es de 0.7145. Dado que no supera el 0.8, desgraciadamente, no podemos decir que nos encontramos ante un buen ajuste en relación al conjunto de datos proporcionado. Esto se debe a que en el intervalo de tiempo que hemos representado, el crecimiento y expansión del virus aumenta drásticamente entre los 40 y 60 dias, perdiendo la forma lineal en cuanto al aumento de casos.

A continuación haremos uso de los comandos de representación gráfica de R para mostrar la siguiente recta:

\[y=-46280+1840x\]

xaj<-seq(0.005,100,0.001)
plot(days,AllCasesUK,ylab="Number of Cases",xlab="Days since 31-1-20", main="Total Cases UK")
abline(LinealModel, col="orange",lwd="3")

Apartado b): Modelar el número acumulado de personas infectadas mediante un ajuste no lineal de los presentados en las transparencias teóricas de la asignatura.

Para este apartado hemos elegido el modelo exponencial, al cual hay que aplicarle un cambio de variable para linealizarlo, obtener los parámetros correspondientes y por último deshacer el cambio para poder representar el grafo. A continuación, la manipulación algebraica del modelo elegido:
\[ \begin{align*} y&=\alpha e^{\beta x}\\ \\ log(y)&=log(\alpha)+\beta x\\ \\ Z&=A+Bx \end{align*} \] \[Z=log(y);\ A=log(\alpha);\ B=\beta \]
Con la forma de esta recta averiguaremos los parámetros \(\alpha\) y \(\beta\) procediendo de la misma manera que en el apartado anterior con los comandos de R. La diferencia es que ahora tendremos que transformar los valores de \(y\) para la nueva variable \(Z\):

Z<-log(AllCasesUK);Z

 [1]  0.6931472  0.6931472  0.6931472  0.6931472  0.6931472  1.0986123  1.0986123  1.0986123  1.3862944  2.0794415  2.0794415  2.0794415
[13]  2.1972246  2.1972246  2.1972246  2.1972246  2.1972246  2.1972246  2.1972246  2.1972246  2.1972246  2.1972246  2.1972246  2.5649494
[25]  2.5649494  2.5649494  2.5649494  2.9444390  3.1354942  3.5553481  3.6888795  3.9318256  4.4426513  4.7361984  5.0751738  5.3278762
[37]  5.6021188  5.7714411  5.9215784  6.1224928  6.3801225  6.6808547  6.9669671  7.2377782  7.3414839  7.5755847  7.8732171  8.0922394
[49]  8.2897906  8.5217826  8.6461140  8.8031237  8.9977658  9.1627245  9.3643482  9.5856899  9.7477106  9.8840499 10.0116691 10.1478053
[61] 10.3124465 10.4523312 10.5788526 10.6757924 10.8096269 10.8899571 10.9599749 11.0571403 11.1283653 11.2089107 11.2773930 11.3421965
[73] 11.3923949 11.4499324 11.4977611 11.5436001 11.5964851 11.6461005 11.6960384 11.7342674 11.7681489 11.8020440 11.8357695 11.8740555
[85] 11.9077003 11.9373105 11.9651151 11.9902460 12.0152450 12.0510841 12.0866414 12.1133815 12.1369101

cor(Z,days)**2

[1] 0.9581646

Como podemos comprobar, el coeficiente de determinación de este conjunto de datos, es bastante mejor que en el caso del ajuste mediante la recta de regresión.

ExpModel<-lm(Z~days)
summary(ExpModel)

Call:
lm(formula = Z ~ days)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.81048 -0.67971  0.06022  0.70234  1.18595 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.034520   0.177367   0.195    0.846    
days        0.149601   0.003277  45.653   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8483 on 91 degrees of freedom
Multiple R-squared:  0.9582,    Adjusted R-squared:  0.9577 
F-statistic:  2084 on 1 and 91 DF,  p-value: < 2.2e-16

Tras haber realizado el “linear model” de los datos Z sobre “days” (datos del eje x), obtenemos los valores de los parámetros del ajuste lineal, ahora toca deshacer el cambio:
\[\alpha=e^A;\ \beta=B \]

A<-ExpModel1$coefficients[1]
b<-ExpModel1$coefficients[2]
a<-exp(A)
a

(Intercept) 
   1.035123 

b

     days 
0.1496007 

Con estos valores para los parámetros \(\alpha\) y \(\beta\) respectivamente dibujaremos la siguiente función:
\[y=1.035123 e^{0.1496007}\]

xaj<-xaj<-seq(0.005,100,0.001)
yaj<-a*exp(b*xaj)
plot(days,AllCasesUK,ylab="Number of Cases",xlab="Days since 31-1-20", main="Total Cases UK")
lines(spline(xaj,yaj),col="blue",lwd="4")

Lo que podemos observar en este ajuste es bastante gráfico. La función exponencial se ajusta mucho mejor que la recta (sobretodo en los 50 primeros dias), pero aún así, después, hay mucha disparidad entre los valores reales y los que la función nos da, ya que esta crece mucho más rápido que el número de infectados.

Apartado c): Mejorar dicho ajuste mediante técnicas de mínimos cuadrados no lineales.

Las técnicas de mínimos cuadrados no lineales que ya hemos ido presentando en los dos apartados anteriores son, a “grosso modo”, un problema de optimización muy complejo. Se trata de aproximar un conjunto de datos de forma lineal y después, refinar los parámetros mediante métodos iterativos:
\[ \min_{\theta} SR(\hat{\theta})= \min_{\theta}\sum_{t=1}^{T}\hat{u_t}(\hat{\theta})= \min_{\theta}\sum_{t=1}^{T}[y_t-f(x_t,\beta)]^2 \]
lo que implica resolver el sistema de ecuaciones,
\[ \left( \dfrac{\partial f(x_t,\beta)}{\partial \beta} \right)^{´}y= \left( \dfrac{\partial f(x_t,\beta)}{\partial \beta} \right)^{'}f(X,\beta) \]
donde el vector gradiente es \(Txk\), y \(f(X,\beta)\) es \(Tk1\). Este sistema puede no tener solución o tener múltiples soluciones. A diferencia del estimador de Mínimos Cuadrados aplicado a un modelo lineal, el estimador no es insesgado. La matriz de covarianzas del estimador resultante es:
\[ Var(\hat{\theta}) = \sigma^{2}_{u} \left[ \left( \dfrac{\partial f(x_t,\beta)}{\partial \beta} \right)^{´} \left( \dfrac{\partial f(x_t,\beta)}{\partial \beta} \right) \right]^{-1} \]
que se reduce a la matriz de covarianzas \(\sigma^{2}_{u}(X^{'}X)^-1\) en el caso de un modelo lineal. Si quisiéramos aplicar Minimos Cuadrados directamente, en el modelo exponencial,
\[ y_t= f(x_t,\theta)+u_t= \alpha + \beta_1e^{\beta_{2}x_t} + u_t \]
con \(\theta = (\alpha,\beta_1,\beta_2)\), tendríamos que resolver el problema,
\[ \min_{\theta} SR(\hat{\theta})= \min_{\theta}\sum_{t=1}^{T} \left[ \hat{u_t}(\hat{\theta}) \right]^2= \min_{\theta}\sum_{t=1}^{T}[y_t-(\alpha + \beta_1e^{\beta_{2}x_t})]^2 \]
que conduce a las condiciones de optimalidad,
\[ \begin{align*} \sum y_t &= \alpha T + \beta_1 \sum e^{\beta_2 x_t}\\ \sum y_te^{\beta_2 x_t} &= \alpha\sum e^{\beta_2 x_t}+\beta_1\sum e^{2\beta_2 x_t}\\ \sum y_tx_te^{\beta_2 x_t} &= \alpha\sum x_te^{2\beta_2 x_t}+\beta_1\sum x_te^{2\beta_2 x_t}\\ \end{align*} \]
que carece de solucion explícita, por lo que debe resolverse por procedimientos numéricos.

Para hacer esto en Rstudio, usaremos el comando nls (“Nonlineal Least Squares”), que nos agilizará este poceso. El comando precisa de una ecuación (el modelo exponencial anterior) y un gérmen inicial para cada parámetro. Estos siempre son muy difíciles de encontrar, pero por suerte, ya tenemos dos muy buenos, los valores de \(\alpha\) y \(\beta\) del anterior apartado.

\[A=\alpha = 1.035123; B=\beta = 0.1496007 \]

NLS<-nls(AllCasesUK~A*exp(B*days),start=list(A=a,B=b))
summary(NLS)

Formula: AllCasesUK ~ A * exp(B * days)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
A 1.381e+03  2.150e+02   6.426  5.9e-09 ***
B 5.429e-02  1.835e-03  29.580  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11320 on 91 degrees of freedom

Number of iterations to convergence: 24 
Achieved convergence tolerance: 8.401e-06

Como podemos observar en el summary, los parámetros han cambiado y se han realizado 24 iteraciones para lograrlo. Vamos a representar la gráfica con los nuevos parámetros para comprobar cuánto ha mejorado el ajuste:

A<-summary(NLS)$coefficients[1];A

[1] 1381.294

B<-summary(NLS)$coefficients[2];B

[1] 0.05429373

\[y=1381.294 e^{0.05429373}\]

xaj<-xaj<-seq(0.005,100,0.001)
yaj<-A*exp(B*xaj)
plot(days,AllCasesUK,ylab="Number of Cases",xlab="Days since 31-1-20", main="Total Cases UK")
lines(spline(xaj,yaj),col="green",lwd="4")

Cabe destacar que el cambio es muy grande. La nueva función se ajusta mucho mejor de forma general a los datos, dejando atrás la gran disparidad que mostraba el ajuste del apartado anterior.

Apartado d): modelar el número acumulado de personas infectadas mediante un ajuste no lineal de tipo asintótico más realista que los anteriores.

Para este apartado hemos decidido usar la “Función de Gompertz”. Esta función es un tipo de modelo matemático para una serie temporal y lleva el nombre de Benjamin Gompertz. Es una función sigmoidea que describe un crecimiento lento al comienzo y el crecimiento del final de un periodo de tiempo dado. Las funciones en base a este modelo, son inicialmente cóncavas y, tras pasar el punto de inflexión, convexas, conocidas también como curvas en forma de S. Benjamin Gompertz utilizó la curva para describir la ley de la naturaleza que rige la mortalidad humana.

Existen diferentes tipos de curvas Gompertz en función de los parámetros que la componen, pero todas tienen una doble exponencial como elemento característistico común.

Benjamin Gompertz introduce esta curva en su libro “On the Nature of the Function Expressive of the Law of Human Mortality, and on a New Mode of Determining the Value of Life Contingencies”. Para ello postula que los incrementos
\[Log(n)-Log(n+m),\ Log(n+m)-Log(n+2m),\ etc\]
son constantes, siendo \(Log(n+km)\) el logaritmo de la población en el instante \(n+km\), siendo indicio de encontrarnos ante un crecimiento geométrico.

La función que aparece en libro previamente citado, se expresa como:
\[L_x=dg^{q^{x}}, _{Ec1}\] donde

\[1\ \ log(g)=\dfrac{mq^{-a}}{1-q^{r}}\]

\[2\ \ q = p^{\dfrac{1}{r}}\]

\[3\ \ m=log(L_a)-log(L_{a+r})\]

\[4\ \ d=\dfrac{L_a}{\xi}\]

\[5\ \ log(\xi)=\dfrac{m}{1-q^{r}}\]

El valor de a denota el instante inicial, el de r la unidad de salto considerada en el tiempo y p es la razón de la progresión geométrica que define el número de personas vivas en el tiempo x.

La curva Gompertz fue objeto de interés en ciencias actuariales, pero con el paso del tiempo su aplicación se trasladó también a otros campos como la Biología o la Economía. P. Winsor, motivado por las diferentes aplicaciones de los modelos Gompertz, realizó en un estudio comparativo entre un modelo de Gompertz y otro de Verhulst. Para ello reescribió “Ec1” cómo:
\[Y = Ke^{-e^{a-bx}},\ k>0,\ b>0\]
deduciendose que \(lim_{x\rightarrow\infty}y=K\) y \(lim_{x\rightarrow-\infty}Y=0\). A partir de esta expresión, se obtiene la ecuación diferencial
\[\dfrac{\partial Y}{\partial x}=bye^{a-bx}\]
donde también se deduce que la pendiente de \(Y(x)\) es siempre positiva para valores finitos de \(x\), tendiendo a 0 cuando la \(x\) tiende a \(+\infty\). Buscando el punto de inflexión de la curva se llega a la expresión
\[\dfrac{\partial^2 y}{\partial x^2}= b^2ye^{a-bx}(e^{a-bx}-1)\]
de donde se extrae que está situado en \(\left(\dfrac{a}{b},\dfrac{k}{e}\right)\), alcanzando aproximadamente en el 36,78% del crecimiento total.

Con esto ya podemos calcular una aproximación de nuestro modelo para los datos que tenemos. Vamos a volver a plantearlos:

plot(days,AllCasesUK,ylab="Number of Cases",xlab="Days since 31-1-20", main="Total Cases UK")

Dado que el número de casos acumulados no presenta un comportamiento decreciente a simple vista, tomaré como punto de inflexión el último dato registrado.

AllCasesUK

 [1]      2      2      2      2      2      3      3      3      4      8      8      8      9      9      9      9      9      9
[19]      9      9      9      9      9     13     13     13     13     19     23     35     40     51     85    114    160    206
[37]    271    321    373    456    590    797   1061   1391   1543   1950   2626   3269   3983   5023   5688   6655   8085   9535
[55]  11665  14555  17115  19615  22285  25535  30105  34625  39295  43295  49495  53635  57525  63395  68075  73785  79015  84305
[73]  88645  93895  98495 103115 108715 114245 120095 124775 129075 133525 138105 143495 148405 152865 157175 161175 165255 171285
[91] 177485 182295 186635

Día 93, con un total de 186635 infectados. Lo siguiente será sacar unos valores de los párametros que necesitamos para:

\[Y = Ke^{-e^{a-bx}}\] Dado que, como hemos mencionado antes, en el punto de inflexión de esta función se cumple:

\[ \left. x=\dfrac{a}{b} \atop y=\dfrac{k}{e} \right\} \]
Cogeremos el punto (93,186635) para averiguar los valores:
\[ \left. x=\dfrac{a}{b} \atop y=\dfrac{k}{e} \right\} \left. 93 \approx\dfrac{a}{b} \atop 186635 \approx\dfrac{k}{e} \right\} \left. \begin{array}{rcl} a = 9.3\\ b = 0.1\\ k = 186635e \end{array} \right\} \]
Con estos datos ya podemos trabajar con el comando nls. Usaremos entonces los valores de a, b y k como gérmenes iniciales:

Gompertz<-nls(AllCasesUK~k*exp(-exp(a-b*days)),star=list(a=9.3,b=0.1,k=507326.5291),control=nls.control(minFactor=2^-24, maxiter=200))
summary(Gompertz)

Formula: AllCasesUK ~ k * exp(-exp(a - b * days))

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
a 4.177e+00  4.098e-02  101.95   <2e-16 ***
b 5.606e-02  7.123e-04   78.70   <2e-16 ***
k 2.609e+05  2.905e+03   89.79   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1022 on 90 degrees of freedom

Number of iterations to convergence: 159 
Achieved convergence tolerance: 2.967e-06

Al principio hemos tenido algún problema con los parámetros iniciales ya que en las iteraciones establecidas por defecto (unas 50), el método no lograba alcanzar la tolerancia establecida. Para ello hemos tenido que añadir un “nls.control” para solucionarlo. Como podemos observar en summary, el método converge tras 159 iteraciones, proporcionándonos una buena aproximación.

A continuación vamos a cargar los valores en unas variables para representar la gráfica.

a<-summary(NLS4)$coefficients[1];a

[1] 4.177245

b<-summary(NLS4)$coefficients[2];b

[1] 0.05605732

k<-summary(NLS4)$coefficients[3];k

[1] 260853.9

La ecuación de Gompertz para nuestro ajuste quedaría tal que así:
\[Y = 260853.9e^{-e^{4.177245-0.05605732x}}\]
Ahora vamos a dibujar el gráfico:

xaj<-seq(1,300,0.005)
yaj<-k*exp(-exp(a-b*xaj))
plot(days,AllCasesUK,xlim=c(0,150),ylim=c(0,300000),ylab="Number of Cases",xlab="Days since 31-1-20", main="Total Cases UK")
lines(spline(xaj,yaj),col="cyan",lwd="4")

Se puede observar que la función se ajusta perfectamente a los puntos e incluso nos predice de manera aproximada cuántos infectados va ha haber en El Reino Unido, \(Infectados\approx 250000\) (en el caso de que a partir de hoy empiecen a bajar el número de nuevos casos). De todas formas no podemos fiarnos de estos datos ya que hemos despreciado una gran cantidad de variables y solo nos hemos fijado en el número acumulado de infectados por día.


Como conclusiones podemos sacar que estos modelos son una gran ayuda a la hora de explicarnos como una epidemia en este caso va a desarrollarse, ya que nos ayudan con las predicciones. En cuanto a las gráficas, esta última se ajusta a una situación social de confinamiento, es decir, obligando a que la gente se quede en sus casas, obligamos también a frenar la expansión del virus. De esta forma la curva se aplana manteniendo constante el número global de casos y evitando una crísis sanitaria. Si estas medidas no se tomaran, la gráfica que mejor ajustaría la expansión del virus, sería el modelo exponencial que hemos hecho en el apartado b y luego mejorado en el apartado c. Una curva sin control y monotonamente creciente, que, hubiera colapsado el sistema.

Datos recogidos hasta el día 03/05/2020.

REFERENCIAS

https://www.google.com/covid19-map/?hl=es
https://www.overleaf.com/learn/latex/
https://masteres.ugr.es/moea/pages/curso201516/tfm1516/simon_minguez_tfm/!
https://es.wikipedia.org/wiki/Funci%C3%B3n_de_Gompertz
https://www.ucm.es/data/cont/media/www/pag-41459/Modelos%20no%20lineales%20corto.pdf
https://es.wikipedia.org/wiki/M%C3%ADnimos_cuadrados_no_lineales
https://www.ucm.es/data/cont/media/www/pag-41459/Modelos%20no%20lineales%20corto.pdf







Ejercicio 2.

Un profesor se interroga por la distribución que sigue el número de asignaturas aprobadas por los alumnos en el cuatrimestre que se imparte la asignatura antes citada de un total de cinco. Olvidando el hecho de que la proporción de aprobar una asignatura varía entre las cinco asignaturas, se aventura a afirmar que esas distribuciones siguen una B(5, 0. 85) en el grupo A y una B(5, 0. 1) en el grupo B. También sabe que a lo largo de los años, el 75 % de los alumnos asisten al grupo A y el 25 % asisten al grupo B. Si se escoge al azar un alumno que ha aprobado más de dos asignaturas en el cuatrimestre, ¿cuál es la probabilidad de que ese alumno haya asistido a clase en el grupo A?

Del enunciado obtenemos lar probabilidades de que un alumno pertenezca a cada uno de los grupos: \[P(A)=0.75\] \[P(B)=0.25\]

Pa<-0.75
Pb<-0.25

Del enunciado también extraemos que la porbabilidad de aprobar una asignatura en el grupo A es de 0.85; mientras que en el grupo B es de 0.1.

Tenemos que hallar la probabilidad de que un alumno que ha aprobado más de dos asignaturas haya asistido a clase en el grupo A. Esta se trata de una probabilidad condicionada. Como el alumno “ha aprobado más de dos asignaturas” establecemos que la variable aleatoria X es el número de asignaturas aprobadas y este tiene que ser estrictamente mayor que dos. Por tanto, tenemos que lo que nos pide el problema viene dado por: \[P(A/X>2)\] Para resolverlo aplicamos la regla de Bayes y posteriormente las propiedades de variables aleatiorias: \[P{(A/X>2)}=\frac{P(A)\cap P(X>2)}{P(X>2)}=\frac{(P(A)\cap (1-P(X\le2))}{1-P(X\le2)}\]

Primero hallamos la probabilidad de que un alumno haya aprobado más de dos asignaturas:

Pa2<-1-pbinom(2,5,0.85);Pa2

[1] 0.9733881

\[P{(X>2)}=(1-P(X\le2)=0.9733881\] A continuación, obtenemos la intersección de la probabilidad de que un alumno pertenezca al grupo A y la probabilidad de que un alumno haya aprobado dos asignaturas. Como se tratan de sucesos independientes, esta es la multiplicación de ambas.

Pa*Pa2

[1] 0.7300411

\[P(A)\cap (1-P(X\le2)=P(A)P(X>2)=0.7300411\]

Ahora, debemos obtener la probabilidad total de que un alumno haya aprobado más de dos asignaturas. Esto es lo mismo que la unión de la probabilidad de que un alumno haya aprobado más de dos asignaturas y vaya al grupo A y la probabilidad de que un alumno haya aprobado más de dos asignaturas y vaya al grupo B. \[P(X>2)=(P(A)\cap (1-P(X\le2))\cup(P(B)\cap(1-P(X\le2))\] Anteriormente, hemos hayado la probabilidad de que esto ocurra en el grupo A. Por lo que ahora debemos obtener la probabilidad de que esto ocurra en el grupo B.

Pb2<-1-pbinom(2,5,0.1);Pb2

[1] 0.00856

Como en el grupo A, se trata de sucesos independientes por lo que la intersección equivale al producto de probabilidadades:

Pb*Pb2

[1] 0.00214

\[P(B)\cap (1-P(X\le2)=P(B)P(X>2)=0.25*0.00856=0.00214\] Ahora ya podemos obtener la probabilidad total de aprobar más de dos asignaturas. Además, como la probabilizar de cumplir este requisito y pertenecer al grupo A es incompatible con la probabilidad de hacerlo y pertenecer al grupo B (ya que un alumno no puede pertenecer a ambos grupos a la vez) la unión de ambos es equivalente a la suma.

P2<-Pa*Pa2+Pb*Pb2;P2

[1] 0.7321811

\[P(X>2)=(P(A)\cap (1-P(X\le2))\cup(P(B)\cap(1-P(X\le2))=0.9733881∗0.75+0.00856∗0.25=0.7321811\] Finalmente, resolvemos:

Pa*Pa2/P2

[1] 0.9970772

\[P({A/X>2})=\frac{P(A)\cap P(X>2)}{P(X>2)}=\frac{(P(A)\cap (1-P(X\le2))}{1-P(X\le2)}=\frac{0.7300411}{0.7321811}=0.9970772\] Por tanto, la probabilidad de que habiendo aprobado más de dos asignaturas un alumno haya ido a clase al grupo A es 0.9970772.








EJERCICIO 3.

La duración de las bombillas de la Marca A siguen una distribución N(14000, 600), mientras que la duración de las bombillas de la Marca B siguen una distribución N(15000, 800). Se establece la siguiente regla para clasificar las bombillas producidas por las Marcas A y B: una bombilla será clasificada como Marca A si su duración es inferior al cuantil teórico 0.05 de la distribución que la gobierna la duración de las bombillas de la Marca B. Las bombillas de la Marca A se catalogan como deficientes si su duración es inferior al primer decil teórico de la distribución que gobierna su duración. Si se clasifica incorrectamente una bombilla de la Marca B, ¿cuál es la probabilidad de que sea clasificada como una bombilla deficiente de la Marca A?

Del enunciado extraemos que la variable aleatoria de duración de las bombillas sigue una distribución normal. Además, este tipo de distribución es característico de las variables aleatorias continuas ya que el conjunto de valores accesibles de duración de las bombillas es infinito no numerable. Tenemos que: -La variable aleatoria que gobierna duración de las bombillas de la Marca A sigue una distribución: \[X_a\sim N(14000,600)\]

-La variable aleatoria que gobierna duración de las bombillas de la Marca B sigue una distribución: \[X_b\sim N(15000,800)\] Así, hacemos las siguientes asignaciones:

mua<-14000
sigmaa<-600
mub<-15000
sigmab<-800

También sabemos que las bombillas se clasifican como Marca A si su duración se encuentra por debajo del cuantil teórico 0.05 de la distribución de las bombillas de la Marca B. Por tanto:

uc<-qnorm(0.05,mub,sigmab);uc

[1] 13684.12

En consecuencia, el umbral de duración de una bombilla por debajo del cual será clasificada como de la Marca A es 13684.12

Además, para que una bombilla sea clasificada como deficiente de la Marca A su duración ha de ser inferior al decil teórico (primer decil–> 0.1) de la distribución de bombillas de la Marca A.

ud<-qnorm(0.1,mua,sigmaa);ud

[1] 13231.07

Por tanto, se determina el umbral de duración de una bombilla tal que, por debajo del cual esta se clasifica como deficiente de la Marca A. Este umbral es 13231.07.

El ejercicio nos pregunta la probabilidad de que una bombilla sea clasificada como deficiente de la Marca A siendo esta una bombilla incorrectamente clasificada de la Marca B. Esto es una probabilidad condicionada: \[P(D/I)\] Siendo:

-** \(D\): Bombilla de la Marca B clasificada como deficiente de la Marca A**

-** \(I\): Bombilla de la Marca B incorrectamente clasificada como de la Marca A**

Para resolverlo aplicamos la regla de Bayes: \[P(D/I)=\frac{P(D)P(I/D)}{P(I)}\]

De acuerdo con el enunciado se tiene que la probabilidad de que una bombilla de la Marca B sea clasificada como de la marca A es de 0.05. \[P(I)=0.05\] La siguiente probabilidad que debemos calcular es P(I/D). Esta es la probabilidad de que una bombilla de la Marca B sea clasificada erróneamente como de la Marca A, habiendo sido previamente clasificada como deficiente de la Marca A. Esta es del 100 % ya que en el enunciado nos asegura que se clasifica incorrectamente una bombilla de la Marca B. Por tanto, tenemos: \[P(I/D)=1\]

Por último, debemos averiguar cual es la probabilidad de que una bombilla B sea clasificada como deficiente de la Marca A, es decir P(D). Antes, hemos calculado que para que esto ocurra su duración ha de estar por debajo del umbral 13231.07 (ud).

PD<-pnorm(ud,mub,sigmab);PD

[1] 0.01351225

Por lo tanto, tenemos que: \[P(D)=0.01351229\] Finalmente, operamos:

PD*1/0.05

[1] 0.270245

\[P(D/I)=\frac{P(D)P(I/D)}{P(I)}=\frac{0.01351229\cdot1}{0.05}=0.2702459\] En consecuencia, la probabilidad de que una bombilla de la Marca B que se ha clasificado incorrectamente como Marca A sea clasificada como deficiente de la marca A es 0.2702459.

0.2702459*100

[1] 27.02459

Esta expresada en porcentaje sería del 27.02459 %

Gráficamente:

curve(dnorm(x,mua,sigmaa),xlim=c(11000,20000),col="blue",lwd="2",ylab="Funcion densidad",xlab="Duracion de la bombilla",main="Clasificaci?n de bombillas")
curve(dnorm(x,mub,sigmab),xlim=c(6500,20000),col="red",lwd= "2",add=T)

lt2<-seq(11000,ud,0.01)
polygon(c(lt2[lt2<ud],ud),c(dnorm(lt2[lt2<ud],mua,sigmaa),0),col="grey")
abline(v=ud,col="black",lwd=1)

lt<-seq(11000,uc,0.01)
polygon(c(lt[lt<uc],uc),c(dnorm(lt[lt<uc],mub,sigmab),0),col="yellow")
abline(v=uc,col="black",lwd=1)

polygon(c(lt[lt<ud],ud),c(dnorm(lt[lt<ud],mub,sigmab),0),col="green")
legend(x="topright",legend=c("Distribucion Marca A","Distribucion Marca b","Clasificada como A","Clasificada como deficiente de A","P(D/I)"),fill=c("blue","red","yellow","grey","green"))