PROGRAM random_walk
! simulation of a random walk in one dimension
DIM xcum(64),x2cum(64)
CALL initial(p,N,xcum(),x2cum(),ntrial)
FOR itrial = 1 to ntrial
    CALL walk(p,N,xcum(),x2cum())
NEXT itrial
CALL output(N,xcum(),x2cum(),ntrial)
END

SUB initial(p,N,xcum(),x2cum(),ntrial)
    RANDOMIZE
    INPUT prompt "maximum number of steps N = ": N
    LET p = 0.5
    LET ntrial = 1000             ! number of trials
    FOR istep = 1 to N
        LET xcum(istep) = 0       ! not necessary in True BASIC
        LET x2cum(istep) = 0
    NEXT istep
END SUB

SUB walk(p,N,xcum(),x2cum())
    LET x = 0       ! initial position of walker for each trial
    FOR istep = 1 to N
        IF rnd <= p then
           LET x = x + 1
        ELSE
           LET x = x - 1
        END IF
        ! collect data after every step
        CALL data(x,xcum(),x2cum(),istep)
    NEXT istep
END SUB

SUB data(x,xcum(),x2cum(),istep)
    LET xcum(istep) = xcum(istep) + x
    LET x2cum(istep) = x2cum(istep) + x*x
END SUB

SUB output(N,xcum(),x2cum(),ntrial)
    PRINT "# steps","<x>","<x^2>","<x^2> - <x>^2"
    PRINT
    FOR istep = 1 to N
        LET xbar = xcum(istep)/ntrial
        LET x2bar = x2cum(istep)/ntrial
        LET variance = x2bar - xbar*xbar
        PRINT istep,xbar,x2bar,variance
    NEXT istep
END SUB