Title                  :SELECT
Keywords               :STATISTICS,MEAN VALUE,SELECTION TEST
Computer               :DEC LSI-11/73
Operating System       :RSX-11M-PLUS
Programming Language   :DEC-FORTRAN IV
Hardware Requirements  :None
Author(s)              :Dr. Pierre TREBBIA
Correspondence Address :Laboratoire de Physique des Solides
                       :Bat. 510, F91405 ORSAY CEDEX FRANCE
Abstract: Let us suppose that we obtained N estimations Vi (i=1,N) of a same 
unknown X. If we know the variance and covariance matrix of the Vi, one can 
compute the mean value (estimation of X) and the confidence limits of this 
mean. One can then check whether the Vi are valid measurements or whether one 
of them (at least) must be removed from the measurement set.
      This program SELECT is of general purpose and can be used as a valuable 
substitute to the general method of the "bad-looking" rejection criterion.
      A complete description of the algorithm is given in "Unbiased 
Method for Signal Estimation in Electron Energy Loss Spectroscopy...", 
P. TREBBIA, Ultramicroscopy (1988).
Title                  :SELECT
Keywords               :STATISTICS,MEAN VALUE,SELECTION TEST
Computer               :DEC LSI-11/73
Operating System       :RSX-11M-PLUS
Programming Language   :DEC-FORTRAN IV
Hardware Requirements  :None
Author(s)              :Dr. Pierre TREBBIA
Correspondence Address :Laboratoire de Physique des Solides
                       :Bat. 510, F91405 ORSAY CEDEX FRANCE
Description : see comments line in the source code text
References : see comments line in the source code text
Compilation procedure : compile the main program SELECT and the 2 subroutines 
                        PERMAT and DPMEAN.
Linking procedure : Link the 3 resulting object files.
Test data : SELECT.DAT is the output that SELECT should give 
            if the appropriate parameters are entered.
General comments : If the compiler is not a DEC-FORTRAN IV compiler, do check 
                   the source code : for example, the ! sign in a line means 
                   that the remaining text is comments, or the LOG function is 
                   implicitely changed to DLOG if the variable is double 
                   precision.
Title                  :SELECT
Keywords               :STATISTICS,MEAN VALUE,SELECTION TEST
Computer               :DEC LSI-11/73
Operating System       :RSX-11M-PLUS
Programming Language   :DEC-FORTRAN IV
Hardware Requirements  :None
Author(s)              :Dr. Pierre TREBBIA
Correspondence Address :Laboratoire de Physique des Solides
                       :Bat. 510, F91405 ORSAY CEDEX FRANCE

       PROGRAM SELECT
c
c       LAST VERSION : 25-JUL-1987 08:43:50
c
c       purpose
c       -------
c       This program estimates the weighted mean & standard deviation from a set
c       of dependent OR independent estimations of the SAME unknown.
c       A selection, by chi-square tests at levels 0.1, 0.05 and 0.01, of the
c       input data can be made upon request.
c
c       subroutines and function subprograms required
c       ---------------------------------------------
c       dpmean,permat
c
c       comments
c       --------
c       1) This program cannot handle more than 22 data pairs (estimation,
c       standard deviation).
c       2) This program estimates the weighted mean at 4 different levels :
c               level=0 => confidence 100%      => risk 0%
c               level=1 => confidence 90%       => risk 10%
c               level=2 => confidence 95%       => risk 5%
c               level=3 => confidence 99%       => risk 1%
c       3) Modify PARAMETER statements in BOTH SELECT and DPMEAN programs to
c          increase these limits, AND enter appropriate reduced chi-square
c          & risk values in arrays CHI and RISK with the DATA instructions.
c       4) The program checks data in increasing level order and stops if the
c          whole set of data is valid at a given level greater or equal to 1.
c       5) After selection at level 1, the whole data set (values, standard
c          deviations and covariances) are reorganized and output in decreasing
c          order of selection priority.
c       6) Enter CTRL Z for exit.
c       7) For building an executive task, one must :
c               a) compile the following FORTRAN programs :
c                       SELECT, DPMEAN, PERMAT
c               b) link all the resulting object files.
c
c       author :
c       --------
c       Pierre TREBBIA
c       US 41 : "Microscopie Electronique Analytique Quantitative"
c       Laboratoire de Physique des Solides, Bat. 510
c       Universite Paris-Sud, F91405 ORSAY Cedex
c       Phone : (33-1) 69 41 53 68
c
        PARAMETER (NMAX=22)
        PARAMETER (LMAX=3)
        double precision sumwgt,wmean,error,vmean,chi2,diff
        dimension data(nmax,2),error(nmax,nmax),risk(0:lmax)
        dimension chi(nmax-2,lmax),wmean(0:lmax)
        dimension sdmean(0:lmax),nvalid(0:lmax),nfree(0:lmax)
        dimension diff(nmax),nord(nmax),vmean(0:lmax),chi2(0:lmax)
        logical pr,ind,sel,check,next
        DATA risk/0.,0.1,0.05,0.01/
