title: Why is my code implementing the Fisher scoring algorithm failing to converge?
body: <pre><code># Set data sets
  data &lt;- MASS::birthwt
  X &lt;- as.matrix(cbind(1, data[, -1]))
  Y &lt;- data[, 1]
  n &lt;- dim(X)[1]
  p &lt;- dim(X)[2]    
  beta &lt;- rep(0, p)  # Initialize beta
  beta1 &lt;- rep(1, p)  # Initialize beta1
  
  
  # fit logit regression using Fisher scoring
  while (sum(abs(beta - beta1))  &gt; 0.01) {
    beta1 &lt;-beta
  
    # Calculate mu and eta
    eta &lt;- X %*% beta1
    mu &lt;- exp(eta) / (1 + exp(eta))
  
    # Calculate the derivatives
    dmu_deta &lt;- exp(eta) / (( 1 + exp(eta)) ^ 2)
    deta_dmu &lt;- 1 / mu + 1 / (1 - mu)
  
    V &lt;- mu * (1 - mu)  # Calculate the variance function
    W &lt;- diag((dmu_deta ^ 2 / V)[,])  # Calculate the weight matrix
    Z &lt;- eta + (Y - mu) * deta_dmu  # Calculate Z
  
    beta &lt;- solve(t(X) %*% W %*% X) %*% (t(X) %*% W %*% Z)
  }
  </code></pre>
  
  <p>I am using Iteratively Reweighted Least Square method. The <code>X</code> and <code>Y</code> come from the built-in dataset <code>birthwt</code>. I do not understand why this method does not converge. It always returns a <code>NaN</code>. But when I remove the intercept, it converges. I know that I can simply use <code>glm</code>, but I would like to understand the implementation.</p>