## 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”.

N <- 500
nmin <- 5;   nmax <- 15
pmin <- 0.2; pmax <- 0.4

## Genera aleatoriamente n y p en las condiciones anteriores y no los
## mires.

sop.n <- nmin:nmax # soporte de n
n <- sample(sop.n, 1)
p <- runif(1, pmin, pmax)

## A continuación genera la muestra de tamaño 500 de la B(n;p)
## para el estudio.

m <- rbinom(N, n, p)

## 2. Calcula y dibuja, como función de n, el máximo en
## p de la función de verosimilitud, max {L(m; n, p) | p},
## y compárala con la de su logaritmo.

emv.p <- structure (pmax(pmin, pmin(pmax, sum(m)/sop.n/N, names=sop.n))) # estima de p en función de n

Ln    <- function (n) prod (dbinom (m, n, emv.p[as.character(n)]))
logLn <- function (n)  sum (dbinom (m, n, emv.p[as.character(n)], log=TRUE))

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

## 3.  Calcula las estimaciones máximo verosímiles usando la verosimilitud
## y su logaritmo. Observa como los problemas de redondeo en los
## cálculos afectan a las estimaciones.

## sin log
imax <- which.max (sapply (x, Ln))
emv.n <- sop.n[imax]
emv.p[as.character(emv.n)]
## con log
imax <- which.max (sapply (x, logLn))
emv.n <- sop.n[imax]
emv.p[as.character(emv.n)]

## por método de los momentos
## mu = n.p; S = n.p.(1-p)  =>  p = (mu-S)/mu; n = mu^2/(mu-S)

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