4pl.html

Sea X una variable aleatoria con distribución binomial B(n;p) de la que se sabe que n ∈ [5;15] y p ∈ [0,2;0,4] y de la que se dispone de una muestra de tamaño 500. Calcular la EMV de (n;p). En este ejercicio se trata de ver la importancia de trabajar con el log de la verosimilitud para evitar los problemas de ”verosimilitud plana”.

Crea las siguentes variables en R:
- N que valga 500
- nmin que valga 5
- nmax que varga 15
- pmin que valga 0,2
- pmax que valga 0,4
```
N <- 500
nmin <- 5;   nmax <- 15
pmin <- 0.2; pmax <- 0.4
```
Crea la variable sop.n (soporte de n) que contenga el vector de posibles valores que puede tomar el parámetro n.
```
sop.n <- nmin:nmax
```
Genera aleatoriamente el valor de p y no lo mires.
```
p <- runif(1, pmin, pmax)
```

Genera aleatoriamente el valor de n y no lo mires.

n <- sample(sop.n, 1) 
## otra posibilidad:
n <- floor(runif(1, nmin, nmax+1)
## según ?runif
## ‘runif’ will not generate either of the extreme values

Genera la muestra de tamaño 500 de la B(n;p) para el estudio.
```
X <- rbinom(N, n, p)
```
Calcula teóricamente la función de verosimilitud de la muestra.
\[ L(n,p) = \mathcal L(x_1,...,x_N\mid n, p) = \prod_{i=1}^N f_{B(n;p)}(x_i) = \prod_{i=1}^N {n \choose x_i} p^{x_i} (1-p)^{n-x_i} \propto \prod_{i=1}^N p^{x_i} (1-p)^{n-x_i} = p^{\sum_{i=1}^N x_i} (1-p)^{N n - \sum_{i=1}^N x_i} \]
Implementa en R la función de verosimilitud L en función de dos argumentos n y p.
```
L <- function (n, p) prod (dbinom (X, n, p))
```

Calcula la verosimilitud para varias combinaciones de los parámetros.

summary(X)  # el máximo da una cota inferior de "n"
L (10, .3)
L (12, .2)
L (12, .4)

Implementa en R la función de logverosimilitud logL en función de dos argumentos n y p.
```
logL <- function (n, p) sum (dbinom (X, n, p, log=TRUE))
```

Calcula la logverosimilitud para varias combinaciones de los parámetros.

summary(X)  # el máximo da una cota inferior de "n"
logL (10, .3)
logL (12, .2)
logL (12, .4)

Para una n fija, calcula teóricamente la estimación máximo-verosímil de p.
\[ L(n) = \max_p L(n,p) = L(n,\hat p_\text{MV}) \]

\[ \hat p_\text{MV} = \arg\max_p L(n,p) = \arg\max_p p^{\sum_{i=1}^N x_i} (1-p)^{N n - \sum_{i=1}^N x_i} = \]

\[ = \arg\max_p \log L(n,p) = \arg\max_p {\sum_{i=1}^N x_i}\log p + (N n - \sum_{i=1}^N x_i)\log(1-p) \]

\[ 0 = {\partial \log L(n,p) \over \partial p} = \frac {\sum_{i=1}^N x_i} {p} - \frac {N n - \sum_{i=1}^N x_i} {1-p} \Longrightarrow \hat p_\text{MV} = \frac {\sum_{i=1}^N x_i} {N n} \]

Por otro lado, hay que tener en cuenta el espacio paramétrico \(0.2 \leq p \leq 0.4\).

Implementa en R una función p.mv que, dado un argumento n, devuelva la estimación máximo-verosímil de p.

p.mv <- function (n) pmax(pmin, pmin(pmax, sum(X) / N / n)) # "pmin" es un escalar, "pmin(" es una función; ídem "pmax"

Implementa en R la función logLn que dependa de un parámetro n y que devuelva la máxima logverosimilitud para tal n.
```
logLn <- function (n) logL (n, p.mv(n))
```

Dibuja en función de n la función implementada en el apartado anterior.

y <- sapply(sop.n, logLn)
plot (sop.n, y, type="b") # los -Inf no se ven

Calcula la estimación máximo-verosímil de n.

imax <- which.max (y)
n.mv <- sop.n[imax]   # EMV de "n"

Calcula la estimación máximo-verosímil de p.
```
p.mv (n.mv)
```

Calcula las estimaciones de n y p mediante el método de los momentos.

## mu = n.p; S = n.p.(1-p)  =>  p = (mu-S)/mu; n = mu^2/(mu-S)

mu   <- mean(X)
S    <-  var(X) * (N-1)/N
p.mm <- (mu-S) / mu
n.mm <- mu^2 / (mu-S)