0

I am trying to replicate the below R code to estimate parameters using Maximum Likelihood method in Python.

I want to obtain the same results using both the codes, but my estimated values are different, I am not sure if both the codes are optimising the same parameters.

R-Code:

ll <- function (prop, numerator, denominator) {
  return(
    lgamma(denominator + 1) - 
      lgamma(  numerator + 1) - 
      lgamma(denominator - numerator + 1) +
      numerator * log(prop) + (denominator - numerator) * log(1 - prop)
  )
}

compLogLike <- function(pvec){
  return(sum(ll(pvec, dat$C, dat$N)))
}


fct_p_ll <- function(a, c0, c1){
  xa_ <- exp(c0 + c1*c(20, a))
  return(1 - exp((xa_[1]-xa_[2])/c1))
}


fct_ll <- function(x){
  pv <- sapply(22.5+5*(0:8), FUN = fct_p_ll, c0 = x[1], c1 = x[2])
  return(compLogLike(pv))
}

opt.res <- optim(par = c(-9.2, 0.07), fn = fct_ll, control = list(fnscale = -1.0), hessian = TRUE)


fisherInfo <-  solve(-opt.res$hessian)
propSigma  <-  sqrt(diag(fisherInfo))
upper      <-  opt.res$par+1.96*propSigma
lower      <-  opt.res$par-1.96*propSigma
interval   <-  data.frame(val = opt.res$par, ci.low=lower, ci.up = upper)

Python Code:

def ll(prop, numerator, denominator):
    print(prop, numerator, denominator)
    if prop > 0:
        value = (math.lgamma(denominator + 1) - 
          math.lgamma(  numerator + 1) - 
          math.lgamma(denominator - numerator + 1) +
          numerator * math.log(prop) + (denominator - numerator) * math.log(1 - prop))
        return  value
    return 0

def compLogLike(pvec):
    p = list(pvec)
    c = list(df["C"])
    n = list(df["N"])
    compLog = 0
    for idx, val in enumerate(p):
        compLog += ll(p[idx],c[idx],n[idx])
    print(compLog)
    return compLog

def fct_p_ll(a,c0,c1):
    val_list = [c1 * val for val in [20, a]]
    xa_ = np.exp([c0 + val for val in val_list])
    return 1 - np.exp((xa_[0] -xa_[1])/c1)

def fct_ll(x):
    ages_1 = np.arange(22.5, 67.5 , 5)
    pv = [fct_p_ll(a=val,c0=x[0],c1=x[1]) for val in ages_1]
    return compLogLike(pv)

opt = minimize(fct_ll, [-9.2, 0.07],  method='Nelder-Mead', hess=Hessian(lambda x: fct_ll(x,a)))

Any inputs would be really helpful.

Akshaya
  • 11

0 Answers0