Numerically solve estimating equations for shift estimator of the Van der Waerden test (normal scores)

#Data from HMC p557.

s1<-c(51.9,56.9,45.2,52.3,59.5,41.4,46.4,45.1,53.9,42.9,41.5,55.2,32.9,54.0,45.0)
s2<-c(59.2,49.1,54.4,47.0,55.9,34.9,62.2,41.6,59.3,32.7,72.1,43.8,56.8,76.7,60.3)

L<-vector("list",21)
L2<-vector("list",21)
L_rank<-vector("list",21) ##21 is from the steps 4 to 6 by 0.1

for(d in seq(from=4, to=6, by=0.1)){
L[[round((d-3.9)/0.1)]]<-(s2-d)
L2[[round((d-3.9)/0.1)]]<-c(s1,L[[round((d-3.9)/0.1)]])
L_rank[[round((d-3.9)/0.1)]]<-rank(c(s1,L[[round((d-3.9)/0.1)]]))
}
L_w<-vector("list",21)
L_ns<-vector("list",21)
L_sum_w<-vector("list",21)

for (i in 1:21){
 for (j in 16:30){

L_w[[i]][j-15]<-L_rank[[i]][j]

L_ns[[i]][j-15]<-qnorm(L_w[[i]][j-15]/31,0,1)

L_sum_w[[i]]<-sum(L_ns[[i]])
 
  }
}


####
[[1]]
[1] 0.6743001

[[2]]
[1] 0.3482201

[[3]]
[1] 0.2545677

[[4]]
[1] 0.2545677

[[5]]
[1] 0.2545677

[[6]]
[1] 0.21268

[[7]]
[1] 0.171218

[[8]]
[1] 0.171218

[[9]]
[1] 0.171218

[[10]]
[1] 0.1301031

[[11]]
[1] 0.08926116

[[12]]
[1] -0.0229044

[[13]]
[1] -0.06941334

[[14]]
[1] -0.1639976

[[15]]
[1] -0.4882569

[[16]]
[1] -0.5844623

[[17]]
[1] -0.6919131

[[18]]
[1] -0.7492664

[[19]]
[1] -0.7492664

[[20]]
[1] -0.7492664

[[21]]
[1] -0.7492664


We can see the 12th  value =-0.0229044 is the closest value to zero, the corresponded delta is 5.1. (form the 4 to 6 by 0.1 step)

Therefore, 5.1 is the solution of the Estimating  Equations.