!********************************************************** !* 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.800000 * !* S2= 54.80000 * !* S3= 73.92002 * !* S4= 1024.976 * !* * !* Error code: 0 * !* * !* ------------------------------------------------------ * !* Ref.: "Journal of Applied Statistics (1972) vol.21, * !* page 226". * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !********************************************************** PROGRAM TMOMNTS integer error,i,K,ndata real, pointer :: Y(:) real S1,S2,S3,S4 print *,' ' write(*,10,advance='no'); read *, ndata allocate(Y(1:ndata),stat=ialloc) do i=1, ndata write(*,20,advance='no') i read *, Y(i) end do do i=1, ndata call MOMNTS(Y(i),4,i-1,S1,S2,S3,S4,error) end do print *,' ' print *,' S1=', S1 print *,' S2=', S2 print *,' S3=', S3 print *,' S4=', S4 print *,' ' print *,' Error code:', error print *,' ' 10 format(' Number of data: ') 20 format(' ',i2,': ') END 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 REAL X, S1, S2, S3, S4, AN, AN1, DX, DX2, ZERO, ONE, TWO, & THREE, FOUR, SIX DATA ZERO, ONE, TWO, THREE, FOUR, SIX /0.0, 1.0, 2.0, 3.0, 4.0, 6.0/ IF (K .GT. 0 .AND. K .LT. 5 .AND. N .GE. 0) GOTO 10 IFAULT = 1 RETURN 10 IFAULT = 0 N = N + 1 IF (N .GT. 1) 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 GOTO (60, 50, 40, 30), K 30 S4 = S4 - DX * (FOUR * S3 - DX * (SIX * S2 + AN1 * & (ONE + AN1 ** 3) * DX2)) 40 S3 = S3 - DX * (THREE * S2 - AN * AN1 * (AN - TWO) * DX2) 50 S2 = S2 + AN * AN1 * DX2 60 S1 = S1 + DX RETURN END !end of file momnts.f90