Ejercicio 1. Modelo de mixturas

En 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\)).

  1. Obtener el EMV para los parámetros e IC de Wald con confianza aproximada de 95% para \(p\) y para \(l\).
# Datos observados
 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
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) # valores coherentes con p y lambda, respectivamente
minimo<-optim(vi,mlv,hessian=TRUE)
## Warning in log(1 - p): Se han producido NaNs
emv<-minimo$par
emv
## [1] 0.4698821 4.6207927
infObs<-minimo$hessian
emvVar<-solve(infObs)
emvVar
##             [,1]        [,2]
## [1,] 0.006365515 0.001203439
## [2,] 0.001203439 0.228362471
# Errores estandar e IC de Wald para los dos parámetros
se<-sqrt(diag(emvVar))
tabla.Wald<-cbind(emv,se,extremo.inferior=emv-qnorm(0.975)*se,
    extremo.superior=emv+qnorm(0.975)*se)
nombres.par<-c("p","l")
rownames(tabla.Wald)<-nombres.par
round(tabla.Wald,4) #usar 4 decimales
##      emv     se extremo.inferior extremo.superior
## p 0.4699 0.0798           0.3135           0.6263
## l 4.6208 0.4779           3.6842           5.5574
  1. Obtener una región de confianza simultánea para (p,l) al 95%, basada en la razón de verosimilitud.
lv<-function(p,lambda){
  return(19*log(p+(1-p)*exp(-lambda))+21*log(1-p)-21*lambda+sum(y)*log(lambda))
}

d2<--minimo$value-qchisq(0.95,2)/2
pSeq<-seq(0,1,length.out=150) # conocemos EMV, para definir mejor este grid
lSeq<-seq(1,8,length.out=150)
#g<-outer(pSeq,lSeq,lv)
g<-outer(pSeq,lSeq,lv)
contour(pSeq,lSeq,g,levels=c(d2),ylim=c(3,6))

emv
## [1] 0.4698821 4.6207927
  1. Obtener IC al 95% para p y l basados en el test de la razón de verosimilitudes.
# Manera 1
# Similar a ejecicio 6 hecho en clase. Expresión compleja
#h<-function(p){
  #return()
#}
#d1<--minimo$value-qchisq(0.95,1)/2
#pSeq<-seq(0,1,length.out=150)
#plot(pSeq,h(pSeq),type="l")
#abline(h=d1)


# Manera 2
d1<--minimo$value-qchisq(0.95,1)/2
pSeq<-seq(0,1,length.out=150)
lSeq<-seq(1,8,length.out=150)
o<-outer(pSeq,lSeq,lv)
contour(pSeq,lSeq,o,levels=c(d1),ylim=c(3,6))# está contenida en la anterior
abline(v=c(0.317,0.625),col=2) # Comparar con los IC calculados para Wald
abline(h=c(3.75,5.645),col=2)

  1. Obtener el p-valor del test basado en el estadístico RV para contrastar la hipótesis \(H_0: \lambda=2p^2+4\) contra \(H_1: \lambda\neq 2p^2+4\).
mlvH0<-function(p){  
  lambda<-2*p^2+4
  return(-(19*log(p+(1-p)*exp(-lambda))+21*log(1-p)-21*lambda+sum(y)*log(lambda)))
}
minimoH0<-optim(par=0.5,fn=mlvH0,method="Brent",lower=0.01,upper=0.99,hessian=TRUE)# Uniparamétrico usamos Brent
emvH0<-minimoH0$par
emvH0
## [1] 0.4780091
Qobs<-2*(-minimo$value+minimoH0$value)
Qobs
## [1] 0.1332147
pvalRV<-1-pchisq(Qobs,1)
pvalRV
## [1] 0.7151219
  1. Obtener el p-valor del test basado en el estadístico de Wald para contrastar la hipótesis \(H_0: \lambda=2p^2+4\) contra \(H_1: \lambda\neq 2p^2+4\). (Repaso)
g<-function(theta){
  p<-theta[1]
  lambda<-theta[2]
  return(lambda-2*p^2-4)
}
G<-cbind(-4*emv[1],1)
uve<-G%*%emvVar%*%t(G)
Wobs<-(g(emv)-0)^2/uve
pvalW<-1-pchisq(Wobs,1)
pvalW
##           [,1]
## [1,] 0.7180305

Ejercicio 2. Modelo Beta

