Club Amiga de Montreal - CAM/CAM_CD_1.iso/files/521.lha/Chi-Square/BASICcode

Text File | 1991-06-10 | 9.4 KB | 425 lines

WINDOW 1,"Chi-Square Test",(0,0)-(640,200),31+256

DEFDBL D,a,b,e-g,p,s,t,x
DEFINT i-k,c,n,r

LIBRARY "requester.library"
  DECLARE FUNCTION DoFileIOWindow&() LIBRARY
  DECLARE FUNCTION GetFileIO&() LIBRARY
  DECLARE FUNCTION AutoFileMessage() LIBRARY

' subroutine to call file requestor from the requester.library.

'      ireqx,ireqy x and y coordinates for the requestor
'      WindowTitle$ a title displayed in the requestor's title bar
'      Values returned by the subroutine:
'           PathName$    full path for selected file
'           FileName$    file name only
'           VolName$     the volume name or drive id
'           DrawerName$  drawer name or subdirectory
'           DiscOrDir%   returns 1 if selection is a drive or directory only
'           NoExist%     returns 1 if file does not exist
'           flag%        returns 0 if request fails and 1 if successful
'           Reason$      short message explaining failure
'
SUB requestor(ireqx,ireqy,WindowTitle$,PathName$,FileName$,VolName$,DrawerName$,DiscOrDir%,NoExist%,flag%,reason$)
    LOCAL i,value
    DEFINT i,value
    DiscOrDir%=0 : NoExist%=0 : flag%=0
    PathName$="" : FileName$="" : VolName$="" : DrawerName$=""
    reason$=""
    PathBuffer$=SPACE$(202)
    ExtPtrBuff$=SPACE$(24)
    FileIO&=GetFileIO&()
    IF FileIO&=0 THEN
        CALL ReleaseFileIO&(FileIO&)
        flag%=0
        reason$="FileIO Failed"
        EXIT SUB
    END IF
    POKEL FileIO&+248,SADD(PathBuffer$)
    POKEL FileIO&+222,SADD(ExtPtrBuff$)
    POKEL FileIO&+244,SADD(WindowTitle$+CHR$(0))
    POKE FileIO&+261,1
    POKE FileIO&+262,1
    POKE FileIO&+263,0
    POKEW FileIO&+232,ireqx
    POKEW FileIO&+234,ireqy
    POKE FileIO&+1,4
    extsize&=LEN(ext$)
    POKE FileIO&+1,128
    result&=DoFileIOWindow&(FileIO&,0&)
    IF result&=-1 THEN
        CALL ResetTitle(FileIO,WINDOW(7))
        CALL ReleaseFileIO(FileIO)
        flag%=0
        reason$="Cancel"
        EXIT SUB
    END IF
    IF result&=0 THEN
        result&=AutoFileMessage&(0&,WINDOW(7))
        INPUT "Type path: ";PathName$
        CALL ParseString&(FileIO&,SADD(PathName$+CHR$(0)))
        GOTO CopyFN
    END IF

    pathname$ = ""
        FOR i = 0 TO 202
            value = PEEK(SADD(PathBuffer$)+i)
            IF value = 0 THEN EXIT FOR
            char$ = CHR$(value)
            pathname$ = pathname$+char$
        NEXT i

    CopyFN: 'Copy out the Filename to Filename$.
        FileName$ = ""
        FOR i = 0 TO 30
            value = PEEK(FileIO&+2+i)
            IF value = 0 THEN EXIT FOR
            char$ = CHR$(value)
            FileName$ = FileName$+char$
        NEXT i

    CopyDrawer: 'Copy out all the drawers
        DrawerName$ = ""
        FOR i = 0 TO 132
            value = PEEK(FileIO&+32+i)
            IF value = 0 THEN EXIT FOR
            char$ = CHR$(value)
            DrawerName$ = DrawerName$+char$
        NEXT i

    CopyVol: 'Copy out the diskname
        VolName$ = ""
        FOR i = 0 TO 30
            value = PEEK(FileIO&+164+i)
            IF value = 0 THEN EXIT FOR
            char$ = CHR$(value)
            VolName$ = VolName$+char$
        NEXT i

    extsize&=LEN(FileName$)
    IF extsize&=0 THEN DiscOrDir%=1
    IF PEEKL(FileIO&+240)=0 THEN NoExist%=1
    CALL ResetTitle&(FileIO&,WINDOW(7))
    CALL ReleaseFileIO&(FileIO&)
    flag%=1
END SUB

'subroutine to print message at line linenum1%, then wait for key press
'key press message is printed at line linenum2%

SUB message (VAL linenum1%, VAL linenum2%, VAL m$)
    LOCATE linenum1%,INT((80-LEN(m$))/2) : ?m$;
    LOCATE linenum2%,27,0 : ?"Press any key to continue";
    DO
        IF INKEY$<>"" THEN EXIT LOOP
        SLEEP
    LOOP
    CLS : LOCATE 1,1,1