c
c       Put in chi array reduced chi-square values for levels 0.10,0.05,0.01
c       and for a number of degree of freedom up to 20
c
        DATA chi/2.706,2.303,2.084,1.945,1.847,1.774,1.717,1.670,1.632,
     &  1.599,1.570,1.546,1.524,1.505,1.487,1.471,1.457,1.444,1.432,
     &  1.421,3.841,2.996,2.605,2.372,2.214,2.099,2.010,1.938,1.880,
     &  1.831,1.789,1.752,1.720,1.692,1.666,1.644,1.623,1.604,1.586,
     &  1.571,6.635,4.605,3.780,3.319,3.017,2.802,2.639,2.511,2.407,
     &  2.321,2.248,2.185,2.130,2.082,2.039,2.000,1.965,1.934,1.905,
     &  1.878/
c
c       Results storage ?
c       -----------------
c
        pr=.true.                       ! enable record of results
        write(5,2)
2       format(' Do you want a report ? [default : Y] ? ',$)
        read(5,1010,end=800) reppr
        if ((reppr.eq.'n').or.(reppr.eq.'N')) pr=.false.
        if(pr) then
         open(4,name='SELECT.DAT',status='new')
         type *,' ** Your report is in the file SELECT.DAT **'
        endif
c
c       Set up
c       ------
c
3       ind=.true.                      ! assumed independent measurements
        check=.true.                    ! enable data checking
        sel=.true.                      ! enable selection of data
        next=.false.                    ! disable test at next level
        do 4 i=0,lmax
4        nvalid(i)=0
c
c       Ask for input values
c       --------------------
c
5       write(5,10) nmax
10      format(' How many measurements of the SAME unknown [max:',
     &  i2,' min:3] ? ',$)
        read(5,*,end=800) ndata
        if((ndata.gt.nmax).or.(ndata.lt.3)) go to 5
        do 40 i=1,ndata
         nord(i)=i                              ! rank of input data
15       write(5,20) i
20       format(' Data # ',i2,' enter value & standard deviation ',$)
         read(5,*,end=800) data(i,1),data(i,2)
         if(data(i,2).le.0.) go to 15           ! need standard deviation !!
40       error(i,i)=data(i,2)*data(i,2)         ! variance of input # i
c
c       Independent measurements ?
c       --------------------------
c
        type 50
50      format(' Are these data INDEPENDENT ? [default : Y] ? ',$)
        read(5,1010,end=800) repin
        if ((repin.eq.'n').or.(repin.eq.'N')) ind=.false.
        if(ind) then
         do 70 i=1,ndata-1
          do 70 j=i+1,ndata
           error(i,j)=0.
70         error(j,i)=error(i,j)
        else
         do 90 i=1,ndata-1
          do 90 j=i+1,ndata
           write(5,80) i,j
80         format(' Input covariance(',i2,',',i2,') : ',$)
           read(5,*,end=800) error(i,j)
90         error(j,i)=error(i,j)
        endif
c
c       Check data ?
c       ------------
c
        write(5,95)
95      format(' Do you want to check your data [default : Y] ? ',$)
        read(5,1010,end=800)repch
        if((repch.eq.'n').or.(repch.eq.'N')) check=.false.
c
c       Estimate weighted mean, standard deviation(mean) & reduced chi-square
c       ---------------------------------------------------------------------
c       on all data
c       -----------
c
        level=0
        nvalid(0)=ndata
        call dpmean(data,error,nvalid,level,wmean,sdmean,diff,nfree,chi2)
        if(check) then                          ! check or output ?
         continue
        else
         go to 500
        endif
c
c       Check data at levels .10,.05,.01
c       --------------------------------
c
        do 300 level=1,lmax
         if((level.eq.1).or.next) then          ! worth to continue checking ?
          continue
         else
          go to 500
         endif
c
c        start with all data
c        -------------------
c
         k=0
         call dpmean(data,error,nvalid,k,wmean,sdmean,diff,nfree,chi2)
         nvalid(level)=nvalid(0)
         wmean(level)=wmean(0)
         sdmean(level)=sdmean(0)
         nfree(level)=nfree(0)
         chi2(level)=chi2(0)
         next=.false.                           ! suppose last check
c
c        chi-square test
c        ---------------
c
100      if(nfree(level).eq.0) go to 300        ! no more data to exclude
         if(chi2(level).gt.chi(nfree(level),level)) then
          sel=.true.                            ! enable data selection
          next=.true.                           ! => must check at next level
         else
          sel=.false.
         endif
         if(sel) then                   ! Exit the data with greatest deviation
          if(level.gt.1) go to 200              ! data already organized
          nexit=1                               ! search data to be excluded
           do 150 i=2,nvalid(level)
            if(diff(i).gt.diff(nexit)) nexit=i  ! index of data to be excluded
