Dúvida em função de maximização (maxLik)

Prezados, estou tentando maximizar uma função de verossimilhança para estimar 5 parâmetros. Contudo, meus resultados estão dando NULL.

Sou iniciante na manipulação do R, por isso posso ter cometido erros primários. Favor alguém me ajude.

Meus dados são:

B <-c(1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,0,1,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,1,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0)

S <- c(0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,1,0,1,1,1,0,1,1,1,0,1,1,1,1,1,1,1,0,1,0,0,0,0,0,1,1,1,1,1,1,1,1,1)

A distribuição é de Poisson:

B <- rpois(50, 1)
S <- rpois(50, 1)

A equação que utilizo é:

loglik <- function(param) {
alpha <- param[1]
mu <- param[2]
Eb <- param[3]
Es <- param[4]
delta <- param[5]
ll <- (1-alpha)*exp(-Eb)*((Eb^B)/factorial(B))*exp(-Es)*((Es^S)/factorial(S)) + alpha*delta*exp(-Eb)*((Eb^B)/factorial(B))*exp(-(mu+Es))*((mu+Es)^S/factorial(S)) + alpha*(1-delta)*exp(-(mu+Eb))*((mu+Eb)^B/factorial(B))*exp(-Es)*((Es^S)/factorial(S))
}
res <- maxLik(loglik, grad=NULL, hess=NULL, start=c(0,1), constraints=NULL)
print(res)
coef(res)

Acho que estou errando nos argumentos dessa equação.

Alguém tem alguma sugestão?

Orleans Martins, seja bem vindo ao fórum!

Estive ausente nos últimos dias por isso não respondi antes, caso continue com a dúvida posso analisar os comandos acima para tentar ajudá-lo.

Prezado gustavoreis,

Nesses últimos dias tenho tentado entender essa função e ainda me restam algumas dúvidas. Sou leigo em R, estou utilizando porque preciso fazer esta maximização para o modelo que vou utilizar na minha tese e não encontrei outro software que faça esse cálculo.

A questão é a seguinte: preciso calcular uma probabilidade a partir de duas variáveis binárias (B e S). Cada uma contém 50 valores (simulados, apenas para calcular os parâmetros do modelo).

Eis a rotina que estou utilizando:

---INICIAR:

library(maxLik)


---DADOS PARA TESTE (50 observações)

B <- c(1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,0,1,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,1,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0)
S <- c(0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,1,0,1,1,1,0,1,1,1,0,1,1,1,1,1,1,1,0,1,0,0,0,0,0,1,1,1,1,1,1,1,1,1)

---EQUAÇÃO DIGITADA:

loglik <- function(param) {
alpha <- param[1]
mu <- param[2]
Eb <- param[3]
Es <- param[4]
delta <- param[5]
ll <- (1-alpha)*exp(-Eb)*((Eb^B)/factorial(B))*exp(-Es)*((Es^S)/factorial(S)) + alpha*delta*exp(-Eb)*((Eb^B)/factorial(B))*exp(-(mu+Es))*((mu+Es)^S/factorial(S)) + alpha*(1-delta)*exp(-(mu+Eb))*((mu+Eb)^B/factorial(B))*exp(-Es)*((Es^S)/factorial(S))
}
res <- maxLik(loglik, grad=NULL, hess=NULL, start=c(1,1,1,1,1), constraints=NULL)
print(res)
coef(res)


DÚVIDAS:

1) Observação: alpha, mu e delta são probabilidades (deveriam estar no intervalo entre 0 e 1?) Eb e Es são contínuos. Um colega que conhece um pouco mais de R disse que eles seriam "irrestritos" e as probabilidades "restritos". Não estou sabendo limitar isso.
2) start=c devem iniciar de 0 ou de 1? Quando coloco tudo 0 o software passa muito tempo calculando e não externa um resultado. Com 1 os parâmetros são calculados mas todos na casa de centena (por exemplo, alpha = 687.2159). Não estou sabendo analisar.
3) Sugeriram-me utilizar maxBFGS usando um espaço paramétrico. Não encontrei como.

Se puderem me ajudar agradecerei muito! Meu tempo está acabando e não estou encontrando solução.

Aceito qualquer negócio!