Interactive PowerBasic Forum

Here's an updated version of the Metropolis-Hastings algorithm for estimating pi that includes graphics:

' Metropolis-Hastings Algorithm for Pi Estimation with Graphics
' Oxygen Basic, ConsoleG, OpenGL

#compact
%filename "metropolis_pi_graphics.exe"
'uses RTL64
% Title "Metropolis-Hastings Algorithm for Pi Estimation with Graphics"

'% WindowStyle WS_OVERLAPPEDWINDOW
'% Animated
'% ScaleUp
% PlaceCentral
% AnchorCentral

% shaders

uses consoleG

BeginScript

const int NUM_STEPS = 1000000;
const float RADIUS = 1.0;
const float PROPOSAL_SD = 0.1;

int acceptCount = 0;

procedure metropolisHastings() {
    float x = 0.0;
    float y = 0.0;
    float proposalX, proposalY, distSq, proposalDistSq, acceptProb;

    for (int i = 0; i < NUM_STEPS; i++) {
        // Generate a proposal point
        proposalX = x + rndn() * PROPOSAL_SD;
        proposalY = y + rndn() * PROPOSAL_SD;
        proposalDistSq = proposalX*proposalX + proposalY*proposalY;

        // Calculate the acceptance probability
        distSq = x*x + y*y;
        if (distSq > RADIUS*RADIUS && proposalDistSq <= RADIUS*RADIUS) {
            acceptProb = 1.0;
        } else if (distSq <= RADIUS*RADIUS && proposalDistSq > RADIUS*RADIUS) {
            acceptProb = (RADIUS*RADIUS / proposalDistSq);
        } else {
            acceptProb = 1.0;
        }

        // Accept or reject the proposal point
        if (rnd < acceptProb) {
            x = proposalX;
            y = proposalY;
            acceptCount++;
        }

        // Draw the current point
        pushstate();
        translate(x, y, 0);
        scale(0.01, 0.01, 0.01);
        color(1, 0, 0);
        go sphere();
        popstate();
    }
}

procedure drawCircle() {
    glBegin(GL_LINE_LOOP);
    for (int i = 0; i <= 360; i++) {
        float angle = i * PI / 180.0;
        float x = RADIUS * cos(angle);
        float y = RADIUS * sin(angle);
        glVertex2f(x, y);
    }
    glEnd();
}

procedure main() {
    cls();
    shading();

    drawCircle();

    metropolisHastings();

    print("Pi estimate: ");
    println(4.0 * acceptCount / NUM_STEPS);

    waitkey();
}

EndScript
In this version, we've added a `drawCircle` procedure that draws a circle with radius `RADIUS` centered at the origin. We've also added code to draw a small red sphere at each accepted point in the `metropolisHastings` procedure.

When you run this program, you should see a circle centered at the origin and a cloud of red points inside the circle. The cloud of points represents the accepted proposals in the Metropolis-Hastings algorithm, and should be roughly uniformly distributed inside the circle. The estimate of pi is calculated as `4.0 * acceptCount / NUM_STEPS`, where `acceptCount` is the number of accepted proposals and `NUM_STEPS` is the total number of steps in the random walk.

You can experiment with different values of `NUM_STEPS` and `PROPOSAL_SD` to see how they affect the accuracy of the estimate and the appearance of the graph. Note that larger values of `NUM_STEPS` will result in a more accurate estimate, but will also take longer to run.