150        continue
          if(nexit.ne.nvalid(level)) then
c
c          exchange last valid data with excluded one
c
           do 180 j=1,2
            ntemp=nord(nexit)                           ! exchange data ranks
            nord(nexit)=nord(nvalid(level))
            nord(nvalid(level))=ntemp
            temp=data(nexit,j)                          ! exchange values & sd
            data(nexit,j)=data(nvalid(level),j)
180         data(nvalid(level),j)=temp
           call permat(error,nmax,nexit,nvalid(level))  ! exchange covariances
          endif
c
c         test again without excluded data
c
200       nvalid(level)=nvalid(level)-1
          call dpmean(data,error,nvalid,level,wmean,sdmean,diff,
     &    nfree,chi2)
          go to 100
         endif
300     continue
c
c       Output results
c       --------------
c
500     write(5,505)
        if(pr) write(4,505)
505     format(//,3x,'Estimation of weighted mean',
     &  /,3x,'---------------------------',/)
        do 520 i=1,ndata
         write(5,510) nord(i),data(nord(i),1),data(nord(i),2)
         if(pr) write(4,510) nord(i),data(nord(i),1),data(nord(i),2)
510      format(' Data # ',i2,' Value : ',g,' Standard deviation : ',g)
520     continue
         if(ind) then
          continue
         else
          do 550 i=1,ndata-1
           do 550 j=i+1,ndata
            write(5,540) nord(i),nord(j),error(nord(i),nord(j))
            if(pr) write(4,540) nord(i),nord(j),error(nord(i),nord(j))
540         format(' Covariance(',i2,',',i2,') : ',e10.3)
550       continue
         endif
        do 600 level=0,lmax
         if(nvalid(level).eq.0) go to 700
         write(5,570) risk(level),nvalid(level)
         if(pr) write(4,570) risk(level),nvalid(level)
570      format(//,' Results at level ',f4.2,
     &   /,' ---------------------',//,
     &   ' Selected data : from 1 to ',i2)
         write(5,580) wmean(level),sdmean(level)
         if(pr) write(4,580) wmean(level),sdmean(level)
580      format(' Mean value : ',e12.5,
     &   /,' Standard deviation(mean) : ',e12.5)
         if(check) then
          continue
         else
          go to 700
         endif
600     continue
c
c       Once more
c       ---------
c
700     type *,' Start again. Input new data. Enter CTRL Z for exit'
        go to 3
c
c       Program exit
c       ------------
c
800     if(pr) then
         close(5)
         open(5,status='new')
         write(5,805)
805      format(10x,'Do you want your report to be printed ? ',
     &   '[default : N] ',$)
         read(5,1010,end=1000) pri
         if((pri.eq.'y').or.(pri.eq.'Y')) then
          close(4,dispose='print')
          type *,'Don''t forget your report !'
         endif
        endif
1000    stop
1010    format(a1)
        end




        SUBROUTINE DPMEAN(DATA,ERROR,NVALID,LEVEL,WMEAN,SDMEAN,DIFF,
     &  NFREE,CHI2)
c
c       LAST VERSION :  6-JAN-1988 11:52:43
c
c       purpose
c       -------
c       computes weighted mean, standard deviation of the mean and reduced
c       chi-square of dependant data with respect to the mean.
c
c       description of parameters
c       -------------------------
c               input :
c               -------
c       data ----> matrix (nmax x 2) containing data & standard deviation
c       error ---> matrix (nmax x nmax) containing variances & covariances
c       nvalid --> range of valid data in DATA and ERROR (nvalid < = nmax)
c       level ---> index for results storage in output variables
c
c               output :
c               --------
c       wmean ---> weighted mean computed on the first NVALID elements
c       sdmean --> standard deviation of the mean
c       diff ----> array (nmax) containing, for each data, the deviation from
c                  the mean
c       nfree ---> number of degree of freedom = NVALID-2
c                  (determination of mean AND standard deviation)
c       chi2 ----> reduced chi-square of data with respect to the mean
c
c       subroutines and function subprograms required
c       ---------------------------------------------
c       none
c
c       comments
c       --------
c       1) This program cannot handle more than 22 data pairs.
c       2) This program stores the output values in arrays ranging from 0 to 3.
c       3) Modify PARAMETER statements to increase these limits in BOTH DPMEAN
c          and SELECT programs.
c       4) This program can handle DEPENDANT data.
c
c       author :
c       --------
c       Pierre TREBBIA
c       US 41 : "Microscopie Electronique Analytique Quantitative"
c       Laboratoire de Physique des Solides, Bat. 510
c       Universite Paris-Sud, F91405 ORSAY Cedex
c       Phone : (33-1) 69 41 53 68
c
        PARAMETER (NMAX=22)
        PARAMETER (LMAX=3)
        double precision sumwgt,wmean,error,vmean,chi2,diff
        dimension data(nmax,2),error(nmax,nmax),wmean(0:lmax)
        dimension sdmean(0:lmax),nvalid(0:lmax),nfree(0:lmax)
        dimension chi2(0:lmax),diff(nmax)
