Q-BASIC Program For 2-State Gibbs

CSU	Hayward
	Statistics	Department

Q-Basic Program For 2-State Gibbs Sampler

Notes: This Q-Basic computer program is similar to the S-Plus script used to produce the simulations in Section 4 of Suess, Fraser, and Trumbo: "Elementary uses of the Gibbs Sampler: Applications to Medical Screening Tests," in STATS #27, Winter 2000. This program is referenced in Session 3 of the Classroom Resources written in support of the STATS article.

The main differences between this Q-Basic program and the S-Plus script are that the use of uniformly distributed pseudorandom numbers in the simulation is more directly shown here, the algorithm may be easier to understand for those not familiar with S-Plus programming, no graphic displays are created, some intermediate results are shown in numerical form, and the program may be run with acceptable speed even on older DOS computers.

Instructions: Cut the code displayed between the next two horizontal rules and paste it into the program editor of Q-Basic.

'Gibbs Sampling -- State Space = {0, 1}
'Q-BASIC program by Bruce E. Trumbo
'Copyright 2000 Bruce E. Trumbo. All rights reserved.

CLS
RANDOMIZE TIMER
eta = .99
theta = .97
gamma = .4024
delta = .9998

Outer = 100      'Cannot exceed 32000
Inner = 1000     'Cannot exceed 32000

'Use of inner and outer loops permits total run M2 > 32000
'Also, intermediate results printed after each inner loop

M2 = Outer * Inner   'Total run length
M1 = M2 / 2          'Burn in set at halfway through

Sum = 0           'Total number runs with D = 1
Steady = 0        'Runs after burn in with D = 1
D = 0             'D(1) = 0 is starting value

FOR i = 1 TO Outer
        FOR j = 1 TO Inner
                m = (i - 1) * Inner + j

                R = RND
                T = 0
                IF D = 1 AND R < eta THEN T = 1
                IF D = 0 AND R > theta THEN T = 1

                R = RND
                D = 0
                IF T = 1 AND R < gamma THEN D = 1
                IF T = 0 AND R > delta THEN D = 1

                Sum = Sum + D
                IF m > M1 THEN Steady = Steady + D
        NEXT j
        PRINT USING "#######: "; m;
        PRINT USING "% all =##.##"; 100 * Sum / m;
        IF m > M1 THEN
                PRINT USING ",  % steady =##.##"; 100 * Steady / (m - M1)
        ELSE
                PRINT
        END IF
NEXT i

END