Logo Search packages:      
Sourcecode: robustbase version File versions  Download package

rf-common.f

c
c--   Routines common to
c--   fastLTS ( ./rfltsreg.f )  and
c--   fastMCD ( ./rffastmcd.f )
c
c
      subroutine rfrangen(n, nsel, index, seed)
c
c     Randomly draws nsel cases out of n cases.
c     Here, index is the index set.
c
      implicit none
      integer n, nsel, index(nsel), seed
c     real uniran
      double precision unifrnd
      integer i,j, num
c
      do 100 i=1,nsel
c     OLD 10        num=int(uniran(seed)*n)+1
 10    num=int(unifrnd()*n)+1
C      if(num .gt. n) then
C         call intpr('** rfrangen(): num > n; num=', -1, num, 1)
C         num=n
C      endif
       if(i.gt.1) then
          do 50 j=1,i-1
             if(index(j).eq.num) goto 10
 50       continue
       endif
       index(i)=num
 100  continue
      return
      end
ccccc
ccccc
cOLD  function uniran(seed)
cOLD cc
cOLD cc  Draws a random number from the uniform distribution on [0,1].
cOLD cc
cOLD  real uniran
cOLD  integer seed
cOLD  integer quot
cOLD cc
cOLD  seed=seed*5761+999
cOLD  quot=seed/65536
cOLD  seed=seed-quot*65536
cOLD  uniran=float(seed)/65536.D0
cOLD  return
cOLD  end
ccccc
ccccc
      subroutine rfgenpn(n,nsel,index)
cc
cc    Constructs all subsets of nsel cases out of n cases.
cc
        implicit none
      integer n,nsel,index(nsel)
cc
        integer k,i

      k=nsel
      index(k)=index(k)+1
c while
 10   if(k.eq.1 .or. index(k).le.(n-(nsel-k))) goto 100
      k=k-1
      index(k)=index(k)+1
      do 50 i=k+1,nsel
        index(i)=index(i-1)+1
 50   continue
      goto 10
c end{while}
 100  return
      end
ccccc
ccccc
      subroutine rfshsort(a,n)
cc
cc  Sorts the array a of length n.
cc
        implicit none
        integer n
      double precision a(n)
c
      double precision t
      integer gap, i,j, nextj

      gap=n
 100  gap=gap/2
      if(gap.eq.0) goto 200
      do 180 i=1,n-gap
        j=i
 120    if(j.lt.1) goto 180
        nextj=j+gap
        if(a(j).gt.a(nextj)) then
          t=a(j)
          a(j)=a(nextj)
          a(nextj)=t
        else
          j=0
        endif
        j=j-gap
        goto 120
 180  continue
      goto 100
 200  return
      end
ccccc
ccccc
      subroutine rfishsort(a,kk)
cc
cc  Sorts the integer array a of length kk.
cc
        implicit none
      integer kk, a(kk)
c
      integer t, gap, i,j, nextj

      gap=kk
 100  gap=gap/2
      if(gap.eq.0) goto 200
      do 180 i=1,kk-gap
        j=i
 120    if(j.lt.1) goto 180
        nextj=j+gap
        if(a(j).gt.a(nextj)) then
          t=a(j)
          a(j)=a(nextj)
          a(nextj)=t
        else
          j=0
        endif
        j=j-gap
        goto 120
 180  continue
      goto 100
 200  return
      end
ccccc
ccccc
      function replow(k)
cc
cc    Find out which combinations of n and p are
cc    small enough in order to perform exaustive search
cc    Returns the maximal n for a given p, for which
cc    exhaustive search is to be done
cc
cc    k is the number of variables (p)
cc
      implicit none
      integer replow, k
c
      integer irep(6)
      data irep/500,50,22,17,15,14/
c
      if(k .le. 6) then
       replow = irep(k)
      else
       replow = 0
      endif
      return
      end
ccccc
ccccc
      function rfncomb(k,n)
cc
cc  Computes the number of combinations of k out of n.
cc  (To avoid integer overflow during the computation,
cc  ratios of reals are multiplied sequentially.)
cc  For comb > 1E+009 the resulting 'comb' may be too large
cc  to be put in the integer 'rfncomb', but the main program
cc  only calls this function for small enough n and k.
cc
      implicit none
      integer rfncomb,k,n
