Sean \(Y_1,\dots, Y_{200}\) v.a.i.i.d. de un modelo de mixtura de dos normales conocidas \(N(0,1)\) y \(N(1,0.8)\) con densidad \(g(y,p)=p*dnorm(y)+(1-p)*dnorm(y,1,0.8),\) es decir, \(g\) es la función de densidad de los datos observados (incompletos). A partir de a muestra dada, obtenga el EMV para \(p\):
# Definir la muestra y de tamaño n=200
set.seed(1234)
x<-rnorm(200) # m.a.s. de tamaño 200 de una N(0,1)
y<-rnorm(200,1,0.8) # m.a.s. de tamaño 200 de una N(1,0.8)
r<-rbinom(200,1,0.3) # m.a.s. de tamaño 200 de una B(0.3)
y[r>0]<-x[r>0] # Reemplazar en y con los valores de x, en aquella posiciones donde r sea 1
# Aunque hemos simulado con p=0.3, trabajamos como desconocido# Menos logverosimilitud de los datos observados (incompletos)
mlv<-function(p,z) { #Depende de p y z, en optim especificaremos los valores de z
verosimilitud<-p*dnorm(z,0,1)+(1-p)*dnorm(z,1,0.8)
return(-sum(log(verosimilitud)))
}
vi<-c(0.4) #valor inicial para el calculo del emv
minimo<-optim(vi,mlv,z=y,method="Brent",lower=0.001,upper=0.999,hessian=T)#se especifica z
EMV_optim<-minimo$par
EMV_optim## [1] 0.2854986# Primera derivada de la logverosimilitud de los datos observados (incompletos)
lprima<-function(p){
num<-dnorm(y)-dnorm(y,1,0.8)
dnom<-p*dnorm(y)+(1-p)*dnorm(y,1,0.8)
return(sum(num/dnom))
}
# Segunda derivada de la logverosimilitud de los datos observados (incompletos)
H<-function(p){
num<--(dnorm(y)-dnorm(y,1,0.8))^2
dnom<-(p*dnorm(y)+(1-p)*dnorm(y,1,0.8))^2
return(sum(num/dnom))
}
# Fórmulas de actualización: theta^(k+1)=theta^(k)-H(theta^(k),X)^-1*lprima(theta^(k),X)
iter<-10 # Basta con iter<-6
p0<-0.2
p<-c(p0,rep(0,iter))
for(k in 2:(iter+1)){
p[k]<-p[k-1]-(lprima(p[k-1])/H(p[k-1]))
}
EMV_NR<-p
EMV_NR## [1] 0.2000000 0.2706448 0.2851369 0.2854984 0.2854986 0.2854986 0.2854986
## [8] 0.2854986 0.2854986 0.2854986 0.2854986# Fórmulas de actualización: theta^(k+1)=1/n*sum((p^(k)*dnorm(y))/(p^(k)*dnorm(y)+(1-p^(k))*dnorm(y,1,0.8)))
iter<-10 # Aumentar 10, 25, 52
p0<-0.2
p<-c(p0,rep(0,iter))
for(k in 2:(iter+1)){
num<-p[k-1]*dnorm(y)
deno<-p[k-1]*dnorm(y)+(1-p[k-1])*dnorm(y,1,0.8)
p[k]<-mean(num/deno)
}
EMV_EM<-p
EMV_EM## [1] 0.2000000 0.2200883 0.2358240 0.2479619 0.2572308 0.2642614 0.2695700
## [8] 0.2735658 0.2765669 0.2788173 0.2805030Se tienen 100 individuos que se clasifican en una de las 6 celdas de una tabla de contingencia \(2\times 3\), en la que se sabe que las filas y columnas son independientes, pero para la que se ha perdido la información de 35 individuos correspondientes a las celdas (1,3) y (2,2). Los datos observados fueron: La tabla observada es:
| 10 | 5 | ? |
| 20 | ? | 30 |Se cuentan con \(X_1,\dots X_{100}\) v.a.i.i.d. de una distribución \(P_\theta\) tal que \(Y_{z}=\sum_{i=1}^{100} \mathbb{I}(X_i=z), z=1,\dots,6\) con probabilidad \(p_1,\dots,p_6,\) y tal que \(\sum_{z=1}^6 p_z=1\), de donde se sigue la siguiente tabla de probabilidad:
| p_1 | p_2 | p_3 | r_1 |
| p_4 | p_5 | p_6 | 1-r_1 |
| c_1 | c_2 | 1-c_1-c_2 | |Además se sabe que no se dispone de datos para \(Y_3\) e \(Y_5\) y que las filas y columnas son independientes, luego \((p_1,\dots,p_6)=(r_1c_1,r_1c_2,\dots,(1-r_1)(1-c_1-c_2)),\) y por tanto \(\theta=(r_1,c_1,c_2).\)
n<-100
y<-c(10,5,2,3,35) # (y1,y2,y4,y6,100-(y1+y2+y4+y6))
w<-c(15,30,5) # (r1,c1,c2)
mlv<-function(tita) {
r1<-tita[1];c1<-tita[2];c2<-tita[3]
return(-sum(w*log(tita))-50*log(1-r1)-30*log(1-c1-c2)-35*log(r1*(1-c1-c2)+(1-r1)*c2)) # sum(y*log(tita)) simplifica
}
vi<-c(1/3,1/3,1/6) # válidos
minimo<-optim(vi,mlv,hessian=T)
EMV_optim<-minimo$par
EMV_optim## [1] 0.3748990 0.2999751 0.1751209# Hemos calculado la probabilidad:
#P(Y3/Y3+Y5=35)=P(Y3=k,Y5=35-k)/P(Y3+Y5=35) = (35!/k!*(35-k)!)*(P3/ P3+P5)^k) * (1-p3/ P3+P5)^(35-k)
n<-100
y1<-10;y2<-5;y4<-20;y6<-30 #datos observados
r1<-1/3; c1<-1/3; c2<-1/6; r2<-1-r1;c3<-1-c1-c2 #valores iniciales de los parámetros
iter<-20
for(k in 1:iter){
#paso E
# Podemos hacerlo como en clase
# Paso E:
y31Star<-35*(r1*(1-c1-c2))/((r1*(1-c1-c2))+(1-r1)*c2)
y51Star<-35-y31Star
# O trabajar en función de p que es más simple
#paso E
#p3<-r1*c3; p5<-r2*c2
#y31Star<-35*p3/(p3+p5); y51Star<-35-y31Star
#paso M
r1<-(15+y31Star)/n
c1<-0.3
c2<-(40-y31Star)/n
# O trabajar en función de p que es más simple
#paso M
#r1<-(y1+y2+y31Star)/n; r2<-1-r1
#c1<-(y1+y4)/n; c2<-(y2+y51Star)/n; c3<-1-c1-c2
cat("\n", "Iter",k,": r1=",r1,",c1=",c1,", c2=",c2,sep="")
}##
## Iter1: r1=0.36,c1=0.3, c2=0.19
## Iter2: r1=0.3605505,c1=0.3, c2=0.1894495
## Iter3: r1=0.3610844,c1=0.3, c2=0.1889156
## Iter4: r1=0.361602,c1=0.3, c2=0.188398
## Iter5: r1=0.3621036,c1=0.3, c2=0.1878964
## Iter6: r1=0.3625895,c1=0.3, c2=0.1874105
## Iter7: r1=0.3630601,c1=0.3, c2=0.1869399
## Iter8: r1=0.3635155,c1=0.3, c2=0.1864845
## Iter9: r1=0.3639562,c1=0.3, c2=0.1860438
## Iter10: r1=0.3643824,c1=0.3, c2=0.1856176
## Iter11: r1=0.3647944,c1=0.3, c2=0.1852056
## Iter12: r1=0.3651926,c1=0.3, c2=0.1848074
## Iter13: r1=0.3655773,c1=0.3, c2=0.1844227
## Iter14: r1=0.3659487,c1=0.3, c2=0.1840513
## Iter15: r1=0.3663073,c1=0.3, c2=0.1836927
## Iter16: r1=0.3666534,c1=0.3, c2=0.1833466
## Iter17: r1=0.3669872,c1=0.3, c2=0.1830128
## Iter18: r1=0.3673091,c1=0.3, c2=0.1826909
## Iter19: r1=0.3676194,c1=0.3, c2=0.1823806
## Iter20: r1=0.3679184,c1=0.3, c2=0.1820816En un estudio de micropropagación de raices, se anotó el número de raices de 40 variedades de manzana. Los datos observados fueron:
| Nº raices | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
|-------------|----|---|---|---|---|---|---|---|---|---|
| Frecuencia | 19 | 2 | 2 | 4 | 3 | 1 | 4 | 3 | 0 | 2 |Se ha establecido como modelo para describir este fenémeno una mixtura de dos distribuciones \(f_1\) y \(f_2\) con proporciones \(p\) y \(1-p\) respectivamente. En concreto, \(f_1\) toma el valor cero con probabilidad 1 (\(f_1(0)=1\)), y \(f_2\) es la función de masa de probabilidad de una distribución de Poisson(\(\lambda\)). Obtener el EMV para los parámetros:
y<-c(rep(0,19),rep(1,2),rep(2,2),rep(3,4),rep(4,3),rep(5,1),rep(6,4),rep(7,3),rep(9,2))
mlv<-function(theta){
p<-theta[1]
lambda<-theta[2]
return(-(19*log(p+(1-p)*exp(-lambda))+21*log(1-p)-21*lambda+sum(y)*log(lambda)))#ojo sum(y), empieza en 0
}
vi<-c(0.5,1)
minimo<-optim(vi,mlv,hessian=TRUE)
EMV_optim<-minimo$par
EMV_optim## [1] 0.4698821 4.6207927n<-40
iter<-8
p0<-0.5
l0<-1
p<-c(p0,rep(0,iter))
l<-c(l0,rep(0,iter))
for(i in 2:(iter+1)){
p[i]<-(19/n)*(p[i-1]/(p[i-1]+(1-p[i-1])*exp(-l[i-1]))) # Para el resto de frecuencias b_i^* es cero
l[i]<-mean(y)*(1/(1-p[i]))
cat("\n", "Iter",i,": p=",p[i],",l=",l[i],sep="")
}##
## Iter2: p=0.3472528,l=3.753367
## Iter3: p=0.4549552,l=4.495044
## Iter4: p=0.4687308,l=4.611598
## Iter5: p=0.4697103,l=4.620116
## Iter6: p=0.4697751,l=4.62068
## Iter7: p=0.4697793,l=4.620717
## Iter8: p=0.4697796,l=4.62072
## Iter9: p=0.4697796,l=4.62072EMV_EM<-c(p[iter+1],l[iter+1])
EMV_EM## [1] 0.4697796 4.6207200Muestras de suelo, seleccionadas aleatoriamente, se combinan en lotes y se examina si en los lotes se detecta o no la presencia de una toxina. Por tanto, las muestras de suelo individuales no se examinan. De \(100\) lotes, de \(10\) muestras de suelo cada uno, se detectó la toxina en \(12\) lotes. Calcular el EMV para \(p=\text{{probabilidad de que la toxina esté presente en una muestra de suelo individual}}\)
#EMV obtenido de forma analítica es EMV_ana<-0.0127 pHat=1-(1-thetaHat)^0.1
# EMV con optim
sY<-12
mlv<-function(theta){
return(-sY*log(theta)-(100-sY)*log(1-theta))
}
vi<-0.1
minimo<-optim(0.1,mlv,method="Brent",lower=0.01,upper=0.99,hessian=TRUE)
EMV_optim<-minimo$par
EMV_optim## [1] 0.12iter<-8
p0<-0.1
p<-c(p0,rep(0,iter))
for(i in 2:(iter+1)){
p[i]<-3*p[i-1]/(25*(1-(1-p[i-1])^10))
}
EMV_EM<-p
EMV_EM## [1] 0.10000000 0.01842408 0.01302880 0.01272047 0.01270302 0.01270204 0.01270198
## [8] 0.01270198 0.01270198Se sabe que la duración de unos componentes eléctricos sigue una distribución exponencial negativa con función de densidad \(f(x,\lambda)=\lambda e^{-\lambda x}, x>0.\) Con el propósito de estudiar la duración media de los componentes se ponen a funcionar \(100\) de ellos y se observa que la duración de todos los componentes fue mayor de \(3\) y menor de \(10\) unidades de tiempo. Obtener un EMV para la duración media de los componentes (\(1/\lambda\)).
# El EMV calculado de forma analitica es
emv<-log(10/3)/7
1/emv## [1] 5.814085#EMV obtenido de forma analítica es EMV_ana<-5.814
mlv<-function(lambda) {
return(-100*log(exp(-3*lambda)-exp(-10*lambda)))
}
vi<-5
minimo<-optim(vi,mlv,hessian=T)
EMV<-minimo$par
EMV_optim<-1/EMV
EMV_optim # Menos preciso## [1] 5.818182iter<-10
l0<-5
l<-c(l0,rep(l0,iter))
for(i in 2:(iter+1)){
num<-exp(-3*l[i-1])*(3+1/l[i-1])-exp(-10*l[i-1])*(10+1/l[i-1])
deno<-exp(-3*l[i-1])-exp(-10*l[i-1])
B<-num/deno
l[i]<-1/B
}
EMV_EM<-1/l
EMV_EM## [1] 0.200000 3.200000 5.315369 5.753114 5.807157 5.813304 5.813997 5.814075
## [9] 5.814084 5.814085 5.814085Se ponen en funcionamiento \(100\) y se anota su tiempo de fallo hasta un instante de tiempo \(t=\text{1 minuto}.\) Se sabe que el tiempo de fallo de estos aparatos se modela según un modelo exponencial negativo con función de densidad \(f(x,\lambda)=\lambda e^{-\lambda x}, x>0.\) Concretamente, se anotaron los siguiente tiempos de fallo:
tFallo<-c(0.59683969, 0.34819508, 0.63876960, 0.03998816, 0.54655953, 0.70678049,
0.43481007, 0.75293272, 0.36484336, 0.30681922)Obtener un EMV para la duración media de los componentes (\(1/\lambda\)):
A partir de los datos observados, tanto analíticamente como de forma aproximada.
A partir de los datos completos, utilizando el algoritmo EM con valor inicial 2 y 8 iteraciones.
Se ponen en funcionamiento dos aparatos cuyos tiempos de fallo se modelan según una exponencial negativa de parámetro \(\lambda.\) Sin embargo, únicamente se dispone de la información del primer aparato que falla. Es decir se tiene que \(X_1\), \(X_2\) i.i.d \(\sim exp(\lambda)\) son los datos completos e \(Y=\min(X_1,X_2)\) los datos observados (incompletos). Calcular el EMV para \(\lambda\):
A partir de los datos observados, tanto analíticamente como de forma aproximada.
A partir de los datos completos, utilizando el algoritmo EM. En el paso E, se requiere utilizar el teorema de cambio de variable multidimensional. Solución \(\hat{\lambda}=0.25\).