Inferencia Estadı́stica. Tutorı́a Grupal 4 (25-02-2016)
Índice
Dada una v.a. X con distribución Geométrica (p), se considera una m.a.s. de tamaño 15 y las hipótesis \(H_0 : p = 0,3\) frente \(H_1 : p \not= 0,3\).
1 Estimador máximo-verosímil
«Calcula el estimador de máxima verosimilitud para p.»
load(distrib) $ f(x,p) := pdf_geometric(x,p) $ /* densidad */ logexpand: all $ /* log(a*b^c) -> log(a)+c*log(b) */ logf(x,p) := log(f(x,p)) $ /* logverosimilitud: */ simpsum: true $ /* sum(cte+...,i,1,n) -> n*cte+... */ logV(X,p) := block ([n], if listp(X) then n:length(X), sum (logf(X[i],p), i, 1, n)) $ derV(X,p) := subst (p, q, diff (logV(X,q), q)) $ /* estimador máximo-verosímil */ EMVp(X) := rhs (solve (derV(X,p) = 0, p) [1]) $ EMVp(X) ; M : [2, 1, 8, 9, 0, 5, 0, 0, 3, 0, 0, 12, 0, 0, 1] $ EMVp(M) ; EMVp(M), numer ;
(%i10) EMVp(X);
n
(%o10) ------------
n
====
\
n + > X
/ i
====
i = 1
(%i11) EMVp(M) ;
15
(%o11) --
56
(%i12) EMVp(M), numer ;
rat: replaced -0.2678571428571428 by -15/56 = -0.2678571428571428
(%o12) 0.2678571428571428
2 Estadígrafo de razón de verosimilitudes.
«Calcula, analiza, representa el estadístico de la razón de verosimilitudes y exprésalo en función de un estadístico suficiente.»
Estadígrafos suficientes podrían ser:
- S = suma de la muestra = \(\sum_{i=1}^n X_i\)
- \(\bar X\) = media de la muestra = \(\frac 1 n \sum_{i=1}^n X_i\)
/* en función del tamaño muestral y de la suma */ kill (all) $ reset () $ p0 : 3/10 $ logV (n, S, p) := if p=1 then if S=0 then 1 else 0 else n * log(p) + S * log(1-p) $ EMV (n, S) := n / (n + S) $ logRV (n, S) := logV(n,S,p0) - logV(n,S,EMV(n,S)) $ derRV (n, S) := subst (S, s, diff (logRV(n,s), s)) $ der2RV (n, S) := subst (S, s, diff (logRV(n,s), s, 2)) $ der2RV (n, S) := subst (S, s, diff (derRV(n,s), s)) $ /* equivalente a la anterior */ M : [2, 1, 8, 9, 0, 5, 0, 0, 3, 0, 0, 12, 0, 0, 1] $ n : length (M) $ ratexpand (derRV (n, S)) ; /* ratexpand no necesario para solve */ solución : solve (derRV(n,S) = 0, S) ; /* [S = 35] */ der2RV (n, S), solución; /* -3/350 < 0 => máximo */ load (distrib) $ maxS (p) := quantile_negative_binomial (0.999999, n, p) $ x : makelist (i, i, 0, maxS(p0)) $ y : map (lambda ([s], logRV(n,s)), x) $ plot2d ([discrete, x, y]); /* S=35 maximiza la RV */ plot2d ([discrete, x, y], [gnuplot_term, dumb]) $ x [sublist_indices (y, lambda ([yi], yi = lmax(y))) [1]] ;
## en función del tamaño muestral y de la suma p0 <- 0.3 # p bajo H0 V <- function (n, S, p) p^n * (1-p)^S # verosimilitud EMV <- function (n, S) n / (n + S) # estimador máximo-verosímil RV <- function (n, S) V(n,S,p0) / V(n,S,EMV(n,S)) # razón de verosimilitudes ## muestra del enunciado M <- c (2, 1, 8, 9, 0, 5, 0, 0, 3, 0, 0, 12, 0, 0, 1) n <- length (M) # tamaño muestral = 15 maxS <- function (p) qnbinom (0.999999, n, p) x <- 0 : maxS(p0) # valores de S y <- RV (n, x) library (txtplot) # install.packages("txtplot") txtplot (x, y, xlab="S", ylab="RV") x [which.max (y)] # S=35 maximiza la RV
+--+------------+------------+------------+-------------+------------+---+
1 + ***** +
| * * |
| ** ** |
0.8 + * * +
| * * |
| * * |
0.6 + * * +
R | * * |
V | * * |
0.4 + * * +
| * * |
| * * |
0.2 + * ** +
| ** ** |
| *** **** |
0 + ********** ********************** +
+--+------------+------------+------------+-------------+------------+---+
0 20 40 60 80 100
S
3 Contraste de la razón de verosimilitudes.
«Calcula la región crítica del test de la razón de verosimilitudes al nivel de significación 0,05.»
S sigue una distribución binomial negativa BN(15,p).
RV tendrá una distribución discreta: contraste aleatorizado.
x <- 0 : maxS(p0) # valores de S y <- dnbinom (x, n, p0) # probabilidades de S txtplot (x, y, xlab="S", ylab="Pr") # distribución de S r <- RV (n, x) # valores de RV stopifnot (max (table (r)) == 1) # comprobar que son distintos txtplot (r, y, xlab="RV", ylab="Pr") # distribución de RV oo <- order (r, decreasing=TRUE) rr <- r[oo] xx <- x[oo] yy <- y[oo] cc <- cumsum (yy) # probabilidades acumuladas alfa <- 0.05 # nivel de significación ii <- which (cc > 1-alfa) [1] # aquí se sobrepasa el 95% cc[(ii-1):ii] # 0.9486096 0.9574614 rr[(ii-1):ii] # 0.1451617 0.1346742 RA <- xx [rr > rr[ii]] # valores de S seguro en región de aceptación ## con 18≤S≤60 se cubre Pr = 0,948.609.6 => faltan 0,001.390.4 p.falta <- 1-alfa - cc[ii-1] ## para la región de aceptación el siguiente sería S= x.alea <- xx[ii] # 17 p.alea <- yy[ii] # Pr[S=17] = 0,008.851.776 p.h0.alea <- p.falta / p.alea # 0,157.075.8 = Pr [elegir H0 | S=17]
Por tanto:
- H0 si \(18\leq S\leq 60\)
- H1 si \(S\leq16\) ó si \(S\geq61\)
- si \(S=17\) entonces H0 con Pr=0,157.075.8 y H1 con Pr=0,842.924.2
0.04 +--+------------+------------+------------+------------+------------+---+
| **** |
| ** ** |
| * * |
0.03 + ** ** +
| * * |
| ** |
| * * |
P 0.02 + * * +
r | * * |
| * * |
| * * |
0.01 + * ** +
| * ** |
| * *** |
| ** **** |
0 + ******** ************************* +
+--+------------+------------+------------+------------+------------+---+
0 20 40 60 80 100
S
0.04 +--+------------+------------+------------+-------------+------------+--+
| * * **** |
| * * ** |
| * * |
0.03 + * * * * +
| * * |
| * * |
| * * |
P 0.02 + * * +
r | * * * |
| * * |
| * * * |
0.01 + * * * +
| * * * * |
| * ** * * |
| *** * *** |
0 + ******* +
+--+------------+------------+------------+-------------+------------+--+
0 0.2 0.4 0.6 0.8 1
RV
4 P-valor.
«Calcula el p-valor asociado a la muestra (2, 1, 8, 9, 0, 5, 0, 0, 3, 0, 0, 12, 0, 0, 1).»
rM <- RV (n, sum (M)) # 0,868.682.3 sum(M)=41 sum (y [r <= rM]) # 0,610.848.4 = p-valor
5 Potencia
«Calcula la función potencia del test y represéntala gráficamente.»
RA. <- union (RA, x.alea) # región de aceptación incluyendo aleatorizado sop <- 0 : maxS(0.01) # soporte de S para menor valor de p representado potencia <- Vectorize (function (p) sum (dnbinom (setdiff(sop,RA.), n, p)) + (1 - p.h0.alea) * dnbinom (x.alea, n, p)) pes <- seq (0.01, 1, .01) txtplot (pes, potencia(pes), xlab="p", ylab="potencia")
+--+------------+------------+-------------+------------+-------------+--+
1 + ********* **************************** +
| ** *** |
| * ** |
0.8 + * ** +
p | * * |
o | ** |
t 0.6 + * * +
e | * * |
n | * |
c | * * |
i 0.4 + * * +
a | * ** |
| * * |
0.2 + * ** +
| ** ** |
| ***** |
+--+------------+------------+-------------+------------+-------------+--+
0 0.2 0.4 0.6 0.8 1
p
6 Unilateral
«Resuelve el mismo problema para las hipótesis H0 : p ≥ 0,3 frente H1 : p < 0,3.»
p0 <- 0.3 # menor valor de p bajo H0 => garantiza alfa EMV0 <- function (n, S) pmax (p0, EMV(n,S)) RV <- function (n, S) V(n,S,EMV0(n,S)) / V(n,S,EMV(n,S)) # razón de verosimilitudes x <- 0:100 # valores de S y <- RV (n, x) txtplot (x, y, xlab="S", ylab="RV") y <- dnbinom (x, n, p0) r <- RV (n, x) stopifnot (identical (r, rev(sort(r)))) # ya están ordenados decrecientemente cc <- cumsum (y) ii <- which (cc > 1-alfa) [1] # aquí se sobrepasa el 95% cc[(ii-1):ii] # 0.9443892 0.9520249 0.95 - 0.9443892 # 0.0056108 x[ii] # 54 0.0056108 / dnbinom (54, n, p0) # 0.7348188 = Pr (elegir H0 | S=54) ## contraste aleatorizado: ## - H0 si \(S\leq 53\) ## - H1 si \(S\geq 55\) ## - si \(S=54\) entonces H0 con Pr=0,734.818.8 y H1 con Pr=0,265.181.2 sum (y [(which(x==sum(M))+1) : length(y)]) # 0.2548661 p-valor potencia <- function (p) 1-pnbinom(54,n,p) + 0.2651812*dnbinom(54,n,p) pes <- seq (0, 1, .01) txtplot (pes, potencia(pes), xlab="p", ylab="potencia")
+--+------------+------------+------------+-------------+------------+---+
1 + ************************** +
| * |
| ** |
0.8 + * +
| * |
| * |
0.6 + * +
R | * |
V | * |
0.4 + * +
| * |
| * |
0.2 + ** +
| ** |
| **** |
0 + ********************** +
+--+------------+------------+------------+-------------+------------+---+
0 20 40 60 80 100
S
+--+------------+------------+-------------+------------+-------------+--+
1 + ********* +
| ** |
| * |
0.8 + * +
p | * |
o | * |
t 0.6 + * +
e | * |
n | * |
c 0.4 + +
i | * |
a | * |
0.2 + ** +
| * |
| *** |
0 + ************************************************ +
+--+------------+------------+-------------+------------+-------------+--+
0 0.2 0.4 0.6 0.8 1
p