c
c       Total weighting of data
c       -----------------------
c
100     sumwgt=0.
        do 110 i=1,nvalid(level)
110      sumwgt=sumwgt+1./error(i,i)
c
c       Estimate weighted mean
c       ----------------------
c
120     wmean(level)=0.
        do 150 i=1,nvalid(level)
150      wmean(level)=wmean(level)+data(i,1)/error(i,i)
        wmean(level)=wmean(level)/sumwgt
c
c       Estimate standard deviation of weighted mean
c       --------------------------------------------
c
        vmean=0.                ! variance of the mean
        do 190 i=1,nvalid(level)
         do 190 j=1,nvalid(level)
190      vmean=vmean+error(i,j)/(error(i,i)*error(j,j))
        xn=float(nvalid(level))
        vmean=(xn/(xn-1.))*vmean/(sumwgt*sumwgt)
        if(vmean.lt.0.) then
         type *,' Variance of the mean is negative !'
         call exit
        endif
        sdmean(level)=sqrt(vmean)
c
c       Compute chi-square
c       ------------------
c
        nfree(level)=nvalid(level)-2            ! number of degree of freedom
        if(nfree(level).eq.0) then              ! give up !
         chi2(level)=0.
         return
        endif
        chi2(level)=0.
        do 200 i=1,nvalid(level)
         diff(i)=abs(data(i,1)-wmean(level))    ! deviation from mean
200      chi2(level)=chi2(level)+diff(i)*diff(i)
        chi2(level)=chi2(level)/(float(nfree(level))*vmean)
300     return
        end
        SUBROUTINE PERMAT(XMAT,NORDER,N1,N2)
c
c       LAST VERSION : 24-JUL-1987 17:20:43
c
c       purpose
c       -------
c       Permutes column & row N2 with column & row N1 in the NORDER by NORDER
c       square matrix XMAT.
c
c       description of parameters
c       -------------------------
c               input :
c               -------
c       xmat ----> square matrix
c       norder --> order of the matrix
c       n1 ------> first selected index
c       n2 ------> second selected index
c
c               output :
c               --------
c       xmat ----> reorganized matrix
c
c       subroutines and function subprograms required
c       ---------------------------------------------
c       none
c
c       comments
c       --------
c       The matrix XMAT can be non symmetric.
c
c       author :
c       --------
c       Pierre TREBBIA
c       US 41 : "Microscopie Electronique Analytique Quantitative"
c       Laboratoire de Physique des Solides, Bat. 510
c       Universite Paris-Sud, F91405 ORSAY Cedex
c       Phone : (33-1) 69 41 53 68
c
        double precision xmat,term
        dimension xmat(norder,norder)
c
c       Reorganize matrix keeping diagonal terms in diagonal positions
c       --------------------------------------------------------------
c
        do 200 i=1,norder
         term=xmat(n2,i)
         xmat(n2,i)=xmat(n1,i)
200      xmat(n1,i)=term
        do 300 j=1,norder
         term=xmat(j,n1)
         xmat(j,n1)=xmat(j,n2)
300      xmat(j,n2)=term
        return
        end


                DEMONSTRATION OF PROGRAM SELECT

  Estimation of weighted mean
  ---------------------------

Data #  1 Value :   0.4249200     Standard deviation :   0.3677000E-02
Data #  2 Value :   0.4267000     Standard deviation :   0.4000000E-02
Data #  3 Value :   0.4325000     Standard deviation :   0.5000000E-02
Data #  4 Value :   0.4128000     Standard deviation :   0.7000000E-02
Data #  5 Value :   0.4518000     Standard deviation :   0.1200000E-01


Results at level 0.00
---------------------

Selected data : from 1 to  5
Mean value :  0.42666E+00
Standard deviation(mean) :  0.24765E-02


Results at level 0.10
---------------------

Selected data : from 1 to  2
Mean value :  0.42574E+00
Standard deviation(mean) :  0.38283E-02


Results at level 0.05
---------------------

Selected data : from 1 to  2
Mean value :  0.42574E+00
Standard deviation(mean) :  0.38283E-02


Results at level 0.01
---------------------

Selected data : from 1 to  3
Mean value :  0.42727E+00
Standard deviation(mean) :  0.29155E-02