title: Maximum likelihood estimation function for 3 parameters
body: <p>I'm attempting to make a maximum likelihood function for 3 parameters in a non-linear dataset. Two of my attempts are below, they are slightly different in approach, as well as the errors they produce. For context, the variable X contains a discrete set of time values and Y is the number of cases for a disease.</p>
  <pre><code>likelihood &lt;- function(params){ 
    params[1]
    params[2]
    params[3]
    top &lt;- (Y + params[1] - 1)
    mu &lt;- params[2] + params[3]*X
    
    result &lt;- sum(lchoose(top, Y)) + 
      sum(params[1]*log(params[1]/(params[1] + exp(mu)))) + 
      sum(Y*log(exp(mu)/(params[1] +exp(mu))))
    return(-result)
  }
  
  #initial param inputs
  param_initial &lt;- c(0.1,0.3,1)   
  
  # minimising the negative log-likelihood
  output &lt;- nlm(likelihood,
                p = beta_inital,
                hessian = T,
                gradtol = 1e-10,
                iterlim = 1000)
  # print model estimates for betas
  betas &lt;- output$estimate
  betas
  </code></pre>
  <p>Yielding: [1] 0.1 1.0</p>
  <p>(which are the same values as input)</p>
  <p>Second attempt:</p>
  <pre><code>log_like &lt;- function(params, Y, X) {
    X &lt;- as.matrix(X) ; Y &lt;- as.matrix(Y)
    params[1]
    params[2]
    params[3]
    
    result &lt;- sum(lchoose((Y + params[1] - 1), Y)) + 
      sum(params[1]*log(params[1]/(params[1] + exp(params[2] + params[3]*X)))) +
      sum(Y*log(exp(params[2] + params[3]*X) / (exp(params[2] + params[3]*X) + params[1])))
    return(-result)
  }
  
  MLE_estimates &lt;- nlm(log_like,
                       par = c(1,1,1),
                       hessian = TRUE, 
                       X,
                       Y)
  MLE_par &lt;- MLE_estimates$par
  MLESE &lt;- sqrt(diag(solve(MLE_estimates$hessian)))
  MLE &lt;- data.table(param = c(&quot;param1&quot;, &quot;param2&quot;, &quot;param3&quot;),
                    estimates = MLE_par,
                    sd = MLESE)
  </code></pre>
  <p>This is what the output under MLESE gives:</p>
  <blockquote>
  <p>Error in h(simpleError(msg, call)) :
  error in evaluating the argument 'x' in selecting a method for function 'diag': Lapack routine dgesv: system is exactly singular: U[1,1] = 0</p>
  </blockquote>
  <p>I also tried to plot this on a 3d graph using the following code:</p>
  <pre><code>log_like_graph &lt;- function(params,...){
    X &lt;- as.matrix(X) ; Y &lt;- as.matrix(Y)
    
    result &lt;- sum(lchoose((Y + params[1] - 1), Y)) + 
      sum(params[1]*log(params[1]/(params[1] + exp(params[2] + params[3]*X)))) +
      sum(Y*log(exp(params[2] + params[3]*X) / (exp(params[2] + params[3]*X) + params[1])))
    return(-result)
  }
  log_like_graph &lt;- Vectorize(log_like_graph)
  
  beta0val &lt;- seq(0,10,by=1)
  beta1val &lt;- seq(0,10,by=1)
  log_vals &lt;- outer(beta0val, beta1val, log_like_graph)
  
  persp(beta1val, beta0val, log_vals, theta = 7, phi = 8, r = 500) # doesnt work for 3 parameters
  </code></pre>
  <p>which yields</p>
  <blockquote>
  <p>Error in persp.default(beta1val, beta0val, log_vals, theta = 7, phi = 8,  :
  invalid 'z' limits
  In addition: Warning messages:
  1: In min(x) : no non-missing arguments to min; returning Inf
  2: In max(x) : no non-missing arguments to max; returning -Inf</p>
  </blockquote>
  <p>Safe to say I have no idea where I am messing this up, I have a feeling it may be in vector/matrix discrepancies.</p>