Create a gist now

Instantly share code, notes, and snippets.

An implementation of the Metropolis-Hastings algorithm and a Bayesian model in Microsoft QBasic 1.1 (but written in a more old-school BASIC style). More info at http://www.sumsar.net
0 CLS
10 REM SETTING UP THE VARIABLES AND DATA
20
30 REM NO OF PUPS PER DEN FROM A SAMPLE OF 16 WOLF DENS
40 DATA 5, 8, 7, 5, 3, 4, 3, 9, 5, 8, 5, 6, 5, 6, 4, 7
50
60 NDATA = 16
70 NSAMPLES = 1000
75 NBURNIN = 250
80 NPARAMS = 2
90 DIM SAMPLES!(NSAMPLES, NPARAMS)
100 DIM PARAMS!(NPARAMS)
110 DIM PROPOSAL!(NPARAMS)
120 REM INITIALIZING THE PARAMETERS TO 0.0
130 FOR J = 0 TO NPARAMS - 1
140 SAMPLES!(1, J) = 0.0
150 NEXT J
160
170 PRINT "RUNNING THE METROPOLIS-HASTINGS ALGORITHM"
180 GOSUB 240
190 PRINT "FINISHED! NOW A SUMMARY OF THE POSTERIOR DISTRIBUTION"
200 GOSUB 690
202 PRINT "NOW SHOWING TRACE PLOTS"
205 GOSUB 890
207 PRINT "WRITING THE SAMPLES TO DISK"
208 GOSUB 1200
210 PRINT "NOW EXITING. HAVE A GOOD DAY!"
220 GOTO 9999
230
240 REM THE METROPOLIS-HASTINGS ALGORITHM
250 FOR I = 1 TO NSAMPLES - 1
260 FOR J = 0 TO NPARAMS - 1
270 PARAMS!(J) = SAMPLES!(I - 1, J)
280 NEXT J
285 REM
290 GOSUB 520
300 CURRDENS! = LOGDENS!
310
320 FOR J = 0 TO NPARAMS - 1
330 PROPOSAL!(J) = SAMPLES!(I - 1, J) + RND - 0.5
340 PARAMS!(J) = PROPOSAL!(J)
350 NEXT J
360 REM
370 GOSUB 520
380 PROPDENS! = LOGDENS!
390
400 IF RND < EXP(PROPDENS! - CURRDENS!) THEN
410 FOR J = 0 TO NPARAMS - 1
420 SAMPLES!(I, J) = PROPOSAL!(J)
430 NEXT J
440 ELSE
450 FOR J = 0 TO NPARAMS - 1
460 SAMPLES!(I, J) = SAMPLES!(I - 1, J)
470 NEXT J
480 END IF
490 NEXT I
500 RETURN
510
520 REM THIS IS THE MODEL
530 MEAN! = PARAMS!(0)
540 SCALE! = EXP(PARAMS!(1))
550
560 REM USING A SLOPPY UNIFORM PRIOR OVER BOTH PARAMETERS
570 LOGPRIOR! = 0
580
590 LOGLIKE! = 0
600 RESTORE
610 FOR IDATA = 1 TO NDATA
620 READ X!
630 REM BELOW IS PROPORTIONAL TO THE LOG-LIKELIHOOD OF A LAPLACE DIST
640 LOGLIKE! = LOGLIKE! + (-ABS(X! - MEAN!) / SCALE!) - LOG(SCALE!)
650 NEXT IDATA
660 LOGDENS! = LOGPRIOR! + LOGLIKE!
670 RETURN
680
690 REM PRODUCING A SUMMARY OF THE POSTERIOR
700 FOR J = 0 TO NPARAMS - 1
710 SUM! = 0
720 FOR I = NBURNIN TO NSAMPLES - 1
730 SUM! = SUM + SAMPLES!(I, J)
740 NEXT I
750 POSTMEAN! = SUM! / (NSAMPLES - NBURNIN)
760 SUMSQUARE! = 0
770 FOR I = NBURNIN TO NSAMPLES - 1
780 SUMSQUARE! = SUMSQUARE! + (POSTMEAN! - SAMPLES!(I, J))^2
790 NEXT I
800 POSTSD! = SQR(SUMSQUARE! / (NSAMPLES - NBURNIN))
810 LOWCI! = POSTMEAN! - 1.96 * POSTSD!
820 HIGHCI! = POSTMEAN! + 1.96 * POSTSD!
830
840 PRINT "PARAMETER"; J + 1
850 PRINT " MEAN: "; POSTMEAN!; ", SD: "; POSTSD!
860 PRINT " 95% CI: ["; LOWCI!; ", "; HIGHCI!; "]"
870 NEXT J
875 RETURN
880
890 REM DRAWING TRACEPLOTS
900 PRINT "PRESS ANY KEY"
910 FOR J = 0 TO NPARAMS - 1
920 INPUT ;TEMP$
930 SCREEN 13
940 CLS
950 MIN! = 99999
960 MAX! = -99999
970 FOR I = 0 TO NSAMPLES - 1
980 IF SAMPLES!(I,J) > MAX! THEN
990 MAX! = SAMPLES!(I,J)
1000 END IF
1010 IF SAMPLES!(I,J) < MIN! THEN
1020 MIN! = SAMPLES!(I,J)
1030 END IF
1040 NEXT I
1050
1060 COL = 2
1070 FOR I = 0 TO NSAMPLES - 2
1080 IF I > NBURNIN THEN
1090 COL = 1
1100 END IF
1110 X1! = I / NSAMPLES * 320
1120 Y1! = 200 - (SAMPLES!(I, J) - MIN!) / (MAX! - MIN!) * 200
1130 X2! = (I + 1) / NSAMPLES * 320
1140 Y2! = 200 -(SAMPLES!(I + 1, J) - MIN!) / (MAX! - MIN!) * 200
1150 LINE (X1!, Y1!)-(X2!, Y2!), COL
1160 NEXT I
1170 NEXT J
1180 RETURN
1190
1200 REM WRING THE SAMPLES TO DISK
1210 OPEN "MCMC.CSV" FOR OUTPUT AS 1
1220 FOR I = 0 TO NSAMPLES - 1
1230 FOR J = 0 TO NPARAMS - 1
1240 PRINT #1, SAMPLES!(I, J);
1250 NEXT J
1260 PRINT #1,
1270 NEXT I
1280 CLOSE 1
1290 RETURN
1300
9999 REM EXIT THE PROGRAM
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment