'**********************************************************
'* Calculate the mean and first third moments of a set of *
'* data Y(i)                                              *
'* ------------------------------------------------------ *
'* Description:                                           *
'*                                                        *
'* Input data:                                            *
'*              ndata:  Integer size of data set          *
'*              Y:      Real Vector of ndata ordinates    *
'* Output data:                                           *
'*              S1:     ndata-Mean of ndata values        *
'*              S2:     Sum (Y(i)-S1)^2 for i=1..ndata    *
'*              S3:     Sum (Y(i)-S1)^3 for i=1..ndata    *
'*              S4:     Sum (Y(i)-S1)^4 for i=1..ndata    *
'* ------------------------------------------------------ *
'* SAMPLE RUN:                                            *
'*  Number of data: 5                                     *
'*  1: 12                                                 *
'*  2: 9                                                  *
'*  3: 7                                                  *
'*  4: 15                                                 *
'*  5: 6                                                  *
'*                                                        *
'*  S1=  9.8                                              *
'*  S2=  54.8                                             *
'*  S3=  73.92001                                         *
'*  S4=  1024.976                                         *
'*                                                        *
'*  Error code: 0                                         *
'*                                                        *
'* ------------------------------------------------------ *
'* Ref.: "Journal of Applied Statistics (1972) vol.21,    *
'*        page 226".                                      *
'*                                                        *
'*                    Basic Release By J-P Moreau, Paris. *
'*                             (www.jpmoreau.fr)          *
'**********************************************************
'PROGRAM TMOMNTS;
DEFINT I-N

  NMAX = 25

  CLS
  PRINT
  INPUT "  Number of data: ", ndata
  DIM Y(ndata)
  FOR i = 1 TO ndata
    PRINT " "; i; ": "; : INPUT "", Y(i)
  NEXT i

  FOR i = 1 TO ndata
    X = Y(i): K = 4: N = i - 1: GOSUB 1000
    'call MOMNTS(X,K,N,S1,S2,S3,S4,IFAULT)
  NEXT i

  PRINT
  PRINT "  S1= "; S1
  PRINT "  S2= "; S2
  PRINT "  S3= "; S3
  PRINT "  S4= "; S4
  PRINT
  PRINT "  Error code: "; IFAULT
  PRINT

END



1000 'Subroutine MOMNTS(X, K, N, S1, S2, S3, S4, IFAULT)
' ------------------------------------------------------------------
'        ALGORITHM AS 52  APPL. STATIST. (1972) VOL.21, P.226

'        ADDS A NEW VALUE, X, WHEN CALCULATING A MEAN, AND SUMS
'        OF POWERS OF DEVIATIONS. N IS THE CURRENT NUMBER OF
'        OBSERVATIONS, AND MUST BE SET TO ZERO BEFORE FIRST ENTRY
' ------------------------------------------------------------------
'     AN, AN1, DX, DX2, ZERO, ONE, TWO, THREE, FOUR, SIX: REAL

      ZERO = 0!: ONE = 1!: TWO = 2!: THREE = 3!: FOUR = 4!: SIX = 6!

      IF K > 0 AND K < 5 AND N >= 0 THEN GOTO 10
      IFAULT = 1
      RETURN
10    IFAULT = 0
      N = N + 1
      IF N > 1 THEN GOTO 20

'     FIRST ENTRY, SO INITIALISE

      S1 = X
      S2 = ZERO
      S3 = ZERO
      S4 = ZERO
      RETURN

'     SUBSEQUENT ENTRY, SO UPDATE

20    AN = N
      AN1 = AN - ONE
      DX = (X - S1) / AN
      DX2 = DX * DX
      IF K = 1 THEN
	GOTO 60      'calculate only mean S1
      ELSEIF K = 2 THEN
	GOTO 50      'calculate S1 and S2
      ELSEIF K = 3 THEN
	GOTO 40      'calculate S1, S2 and S3
      ELSEIF K = 4 THEN
	GOTO 30      'calculate S1..S4
      END IF
30    S4 = S4 - DX * (FOUR * S3 - DX * (SIX * S2 + AN1 * (ONE + AN1 * AN1 * AN1) * DX2))
40    S3 = S3 - DX * (THREE * S2 - AN * AN1 * (AN - TWO) * DX2)
50    S2 = S2 + AN * AN1 * DX2
60    S1 = S1 + DX

RETURN

'end of file momnts.bas