En una muestra aleatoria de tamaño 30 de una población \(B(a,b)\) de parámetros \(a=b=2\):
(https://mathworld.wolfram.com/BetaDistribution.html https://www.wolframalpha.com/input/?i=beta+distribution)

  1. Obtener EMV de a y b, su distribución asintótica e IC Wals con confienza aproximada 95% para a y b.
n<-30
set.seed(1234)
x<-rbeta(n,2,2) # Datos

A<-sum(log(x))
B<-sum(log(1-x))

mlv<-function(tita){
    a<-tita[1]
    b<-tita[2]
    return(30*lbeta(a,b)-A*(a-1)-B*(b-1)) #lbeta(a,b) calcula log(beta(a,b))
}


vi<-c(1,1) # ya que en la muetra proviene de B(2,2)
minimo<-optim(vi,mlv,hessian=TRUE)
emv<-minimo$par
emv
## [1] 2.25584 3.18622
infobs<-minimo$hessian
emvVar<-solve(infobs)
emvVar
##           [,1]      [,2]
## [1,] 0.3040330 0.3682515
## [2,] 0.3682515 0.6462389
# Errores estandar e IC de Wald para los dos parámetros
se<-sqrt(diag(emvVar))
tabla.Wald<-cbind(emv,se,extremo.inferior=emv-qnorm(0.975)*se,
    extremo.superior=emv+qnorm(0.975)*se)
nombres.par<-c("a","b")
rownames(tabla.Wald)<-nombres.par
round(tabla.Wald,4) #usar 4 decimales
##      emv     se extremo.inferior extremo.superior
## a 2.2558 0.5514           1.1751           3.3365
## b 3.1862 0.8039           1.6106           4.7618
  1. Obtener el p-valor basado en el estadístico de RV para contrastar la hipótesis (compuesta) \(H_0:a=2\) vs \(H_1:a\neq 2\)
mlvHO<-function(b){
    return(30*lbeta(2,b)-A*(2-1)-B*(b-1))
}
minimoH0<-optim(1,mlvHO,method="Brent",hessian=TRUE,lower=0.1,upper=4)
emvH0<-minimoH0$par
emvH0
## [1] 2.876028
T1<- 2*(-minimo$value+minimoH0$value)
T1
## [1] 0.2340296
pvalorT1<-1-pchisq(T1,1)
pvalorT1
## [1] 0.6285519
  1. Obtener el p-valor basado en el estadístico de RV para contrastar la hipótesis (compuesta) \(H_0:a=b\) vs \(H_1:a\neq b\)
mlvH2<-function(b){
    return(30*lbeta(b,b)-A*(b-1)-B*(b-1))
}
minimoH02<-optim(1,mlvH2,method="Brent",hessian=TRUE,lower=0.1,upper=4)
emvH02<-minimoH02$par
emvH02
## [1] 2.277675
T2<- 2*(-minimo$value+minimoH02$value)
T2
## [1] 5.269179
pvalorT2<-1-pchisq(T2,1)#cuando a=b tenemos un solo parametro
pvalorT2
## [1] 0.02170625
  1. Obtener la region de confianza simultane para (a,b) con confianza aproximada 0.95 basada en el test de la razón de verosimilitudes.
d2<--minimo$value-qchisq(0.95,2)/2
aSeq<-seq(0,6,length.out=100) # Definir una secuencia con sentido para los parametros
bSeq<-seq(0,6,length.out=100)
o1<-outer(aSeq,bSeq,function(a,b)-30*lbeta(a,b)+A*(a-1)+B*(b-1)) # o definir lv
contour(aSeq,bSeq,o1,levels=c(d2))

#Comprobamos que esta EMV
emv
## [1] 2.25584 3.18622
  1. Obtener intervalos de confianza para a y b con confianza aproximada 0.95. para a y b
d1<--minimo$value-qchisq(0.95,1)/2
aSeq<-seq(0,6,length.out=100)# Definir una secuencia con sentido para los parametros
bSeq<-seq(0,6,length.out=100)
#definimos la cantidad f
o2<-outer(aSeq,bSeq,function(a,b)-30*lbeta(a,b)+A*(a-1)+B*(b-1)) # o definir lv
contour(aSeq,bSeq,o1,levels=c(d2))
contour(aSeq,bSeq,o2,levels=c(d1),add=TRUE)

# IC para a 
abline(v=1.35,col=2)#extremo inferior 
abline(v=3.55,col=2)#extremo superior 
# IC para a 
abline(h=1.85,col=2)#extremo inferior
abline(h=5,col=2)#extremo superior

# Se pueden calcular con los intervalos para Wald
  1. Obtener el p-valor basado en el estadístico de Wald para contrastar la hipótesis \(H_0:a=2\) vs \(H_1:a\neq 2\)
#Se trata del valor de a con el que ha generado la muestra. No debería rechazarse
W1<-(emv[1]-2)^2/emvVar[1,1]
W1
## [1] 0.2152854
pvalorW1<-1-pchisq(W1,1)
pvalorW1
## [1] 0.6426559
  1. Obtener el p-valor basado en el estadístico de Wald para contrastar la hipótesis \(H_0:a=b\) vs \(H_1:a\neq b\)
g<-function(theta){
  a<-theta[1]
  b<-theta[2]
  return(a-b)
}
G<-cbind(1,-1)
uve<-G%*%emvVar%*%t(G)
W2<-((g(emv)-0)^2)/uve
W2
##          [,1]
## [1,] 4.049267
pvalorW2<-1-pchisq(W2,1)
pvalorW2
##            [,1]
## [1,] 0.04419053

Ejercicio propuesto

Sean \(X_1\),,\(X_{100}\) v.a.i.i.d. distribución \(N(\mu,\sigma)\). Los datos observados de los que se dispone son: \(\sum_{i=i}^n X_i=134.94\) y \(\sum_{i=1}^n X_i^2 = 547.43\).

  1. Obtener el p-valor del test de Wald para contrastar \(H_0: \mu=\sigma\) contra \(H_1: \mu\neq\sigma\).
  2. Obtener el p-valor del test de RV para contrastar la hipótesis nula \(H_0: \mu=\sigma\) contra \(H_1: \mu\neq\sigma\).
  3. Obtener la región de confianza simultánea para \((\mu,\sigma)\), con confianza aproximada 0.95.
  4. Obtener intervalos de confianza de Wald, con confianza aproximada 0.95, para ambos parámetros \(\mu\) y \(\sigma\).
  5. Obtener intervalos de confianza basados en la RV, con confianza 0.95, para ambos parámetros \(\mu\) y \(\sigma\).