title: L-BFGS-B needs finite values of 'fn' for a GTR model
body: <p>I am trying to estimate the parameters of the Generalized Time Reversible Model. My log likelihood function for the GTR model looks like this:</p>
  
  <pre><code>GTRLL &lt;- function(pars){
      Cij &lt;- countsMatrix()
      BL &lt;- pars[1]
      piA &lt;- pars[2]
      piC &lt;- pars[3]
      piG &lt;- pars[4]
      A &lt;- pars[5]
      B &lt;- pars[6]
      C &lt;- pars[7]
      D &lt;- pars[8]
      E &lt;- pars[9]
      F &lt;- 1
      piT &lt;- 1 - (piA+piC+piG)
      matSum &lt;- 0
  
      TT &lt;- (piC * A) + (piA*B) + (piG*C) * -1
      CC &lt;- (piT * A) + (piA*D) + (piG*E) * -1
      AA &lt;- (piT * B) + (piC*D) + (piG*F) * -1
      GG &lt;- (piT * C) + (piC*E) + (piA*F) * -1
  
      sigma &lt;- ((piA * AA) + (piT * TT) + (piC * CC) + (piG * GG)) * -1
  
      qMat = matrix(c((TT * sigma), (piT * A), (piT*B), (piT*C),(piC*A),(CC* sigma),(piC*D),(piC*E),(piA*B),(piA*D),(AA*sigma),(piA*F),(piG*C),(piG*E),(piG*F), (GG*sigma)),nrow = 4, ncol =4)
  
      Pij&lt;-expm(qMat*BL)
  
      for(j in 1:4){
          for(i in 1:4){
              matSum &lt;- matSum + (Cij[i,j] * log(Pij[i,j]))
          }
      }
      return(-matSum)
  }
  </code></pre>
  
  <p>I then attempt to run optim on this function using this:</p>
  
  <pre><code>optim(c(1,0.2,0.2,0.2,2,2,2,2,2), GTRLL, method = "L-BFGS-B", upper = c(Inf,0.25,0.25,0.25,Inf,Inf,Inf,Inf,Inf), lower = c(1,0.0001,0.0001,0.0001,1,1,1,1,1), control = list(trace = 5)
  </code></pre>
  
  <p>However, the return is:</p>
  
  <pre><code>Error in optim(c(1, 0.2, 0.2, 0.2, 2, 2, 2, 2, 2), GTRLL, method = "L-BFGS-B",  : 
    L-BFGS-B needs finite values of 'fn'
  </code></pre>
  
  <p>I have been trying to figure out a way to get this function to converge but it does not seem to be accomplishing that. I have attempted to play with the upper and lower bounds, but that does not seem to be working either.</p>
  
  <p><strong>Edit</strong>
  countsMatrix() returns a matrix of counts based on the sequence. In my example the returned matrix looks like this:</p>
  
  <pre><code>     [,1] [,2] [,3] [,4]
  [1,]    3    1    1    1
  [2,]    1    4    1    2
  [3,]    1    1    4    0
  [4,]    1    2    0    2
  </code></pre>
  
  <p>countsMatrix() looks like this:</p>
  
  <pre><code>countsMatrix &lt;- function(){
      S1 &lt;- "TTCAGTCACCTCTGAGTAA"
      S2 &lt;- "TACAGGCACGTCTCAGCAC"
      S1mat &lt;- strsplit(S1, "")[[1]]
      S2mat &lt;- strsplit(S2, "")[[1]]
      aa &lt;- 0
      at &lt;- 0
      ac &lt;- 0
      ag &lt;- 0
      ta &lt;- 0
      tt &lt;- 0
      tc &lt;- 0
      tg &lt;- 0
      ca &lt;- 0
      ct &lt;- 0
      cc &lt;- 0
      cg &lt;- 0
      ga &lt;- 0
      gt &lt;- 0
      GC &lt;- 0
      gg &lt;- 0
  
      for (i in 1:length(S1mat)){ 
          x &lt;- S1mat[i]
          y &lt;- S2mat[i]
          if (x == "A" &amp; y == "A"){
              aa &lt;- aa + 1
          }
          if (x == "A" &amp; y == "T"){
              at &lt;- at + 1
          }
          if (x == "A" &amp; y == "C"){
              ac &lt;- ac + 1
          }
          if (x == "A" &amp; y == "G"){
              ag &lt;- ag + 1
          }
          if (x == "T" &amp; y == "A"){
              ta &lt;- ta + 1
          }
          if (x == "T" &amp; y == "T"){
              tt &lt;- tt + 1
          }
          if (x == "T" &amp; y == "C"){
              tc &lt;- tc + 1
          }
          if (x == "T" &amp; y == "G"){
              tg &lt;- tg + 1
          }
          if (x == "C" &amp; y == "A"){
              ca &lt;- ca + 1
          }
          if (x == "C" &amp; y == "T"){
              ct &lt;- ct + 1
          }
          if (x == "C" &amp; y == "C"){
              cc &lt;- cc + 1
          }
          if (x == "C" &amp; y == "G"){
              cg &lt;- cg + 1
          }
          if (x == "G" &amp; y == "A"){
              ga &lt;- ga + 1
          }
          if (x == "G" &amp; y == "T"){
              gt &lt;- gt + 1
          }
          if (x == "G" &amp; y == "C"){
              GC &lt;- GC + 1
          }
          if (x == "G" &amp; y == "G"){
              gg &lt;- gg + 1
          }
      }
  
  
      Cij = matrix( c(tt, ct+ tc, at + ta, ct + tc , tc + ct, cc, ca + ac, cg + GC, ta + at, ac + ca, aa, ag+ga, tg+gt,GC+cg,ga+ag,gg), nrow = 4, ncol = 4)
  
      return(Cij)
  }
  </code></pre>