0
$\begingroup$

I am running a simple Newton-Raphson problem in R as given below. However I am getting very weird results. Can anyone tell me where I am going wrong?

    N=100
x=rnorm(N)
y=rnorm(N)
mu1=1
sigma1=1

maxiter=100
   out2 <- matrix(NA, nrow=maxiter+1,ncol=2)
   f_mu1=sum(((x-mu1)/sigma1^2))
   f_mu1[!is.finite(f_mu1)]=0

   out2[1,] <- c(1,f_mu1)
   i <- 1
  continue <- T

while(continue) {

i=i+1
mu1.old=mu1

df_mu1=sum((-1/sigma1^2))

mu1=mu1-(f_mu1)*(df_mu1)^(-1)
f_mu1=sum(((x-mu1)/sigma1^2))
f_mu1[!is.finite(f_mu1)]=0
out2[i,] <- c(mu1,f_mu1)

continue <- (abs(mu1-mu1.old) > 0.00001) &&
                   (i <= maxiter)
#print(mu1)
   }

One set of results that I got is as under, just to give an idea how weird it is

               [,1]           [,2]
  [1,]   1.000000e+00  -8.645455e+01
  [2,]  -8.545455e+01   8.559001e+03
  [3,]   8.473546e+03  -8.473411e+05
  [4,]  -8.388675e+05   8.388677e+07
  [5,]   8.304790e+07  -8.304790e+09
  [6,]  -8.221742e+09   8.221742e+11
  [7,]   8.139525e+11  -8.139525e+13
  [8,]  -8.058129e+13   8.058129e+15
  [9,]   7.977548e+15  -7.977548e+17
 [10,]  -7.897773e+17   7.897773e+19
 [11,]   7.818795e+19  -7.818795e+21
 [12,]  -7.740607e+21   7.740607e+23
 [13,]   7.663201e+23  -7.663201e+25
 [14,]  -7.586569e+25   7.586569e+27
 [15,]   7.510703e+27  -7.510703e+29
 [16,]  -7.435596e+29   7.435596e+31
 [17,]   7.361240e+31  -7.361240e+33
 [18,]  -7.287628e+33   7.287628e+35
 [19,]   7.214752e+35  -7.214752e+37
 [20,]  -7.142604e+37   7.142604e+39
 [21,]   7.071178e+39  -7.071178e+41
 [22,]  -7.000466e+41   7.000466e+43
 [23,]   6.930462e+43  -6.930462e+45
 [24,]  -6.861157e+45   6.861157e+47
 [25,]   6.792545e+47  -6.792545e+49
 [26,]  -6.724620e+49   6.724620e+51
 [27,]   6.657374e+51  -6.657374e+53
 [28,]  -6.590800e+53   6.590800e+55
 [29,]   6.524892e+55  -6.524892e+57
 [30,]  -6.459643e+57   6.459643e+59
 [31,]   6.395047e+59  -6.395047e+61
 [32,]  -6.331096e+61   6.331096e+63
 [33,]   6.267785e+63  -6.267785e+65
 [34,]  -6.205107e+65   6.205107e+67
 [35,]   6.143056e+67  -6.143056e+69
 [36,]  -6.081626e+69   6.081626e+71
 [37,]   6.020809e+71  -6.020809e+73
 [38,]  -5.960601e+73   5.960601e+75
 [39,]   5.900995e+75  -5.900995e+77
 [40,]  -5.841985e+77   5.841985e+79
 [41,]   5.783566e+79  -5.783566e+81
 [42,]  -5.725730e+81   5.725730e+83
 [43,]   5.668473e+83  -5.668473e+85
 [44,]  -5.611788e+85   5.611788e+87
 [45,]   5.555670e+87  -5.555670e+89
 [46,]  -5.500113e+89   5.500113e+91
 [47,]   5.445112e+91  -5.445112e+93
 [48,]  -5.390661e+93   5.390661e+95
 [49,]   5.336754e+95  -5.336754e+97
 [50,]  -5.283387e+97   5.283387e+99
 [51,]   5.230553e+99 -5.230553e+101
 [52,] -5.178247e+101  5.178247e+103
 [53,]  5.126465e+103 -5.126465e+105
 [54,] -5.075200e+105  5.075200e+107
 [55,]  5.024448e+107 -5.024448e+109
 [56,] -4.974204e+109  4.974204e+111
 [57,]  4.924462e+111 -4.924462e+113
 [58,] -4.875217e+113  4.875217e+115
 [59,]  4.826465e+115 -4.826465e+117
 [60,] -4.778200e+117  4.778200e+119
 [61,]  4.730418e+119 -4.730418e+121
 [62,] -4.683114e+121  4.683114e+123
 [63,]  4.636283e+123 -4.636283e+125
 [64,] -4.589920e+125  4.589920e+127
 [65,]  4.544021e+127 -4.544021e+129
 [66,] -4.498581e+129  4.498581e+131
 [67,]  4.453595e+131 -4.453595e+133
 [68,] -4.409059e+133  4.409059e+135
 [69,]  4.364968e+135 -4.364968e+137
 [70,] -4.321319e+137  4.321319e+139
 [71,]  4.278106e+139 -4.278106e+141
 [72,] -4.235325e+141  4.235325e+143
 [73,]  4.192971e+143 -4.192971e+145
 [74,] -4.151042e+145  4.151042e+147
 [75,]  4.109531e+147 -4.109531e+149
 [76,] -4.068436e+149  4.068436e+151
 [77,]  4.027751e+151 -4.027751e+153
 [78,] -3.987474e+153  3.987474e+155
 [79,]  3.947599e+155 -3.947599e+157
 [80,] -3.908123e+157  3.908123e+159
 [81,]  3.869042e+159 -3.869042e+161
 [82,] -3.830352e+161  3.830352e+163
 [83,]  3.792048e+163 -3.792048e+165
 [84,] -3.754128e+165  3.754128e+167
 [85,]  3.716586e+167 -3.716586e+169
 [86,] -3.679420e+169  3.679420e+171
 [87,]  3.642626e+171 -3.642626e+173
 [88,] -3.606200e+173  3.606200e+175
 [89,]  3.570138e+175 -3.570138e+177
 [90,] -3.534437e+177  3.534437e+179
 [91,]  3.499092e+179 -3.499092e+181
 [92,] -3.464101e+181  3.464101e+183
 [93,]  3.429460e+183 -3.429460e+185
 [94,] -3.395166e+185  3.395166e+187
 [95,]  3.361214e+187 -3.361214e+189
 [96,] -3.327602e+189  3.327602e+191
 [97,]  3.294326e+191 -3.294326e+193
 [98,] -3.261383e+193  3.261383e+195
 [99,]  3.228769e+195 -3.228769e+197
