error message saying missing value where true/false is needed, but i don't see where the missing value comes from

Question

i keep getting the error message for my while statement at the end, but I don't see why there is a missing value in this statement. Should be a simple statement,but I am not sure why it is not working. Can someone please help me on that?

y = c(0.9683, 0.4515, 17.4488, 0.6287, 2.2330, 2.6467, 3.9589, 0.0782, 5.4717, 4.1161, 0.6715, 1.6350, 0.1640, 0.3331, 0.7501, 3.0846, 0.6889, 6.3826, 7.0869, 0.7967, 3.2684,0.1373, 2.8698, 1.5126, 0.9055)

x = c(0.1036, 2.1824, 0.1745, 2.0089, 1.2317, 0.6166, 0.4675, 3.2074, 0.0277, 1.2962, 0.6812, 0.1946, 1.3291, 0.4381, 0.2984, 0.3018, 0.7928, 0.2021, 1.0280, 0.0121, 1.2043, 2.9322, 1.4526, 0.6444, 0.1849)


score = function(b0,b1) {
i = 1:25
score1 = sum(-y[i]+(b0+b1*x[i])^-1)
score2 = sum(-y[i]*x[i]+x[i]*(b0+b1*x[i])^-1)

scorevector = c(score1,score2)
scorematrix = cbind(c(score1,score2))

scorematrix
}


jb = function(b0,b1) {
i=1:25

jb11 = sum((b0+b1*x[i])^{-2})
jb12 = sum(x[i]/(b0+b1*x[i])^{2})
jb21 = sum(x[i]/(b0+b1*x[i])^{2})
jb22 = sum(x[i]^{2}/(b0+b1*x[i])^{2})


jbmatrix = cbind(c(jb11,jb21),c(jb12,jb22))
jbmatrixinv = (1/(jb11*jb22-jb21*jb12))*cbind(c(jb22,-jb21),c(-jb12,jb11))

jbmatrixinv
}

bold = cbind(c(1,1)) ; bnew = cbind(c(2,2))
track = c(bnew[1,1],bnew[2,1])

while ((bnew[1,1]-bold[1,1])^2+(bnew[2,1]-bold[2,1])^2 > 10^{-6}) {
  
  bold = bnew
  bnew = bold + jb(bold[1,1],bold[2,1]) %*% score(bold[1,1],bold[2,1])
  track = rbind(track, c(bnew[1,1],bnew[2,1]))
}

Error in while ((bnew[1, 1] - bold[1, 1])^2 + (bnew[2, 1] - bold[2, 1])^2 >  : 
  missing value where TRUE/FALSE needed

Answer 1

Unlike what is said in comment, the left hand side is not constant, both bold and bnew will vary in the while loop. The LHS grows until it becomes NaN, and the condition NaN > 10^(-6) is now NA. The last value before the code breaks is 1.442032e+213.

I have rewritten the loop's first few instructions to have the LHS value be printed. As you can see it grows very quickly and breaks the code at the 9th iteration. And the growth comes from function jb. That's what you have to solve.

iter <- 0L
while (TRUE) {
  
  iter <- iter + 1L
  aux <- (bnew[1,1]-bold[1,1])^2 + (bnew[2,1]-bold[2,1])^2
  cat("iter:", iter, "\taux:", aux, "\n")
  
  if(!(aux > 10^{-6})) break
  
  bold = bnew
  bnew = bold + jb(bold[1,1],bold[2,1]) %*% score(bold[1,1],bold[2,1])
  track = rbind(track, c(bnew[1,1],bnew[2,1]))
}
#> iter: 1  aux: 2 
#> iter: 2  aux: 436.0846 
#> iter: 3  aux: 1754199 
#> iter: 4  aux: 5.234798e+13 
#> iter: 5  aux: 1.291931e+28 
#> iter: 6  aux: 7.904304e+56 
#> iter: 7  aux: 2.959099e+114 
#> iter: 8  aux: 1.442032e+213 
#> iter: 9  aux: NaN
#> Error in if (!(aux > 10^{: missing value where TRUE/FALSE needed

^{Created on 2022-01-31 by the reprex package (v2.0.1)}