END SUB

'subroutine to ask if user wants to continue
'returns flag%=0 if no, flag%=1 if yes

SUB continue (flag%)
    LOCAL t$ : t$="" : flag% = 0
    WHILE t$<>"YES" AND t$<>"NO" AND t$<>"Y" AND t$<>"N"
        CLS : LOCATE 12,32 : INPUT "Try again? ",t$
        t$=UCASE$(t$)
        IF t$="YES" OR t$="Y" THEN
            flag%=1 : EXIT SUB
        END IF
        IF t$="NO" OR t$="N" THEN
            flag%=0 : EXIT SUB
        END IF
    WEND
END SUB

'ln(Gamma) function

FUNCTION gammln(VAL xx)
    LOCAL stp,x,tmp,ser
    stp=2.50662827465
    x=xx-1.0
    tmp=x+5.5
    tmp=(x+0.5)*LOG(tmp)-tmp
    ser=1.0+76.18009173/(x+1.0)-86.50532033/(x+2.0)+24.01409822/(x+3.0)
    ser=ser-1.231739516/(x+4.0)+0.120858003e-2/(x+5.0)
    ser=ser-0.536382e-5/(x+6.0)
    gammln=tmp+LOG(stp*ser)
END FUNCTION

'incomplete Gamma function by series representation

SUB gser (VAL a,VAL x,gamser,gln)
    LOCAL itmax,eps,n,sum,del,ap
    itmax=100
    eps=3.0e-7
    gln=gammln(a)
    IF x<=0.0 THEN
        IF x<0.0 THEN ? "**note** gser: x<0"
        gamser=0.0
    ELSE
        ap=a
        sum=1.0/a
        del=sum
        FOR n=1 TO itmax
            ap=ap+1.0
            del=del*x/ap
            sum=sum+del
            IF ABS(del)<ABS(sum)*eps THEN GOTO gserlabel
        NEXT n
        ? "**note** gser: a too large, itmax too small"
gserlabel:    gamser=sum*EXP(-x+a*LOG(x)-gln)
    END IF
END SUB

'incomplete gamma function by continued fraction representation

SUB gcf(VAL a,VAL x, gammcf,gln)
    LOCAL itmax,eps,n,gold,g,fac,b1,b0,anf,ana,an,a1,a0
    itmax=100
    eps=3.0e-7
    gln=gammln(a)
    gold=0.0
    a0=1.0
    a1=x
    b0=0.0
    b1=1.0
    fac=1.0
    FOR n=1 TO itmax
        an=1.0*n
        ana=an-a
        a0=(a1+a0*ana)*fac
        b0=(b1+b0*ana)*fac
        anf=an*fac
        a1=x*a0+anf*a1
        b1=x*b0+anf*b1
        IF a1<>0.0 THEN
            fac=1.0/a1
            g=b1*fac
            IF ABS((g-gold)/g)<eps THEN GOTO gcflabel
            gold=g
        END IF
    NEXT n
    ? "**note** gcf: a too large, itmax too small"
gcflabel:  gammcf=EXP(-x+a*LOG(x)-gln)*g
END SUB

'incomplete Gamma function Q(a,x)=1-P(a,x)

FUNCTION gammq(VAL a,VAL x)
    LOCAL gamser,gammcf,gln
    IF (x<0.0) OR (a<=0.0) THEN
        ? "**OOPS!** gammq: invalid arguments"
    END IF
    IF x<a+1.0 THEN
        gser a,x,gamser,gln
        gammq=1.0-gamser
    ELSE
        gcf a,x,gammcf,gln
        gammq=gammcf
    END IF
END FUNCTION

'print instructions

CLS
? : ? : ?
?"This program will allow you to select a data file, and then calculate"
?"  the Chi-square parameter and p value from that file.  The data file"
?"   MUST be of a certain format."
?
?"  WARNING: This program will NOT accept files that do not conform to"
?"           the expected format!"
?
?"The format is the number of columns (MUST be an integer), followed"
?"  by the number of rows (MUST be an integer), followed by the data"
?"  table itself.  For example, an acceptable 3 x 2 table file would"
?"  be:"
?
?"          3   2"
?"          24  11  15"
?"          12  14   3"
message 23,23,""
?
?"You may use your favorite word processor, spreadsheet, or whatever to"
?"  generate the data file.  However, the file MUST be in ASCII.  This"
?"  is actually pretty simple--it usually just means selecting the"
?"  ""general"" option when saving the file from your word processor, or"
?"  the ""print to file"" option when saving the file from your"
?"  spreadsheet."
?
?"Note:  the program does not care about spaces in a file.  It DOES however"
?"  care about carriage returns, and will treat them as zeros.  IF YOU"
?"  HAVE EXTRA CARRIAGE RETURNS IN YOUR DATA FILE, DELETE THEM PRIOR TO"
?"  RUNNING THIS PROGRAM!  Only put carriage returns at the end of a"
?"  line of data, and no extra lines."
message 19,23,"Happy Analyzing!     C. Niederberger 2/91"

start:    CLS

CALL requestor (10,10,"Select Data File",Name$,f$,v$,d$,dod%,ne%,f%,r$)

IF dod%=1 OR ne%=1 OR f%=0 THEN
    CLS : message 11,13,"File lookup failed: "+r$
    continue flg%
    IF flg%=0 THEN GOTO finish ELSE GOTO start
END IF

OPEN Name$ FOR INPUT AS #1

'load data array with data from file

INPUT #1,Dnc    'nc = Dnc = number of columns

nc=INT(Dnc)
IF Dnc<>nc THEN
    CLOSE #1
    CLS : message 11,13,"number of columns is not an integer!"
    continue flg%
    IF flg%=0 THEN
        GOTO finish
    ELSE
        ERASE DObs, DExd, DSumRow, DSumCol
        GOTO start
    END IF
END IF

INPUT #1,Dnr    'nr = Dnr = number of rows

nr=INT(Dnr)
IF Dnr<>nr THEN
    CLOSE #1
    CLS : message 11,13,"number of rows is not an integer!"
    continue flg%
    IF flg%=0 THEN
        GOTO finish
    ELSE
        ERASE DObs, DExd, DSumRow, DSumCol
        GOTO start
    END IF
END IF

DIM DObs(nc,nr)
DIM DExd(nc,nr)
DIM DSumRow(nr)
DIM DSumCol(nc)

DSum=0
FOR r=1 TO nr
    DSumRow(r)=0
    FOR c=1 TO nc
    IF EOF(1) THEN
        CLOSE #1
        CLS : message 11,13,"Not enough data in table!"
        continue flg%
        IF flg%=0 THEN
            GOTO finish
        ELSE
            ERASE DObs, DExd, DSumRow, DSumCol
            GOTO start
        END IF
    END IF
        INPUT #1,DObs(c,r)
        DSumRow(r)=DSumRow(r)+DObs(c,r)
        DSum=DSum+DObs(c,r)
    NEXT c
NEXT r

IF NOT EOF(1) THEN
    CLOSE #1
    CLS : message 11,13,"Too much data in table!"
    continue flg%
    IF flg%=0 THEN
        GOTO finish
    ELSE
        ERASE DObs, DExd, DSumRow, DSumCol
        GOTO start
    END IF
END IF

CLOSE #1

'calculate column sums

FOR c=1 TO nc
    DSumCol(c)=0
    FOR r=1 TO nr
        DSumCol(c)=DSumCol(c)+DObs(c,r)
    NEXT r
NEXT c

'calculate expected table and Chi-Square

DChiSqr=0
FOR c=1 TO nc
    FOR r=1 TO nr
        DExd(c,r)=DSumRow(r)*DSumCol(c)/DSum
        DChiSqr=DChiSqr+(((DObs(c,r)-DExd(c,r))^2)/DExd(c,r))
    NEXT r
NEXT c

'print Chi-square results
'calculate and print significance level

? : ? : ? : ? : ? : ?
ndf=(nr-1)*(nc-1)
? "     For"ndf"degree(s) of freedom, Chi-Square ="DChiSqr
?
prob=gammq(0.5*ndf,0.5*DChiSqr)
IF prob<1E-15 THEN
    ?"     significance level p="prob
    ?
    ?"     !! significance level p <.000000000000001 !!"
    GOTO waituser
END IF
nzeros=INT(ABS(LOG10(prob))) : ans$="."
nmant=INT(10^(nzeros+1)*prob+1) : nmant$=RIGHT$(STR$(nmant),1)
IF nmant$="0" THEN
    nmant$="1"
    nzeros=nzeros-1
END IF
IF prob>=1.0 THEN
    ? "     **OOPS!** significance level >= 1.0"
ELSE
    ? "     significance level p ="prob
    IF nzeros=0 THEN
        ans$=ans$+nmant$
    ELSE
        FOR i=1 TO nzeros
            ans$=ans$+"0"
        NEXT i
        ans$=ans$+nmant$
    END IF
    IF prob<=.9 THEN
        ? "                     or p < "ans$
    ELSE
        ? "                     or p < 1.0"
    END IF
END IF

'after results display, wait for user entry

waituser:
LOCATE 15,27,0 : ?"Press any key to continue"
DO
    IF INKEY$<>"" THEN EXIT LOOP
    SLEEP
LOOP

continue flg%
IF flg%=0 THEN
    GOTO finish
ELSE
    ERASE DObs, DExd, DSumRow, DSumCol
    GOTO start
END IF

finish: WINDOW CLOSE 1
    STOP