[100,] -3.196481e+197  3.196481e+199
[101,]  3.164516e+199 -3.164516e+201
> 
$\endgroup$
  • 3
    $\begingroup$ As posted this is really hard to follow, but i'd love to help. Before I do could you: 1) State exactly what function you are trying to minimize, maximize or determine the roots of with this newton iteration, this will help set context for my reading. 2) Properly indent your code. 3) Clean up your use of = vs <-? $\endgroup$ – Matthew Drury May 26 '15 at 22:55
  • 2
    $\begingroup$ Newton-Raphson, aka Newtons method, has something called a "radius of convergence". If your initial conditions start within that radius then it is guaranteed to converge, and do so quadratically. (yay). If your initial condition is outside then one of three things occur: 1) it rapidly converges toward an infinite value (you are getting this), 2) it wanders in finite-value land forever without converging, or 3) it accidentally jumps into the radius of convergence and converges. $\endgroup$ – EngrStudent - Reinstate Monica May 26 '15 at 23:13
  • 2
    $\begingroup$ This might possibly be on topic, but either way please explain the problem you're implementing, not in R-code but in words with algebra and (briefly) what the algorithmic steps are you're implementing in simple terms understandable to a user of any software/programming language. $\endgroup$ – Glen_b -Reinstate Monica May 27 '15 at 0:17
3
$\begingroup$

If I understand to code correctly, this looks like only a bug. It is trying to find the roots for

$$ f(\mu_1) = \sum_{i=1}^{N}\frac{x-\mu_1}{\sigma_1^2} $$

with $x\sim N(0,1)$, $N=100$, and $\sigma_1=1$. The derivative is

$$ \frac{\partial f}{\partial \mu_1} = \sum_{i=1}^{N}\frac{-1}{\sigma_1^2} $$

The R code is df_mu1=sum((-1/sigma1^2)) which looks correct. Unfortunately, sigma1 is just a single number and not a vector so the sum only returns -1/sigma1^2. As the comments point out, the correct fix depends on your original intentions. One way to fix this is to set df_mu1=N * (-1/sigma1^2) instead. When I do this, it converges in 2 iterations and mu1 converges to mean(x)

$\endgroup$
  • $\begingroup$ thats exactly the problem. Why does sum (1/sigma^2) does not work out to N*(1/sigma^2) in R. $\endgroup$ – indu mann May 27 '15 at 2:55
  • $\begingroup$ @Glen I am doing this newton raphson step as a part of a bigger EM algorithm where the M step is Newton Raphson. so my log likelihood is Q=sum(sum(log(f)))+sum(log(lambdacopula density))-Nsum(lambda^2)+langrange multiplier*(sum(lambda's)-1) where f is marginal density. The E step is estimation of lambda's and M step is getting the marginal parameters mu1, sigma1 which then is again used in the estimation of lambda inside the loop till the value of loglik converges. $\endgroup$ – indu mann May 27 '15 at 2:56
  • $\begingroup$ @indumann Why would it? In sigma1=1 you explicitly told R that sigma is a vector of length 1. So 1/sigma^2 is also a vector of length 1, and sum(1/sigma^2) is a sum of a single number, so it is just the number. For example sum(2) = 2. $\endgroup$ – Matthew Drury May 27 '15 at 3:10
  • $\begingroup$ As a side comment, and sorry to be slightly rude, but if this is a small part of a bigger EM algorithm you are going to have some serious issues with maintainability and understandability down the line with the state of your code. Please clean it up for the sake of your future self. $\endgroup$ – Matthew Drury May 27 '15 at 3:11
  • $\begingroup$ @MatthewDrury Thanks for the guidance. As I am in the initial stages of learning a programming language, my code appears to be messy. It would be very nice of you if you could tell me what exactly is meant by 'cleaning up'. what is so wrong in my code that may cost me down the line. $\endgroup$ – indu mann May 27 '15 at 4:33
0
$\begingroup$

This generally happens when the method hits a near-stationary point. Notice that if the derivative is close to 0, the inverse of the derivative will blow up and that is part of your update term for mu_1

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.