' ******************************************************************
' Statistical distribution
' ******************************************************************

#INCLUDE "fbmath.bi"
#INCLUDE "fbgfx.bi"

' ******************************************************************
' General declarations
' ******************************************************************

' Number of values
CONST N = 30

' Data
DIM AS DOUBLE X(1 TO N) = _
{141, 144, 146, 148, 149, 150, 150, 151, 153, 153, _
 153, 154, 155, 156, 156, 160, 160, 160, 163, 164, _
 164, 165, 166, 168, 168, 170, 172, 172, 176, 179}

' Mean and standard deviation
DIM SHARED AS DOUBLE M, S

' Significance level
CONST Alpha = 0.05

' ******************************************************************
' Subprograms
' ******************************************************************

SUB WriteResults(C()    AS StatClass, _
                 Calc() AS DOUBLE,    _
                 Ntot   AS INTEGER,   _
                 Khi2   AS DOUBLE,    _
                 G      AS DOUBLE,    _
                 DoF    AS INTEGER,   _
                 K2c    AS DOUBLE)
' ------------------------------------------------------------------
' Writes results to screen
' ------------------------------------------------------------------

  DIM AS STRING  Line1 = STRING$(28, "-")
  DIM AS DOUBLE  Sum
  DIM AS INTEGER I

  PRINT "Statistical distribution"

  PRINT Line1
  PRINT "    Inf    Sup      N  Ncalc"
  PRINT Line1

  Sum = 0
  FOR I = LBOUND(C) TO UBOUND(C)
    PRINT USING "#######"; C(I).Inf; C(I).Sup; C(I).N;
    PRINT USING "####.##"; Calc(I)
    Sum = Sum + Calc(I)
  NEXT I

  PRINT Line1
  PRINT USING "Total            #### ###.##"; Ntot; Sum
  PRINT Line1

  PRINT
  PRINT "Comparison with normal distribution:"
  PRINT

  PRINT USING "Pearson's Khi-2            = ####.####  (### DoF)"; Khi2; DoF
  PRINT USING "Woolf's G                  = ####.####  (### DoF)"; G; DoF
  PRINT USING "Critical value (p = #.###) = ####.####"; Alpha; K2c
END SUB

FUNCTION PltFunc(X AS DOUBLE) AS DOUBLE
' ------------------------------------------------------------------
' Function to be plotted (density of normal distribution)
' ------------------------------------------------------------------

  PltFunc = DNorm((X - M) / S) / S
END FUNCTION

SUB PlotGraph(C()   AS StatClass, _
              Xmin  AS DOUBLE,    _
              Xmax  AS DOUBLE,    _
              Xstep AS DOUBLE)
' ------------------------------------------------------------------
' Plots histogram and normal curve
' ------------------------------------------------------------------

  DIM AS DOUBLE  Ymin, Ymax, Ystep  ' Oy scale
  DIM AS INTEGER I, Lb, Ub

  IF NOT InitScreen THEN
    PRINT "Unable to set graphic mode!"
    EXIT SUB
  END IF

  Lb = LBOUND(C)
  Ub = UBOUND(C)

  DIM AS DOUBLE Y(Lb TO Ub)

  FOR I = Lb TO Ub
    Y(I) = C(I).D
  NEXT I

  AutoScale Y(), LinScale, Ymin, Ymax, Ystep

  if Ymin <> 0 then Ymin = 0
  Ymax = Ymax + Ystep

  InitGraph

  SetOxScale LinScale, Xmin, Xmax, Xstep
  SetOyScale LinScale, Ymin, Ymax, Ystep

  PlotOxAxis Ymin, 4
  PlotOyAxis Xmin, 4

  ' PlotGrid

  WriteTitles "Statistical Distribution", _
              "X", "Frequency Density"

  SetClipping

  FOR I = Lb TO Ub
    IF Y(I) > 0 THEN
      LINE (Xpixel(C(I).Inf), Ypixel(0)) - _
           (Xpixel(C(I).Sup), Ypixel(Y(I))),  _
           5, BF
    END IF
  NEXT I

  COLOR 12

  PlotFunc @PltFunc

  CloseGraph
END SUB

' ******************************************************************
' Main program
' ******************************************************************

DIM AS INTEGER I, Nc, Ntot, DoF
DIM AS DOUBLE  XMin, XMax, XStep, T, F0, F, Sum, Khi2, K2c, G

' Sort array X() if necessary
' QuickSort X(), 1, N

' Compute an appropriate interval for the set of values
Interval X(1), X(N), 5, 10, XMin, XMax, XStep

' Compute number of classes and dimension arrays
Nc = (Xmax - Xmin) / XStep

DIM AS StatClass C(1 TO Nc)
DIM AS INTEGER   Obs(1 TO Nc)
DIM AS DOUBLE    Calc(1 TO Nc)

' Compute distribution
Distrib X(), Xmin, XStep, Ntot, C()

' Compute mean and S.D.
M = Mean(X())
S = StDev(X(), M)

' Compute theoretical values
F0 = 0
FOR I = 1 TO Nc
  IF I = Nc THEN
    F = 1
  ELSE
    T = (C(I).Sup - M) / S
    F = FNorm(T)
  END IF
  Calc(I) = Ntot * (F - F0)
  Obs(I) = C(I).N
  F0 = F
NEXT I

' Perform Khi-2 and Woolf tests
Khi2_Conform Obs(), Calc(), 2, Khi2, DoF
Woolf_Conform Obs(), Calc(), 2, G, DoF

' Compute critical value
K2c = InvKhi2(DoF, 1 - Alpha)

' Print results
WriteResults C(), Calc(), Ntot, Khi2, G, DoF, K2c

' Plot histogram
PlotGraph C(), Xmin - Xstep, Xmax + Xstep, Xstep

SLEEP
