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

MCMCDEMO.BAS

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

Thanks for magnificent example!