c
      double precision comb,fact
      integer j
c
      comb=dble(1.0)
      do 10 j=1,k
         fact=(dble(n-j+1.0))/(dble(k-j+1.0))
         comb=comb*fact
 10   continue
      rfncomb=int(comb+0.5D0)
      return
      end
ccccc
ccccc
      subroutine rfcovcopy(a,b,n1,n2)
cc
cc  Copies matrix a to matrix b.
cc
      double precision a(n1,n2)
      double precision b(n1,n2)
c
      do 100 i=1,n1
        do 90 j=1,n2
          b(i,j)=a(i,j)
 90     continue
 100  continue
      return
      end
ccccc
ccccc
      function rffindq(aw,ncas,k,index)
cc
cc  Finds the k-th order statistic of the array aw of length ncas.
cc
cc MM{FIXME}: "rather" use R
































































's C API   rPsort (double* X, int N, int K)        implicit none        integer ncas,k,index(ncas)    double precision rffindq, aw(ncas)c double precision ax,wa        integer i,j,l,lr,jnccc  do 10 j=1,ncas      index(j)=j 10   continue    l=1   lr=ncasc    while(l < lr) 20  if(l.ge.lr) goto 90     ax=aw(k)    jnc=l j=lr 30     if(jnc.gt.j) goto 80 40 if(aw(jnc).ge.ax) goto 50     jnc=jnc+1   goto 40 50  if(aw(j).le.ax) goto 60 j=j-1 goto 50 60  if(jnc .le. j) then           i=index(jnc)           index(jnc)=index(j)           index(j)=i           wa=aw(jnc)           aw(jnc)=aw(j)           aw(j)=wa           jnc=jnc+1           j=j-1        endif    goto 30 80  if(j.lt.k) l=jnc  if(k.lt.jnc) lr=j goto 20 90  rffindq=aw(k)     return      endcccccccccc     subroutine rfrdraw(a,n,seed,ntot,mini,ngroup,kmini)cccc  Draws ngroup nonoverlapping subdatasets out of a dataset of size n,cc  such that the selected case numbers are uniformly distributed from 1 to n.cc        implicit none   integer ntot, kmini, a(2,ntot), n, seed, mini(kmini), ngroupc        double precision unifrndc        integer jndex, nrand, k,m,i,jcc     jndex=0     do 10 k=1,ngroup    do 20 m=1,mini(k)cOLD           nrand=int(uniran(seed)*(n-jndex))+1       nrand=int(unifrnd()*(n-jndex))+1C         if(nrand .gt. n-jndex) thenC           call intpr(C      1         '** rfrdraw(): need to correct nrand > n-jndex; nrand=





























































',C      2                    -1, nrand, 1)C           nrand=n-jndexC       endif       jndex=jndex+1     if(jndex.eq.1) then       a(1,jndex)=nrand        a(2,jndex)=k          else          a(1,jndex)=nrand+jndex-1            a(2,jndex)=k            do 5,i=1,jndex-1          if(a(1,i).gt.nrand+i-1) then                do 6, j=jndex,i+1,-1                  a(1,j)=a(1,j-1)               a(2,j)=a(2,j-1) 6           continue                a(1,i)=nrand+i-1              a(2,i)=k                goto 20c                   ------- break            endif 5         continue        endif 20        continue 10     continue    return      endcccccccccc     function rfodd(n)cc     logical rfoddcc   rfodd=.true.      if(2*(n/2).eq.n) rfodd=.false.      return      endcccccc unused  function rfnbreak(nhalf,n,nvar)c unused ccc unused cc  Computes the breakdown value - in percent! - of the MCD estimatorc unused ccc unused         implicit nonec unused     integer rfnbreak, nhalf, n, nvarc unusedc unused      if (nhalf.le.(n+nvar+1)/2) thenc unused     rfnbreak=(nhalf-nvar)*100/nc unused     elsec unused        rfnbreak=(n-nhalf+1)*100/nc unused      endifc unused     returnc unused    endccccc    subroutine rfmcduni(w,ncas,jqu,slutn,bstd,aw,aw2,factor,len)cccc  rfmcduni : calculates the MCD in the univariate case.cc          w contains the ordered observationsccc This version returns the index (jint) in 'len















































Generated by  Doxygen 1.6.0   Back to index