R/nelder.R

In musto101/wilcox_R: What the package does (short line)

nelder <-
function(x,N,FN,START=c(rep(1,N)),STEP=c(rep(1,N)),
XMIN=c(rep(0,N)),XSEC=c(rep(0,N))){
#     NELDER-MEAD method for minimzing a function
#
#     TAKEN FROM OLSSON, J QUALITY TECHNOLOGY, 1974, 6, 56.
#
#     x= n by p matrix containing data; it is used by
#        function to be minimized.
#     N= number of parameters
#
#     FN=the function to be minimized
#     FORM: FN(x,theta), theta is vector containing
#     values for N parameters.
#
#     START = starting values.
#     STEP=initial step.
#     This function returns the N values for theta that minimize FN
#
      ICOUNT<-500
      REQMIN<-.0000001
      NN<-N+1
      P<-matrix(NA,nrow=N,ncol=NN)
      P[,NN]<-START
      PBAR<-NA
      RCOEFF<-1
      ECOEFF<-2
      CCOEFF<-.5
      KCOUNT<-ICOUNT
      ICOUNT<-0
      DABIT<-2.04067e-35
      BIGNUM<-1.e38
      KONVGE<-5
      XN<-N
      DN<-N
      Y<-rep(0,NN)
      Y[NN]<-FN(x,START)
      ICOUNT<-ICOUNT+1
      for(J in 1:N){
      DCHK<-START[J]
      START[J]<-DCHK+STEP[J]
      for(I in 1:N){
      P[I,J]<-START[I]
}
      Y[J]<-FN(x,START)
      ICOUNT<-ICOUNT+1
      START[J]<-DCHK
}
      I1000<-T
       while(I1000){
      YLO<-Y[1]
      YNEWLO<-YLO
      ILO<-1
      IHI<-1
      for(I in 2:NN){
      if(Y[I] <  YLO){
      YLO<-Y[I]
      ILO<-I}
      if(Y[I] > YNEWLO){
      YNEWLO<-Y[I]
      IHI<-I}
}
      DCHK<-(YNEWLO+DABIT)/(YLO+DABIT)-1
      if(abs(DCHK) < REQMIN){
      I1000<-F
      next
}
      KONVGE<-KONVGE-1
      if(KONVGE == 0){
      KONVGE<-5
      for(I in 1:N){
      COORD1<-P[I,1]
      COORD2<-COORD1
      for(J in 2:NN){
      if(P[I,J] < COORD1)COORD1<-P[I,J]
      if(P[I,J] > COORD2)COORD2<-P[I,J]
}     # 2010 CONTINUE
      DCHK<-(COORD2+DABIT)/(COORD1+DABIT)-1
      if(abs(DCHK) > REQMIN)break
}
}
     if(ICOUNT >= KCOUNT){
      I1000<-F
      next
}
      for(I in 1:N){
      Z<-0.0
      Z<-sum(P[I,1:NN]) # 6
      Z<-Z-P[I,IHI]
  PBAR[I]<-Z/DN
}
    PSTAR<-(1.+RCOEFF)*PBAR-RCOEFF*P[,IHI]
      YSTAR<-FN(x,PSTAR)
      ICOUNT<-ICOUNT+1
      if(YSTAR < YLO && ICOUNT >= KCOUNT){
       P[,IHI]<-PSTAR
       Y[IHI]<-YSTAR
       next
}
  IFLAG<-T
      if(YSTAR < YLO){
    P2STAR<-ECOEFF*PSTAR+(1-ECOEFF)*PBAR
      Y2STAR<-FN(x,P2STAR)
      ICOUNT<-ICOUNT+1
      if(Y2STAR >= YSTAR){
       P[,IHI]<-PSTAR
       Y[IHI]<-YSTAR
       next #In essence, go to 19 which goes to 1000
}
      IFLAG<-T
      while(YSTAR < Y[IHI]){
      P[,IHI]<-P2STAR
      Y[IHI]<-Y2STAR
      IFLAG<-F
     break
     L<-sum(Y[1:NN] > YSTAR)
      if(L > 1){
       P[,IHI]<-PSTAR
       Y[IHI]<-YSTAR
       IFLAG<-T
       break
}
       if(L > 1)break # go to 19
      if(L != 0){
      P[1:N,IHI]<-PSTAR[1:N]
      Y[IHI]<-YSTAR
}
I1000<-F
break
  if(ICOUNT >= KCOUNT){
      I1000<-F
      next
}
   P2STAR[1:N]<-CCOEFF*P[1:N,IHI]+(1-CCOEFF)*PBAR[1:N]
      Y2STAR<-FN(x,P2STAR)
      ICOUNT<-ICOUNT+1
}   # END WHILE
}
if(IFLAG){
for(J in 1:NN){
P[,J]<-(P[,J]+P[,ILO])*.5
   XMIN<-P[,J]
      Y[J]<-FN(x,XMIN)
}
      ICOUNT<-ICOUNT+NN
      if(ICOUNT < KCOUNT)next
      I1000<-F
next
}
      P[1:N,IHI]<-PSTAR[1:N]
      Y[IHI]<-YSTAR
}
    for(J in 1:NN){
      XMIN[1:N]<-P[1:N,J]
}
      Y[J]<-FN(x,XMIN)
      YNEWLO<-BIGNUM
  for(J in 1:NN){
      if (Y[J] < YNEWLO){
      YNEWLO<-Y[J]
      IBEST<-J
}}
      Y[IBEST]<-BIGNUM
      YSEC<-BIGNUM
for(J in 1:NN){
if(Y[J] < YSEC){
      YSEC<-Y[J]
      ISEC<-J
}}
      XMIN[1:N]<-P[1:N,IBEST]
      XSEC[1:N]<-P[1:N,ISEC]
XMIN
}