Use Monte Carlo Simulation for Probability in C

Introduction

In this lab, we will explore how to use Monte Carlo simulation to estimate probabilities in C programming. We will start by defining a simple random experiment, such as a coin flip, and then run multiple trials to count the number of successful outcomes. Finally, we will calculate the estimated probability by dividing the number of successes by the total number of trials. This lab provides a practical introduction to the fundamental concepts of probability and combinatorics using C, which are essential skills for data analysis and decision-making.

Skills Graph

Define a Random Experiment

In this step, we'll explore how to define a random experiment using Monte Carlo simulation in C. A random experiment is a process with uncertain outcomes that can be simulated using probability techniques.

Understanding Random Experiments

Let's create a simple C program to demonstrate a basic random experiment: coin flipping simulation.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define NUM_TRIALS 1000

int main() {
    // Seed the random number generator
    srand(time(NULL));

    // Count of heads
    int heads_count = 0;

    // Simulate coin flips
    for (int i = 0; i < NUM_TRIALS; i++) {
        // Generate random number 0 or 1
        int flip = rand() % 2;

        // Count heads
        if (flip == 0) {
            heads_count++;
        }
    }

    // Calculate probability of heads
    double probability = (double)heads_count / NUM_TRIALS;

    printf("Coin Flip Experiment:\n");
    printf("Total Trials: %d\n", NUM_TRIALS);
    printf("Heads Count: %d\n", heads_count);
    printf("Estimated Probability of Heads: %.2f\n", probability);

    return 0;
}

Example output:

Coin Flip Experiment:
Total Trials: 1000
Heads Count: 502
Estimated Probability of Heads: 0.50

Key Concepts Explained

1. Random Number Generation:

   • srand(time(NULL)) seeds the random number generator
   • rand() % 2 generates either 0 or 1 with equal probability

2. Experiment Design:

   • We define a coin flip as our random experiment
   • Run multiple trials (1000 in this case)
   • Count the number of successful outcomes (heads)

3. Probability Estimation:

   • Probability = (Number of Successful Outcomes) / (Total Number of Trials)
   • In this case, we expect close to 0.5 probability for heads

Compile and Run the Program

## Create the source file
nano ~/project/coin_flip_experiment.c

## Compile the program
gcc ~/project/coin_flip_experiment.c -o ~/project/coin_flip_experiment

## Run the experiment
~/project/coin_flip_experiment

Example compilation and run output:

## Compilation
gcc ~/project/coin_flip_experiment.c -o ~/project/coin_flip_experiment

## Execution
~/project/coin_flip_experiment
Coin Flip Experiment:
Total Trials: 1000
Heads Count: 502
Estimated Probability of Heads: 0.50

Run Many Random Trials and Count Successes

In this step, we'll expand our Monte Carlo simulation by running multiple random trials and accurately counting successful outcomes. We'll demonstrate this through a more complex probability experiment: estimating π (pi) using random point generation.

Monte Carlo Method for π Estimation

We'll use a geometric approach to estimate π by generating random points within a square that contains a quarter circle.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define NUM_TRIALS 100000

int main() {
    // Seed random number generator
    srand(time(NULL));

    // Counters for total and inside points
    int total_points = 0;
    int points_inside_circle = 0;

    // Run multiple trials
    for (int i = 0; i < NUM_TRIALS; i++) {
        // Generate random x and y coordinates between 0 and 1
        double x = (double)rand() / RAND_MAX;
        double y = (double)rand() / RAND_MAX;

        // Check if point is inside quarter circle
        if (sqrt(x*x + y*y) <= 1.0) {
            points_inside_circle++;
        }
        total_points++;
    }

    // Estimate π
    double pi_estimate = 4.0 * points_inside_circle / total_points;

    printf("π Estimation Experiment:\n");
    printf("Total Points: %d\n", total_points);
    printf("Points Inside Circle: %d\n", points_inside_circle);
    printf("Estimated π: %.6f\n", pi_estimate);
    printf("Actual π:    %.6f\n", M_PI);

    return 0;
}

Compile and Run the Experiment

## Create the source file
nano ~/project/pi_estimation.c

## Compile the program (note the math library link)
gcc ~/project/pi_estimation.c -o ~/project/pi_estimation -lm

## Run the experiment
~/project/pi_estimation

Example output:

π Estimation Experiment:
Total Points: 100000
Points Inside Circle: 78540
Estimated π: 3.141600
Actual π:    3.141593

Key Concepts Explained

1. Multiple Trials:

   • We run a large number of random trials (100,000)
   • Each trial generates a random point in a 1x1 square

2. Counting Successes:

   • Track total points and points inside the quarter circle
   • Success is defined as a point falling within the circle

3. Probability Estimation:

   • Probability = (Successful Points) / (Total Points)
   • Multiply by 4 to estimate π due to quarter circle method

Important Techniques

• (double)rand() / RAND_MAX generates random decimal between 0 and 1
• sqrt(x*x + y*y) calculates distance from origin
• Large number of trials improves estimation accuracy

Estimate Probability = Successes/Trials

In this final step, we'll demonstrate how to calculate probability by analyzing the ratio of successful outcomes to total trials in a more practical scenario.

Dice Rolling Probability Experiment

We'll simulate rolling two dice and calculate the probability of getting a specific sum (e.g., sum of 7).

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define NUM_TRIALS 100000
#define TARGET_SUM 7

int roll_die() {
    return rand() % 6 + 1;
}

int main() {
    // Seed random number generator
    srand(time(NULL));

    // Counters for total rolls and successful rolls
    int total_rolls = 0;
    int successful_rolls = 0;

    // Run multiple trials
    for (int i = 0; i < NUM_TRIALS; i++) {
        // Roll two dice
        int die1 = roll_die();
        int die2 = roll_die();

        // Check if sum matches target
        if (die1 + die2 == TARGET_SUM) {
            successful_rolls++;
        }
        total_rolls++;
    }

    // Calculate probability
    double probability = (double)successful_rolls / total_rolls;

    printf("Dice Rolling Probability Experiment:\n");
    printf("Target Sum: %d\n", TARGET_SUM);
    printf("Total Rolls: %d\n", total_rolls);
    printf("Successful Rolls: %d\n", successful_rolls);
    printf("Estimated Probability: %.4f\n", probability);

    // Theoretical probability for comparison
    printf("Theoretical Probability: %.4f\n", 1.0/6);

    return 0;
}

Compile and Run the Experiment

## Create the source file
nano ~/project/dice_probability.c

## Compile the program
gcc ~/project/dice_probability.c -o ~/project/dice_probability

## Run the experiment
~/project/dice_probability

Example output:

Dice Rolling Probability Experiment:
Target Sum: 7
Total Rolls: 100000
Successful Rolls: 16644
Estimated Probability: 0.1664
Theoretical Probability: 0.1667

Key Concepts Explained

1. Probability Calculation:

   • Probability = (Number of Successful Outcomes) / (Total Number of Trials)
   • In this case: Successful Rolls / Total Rolls

2. Monte Carlo Simulation:

   • Large number of trials (100,000) provides accurate estimation
   • Simulated probability closely matches theoretical probability

3. Randomness and Accuracy:

   • srand(time(NULL)) ensures different random sequences
   • More trials increase estimation accuracy

Probability Interpretation

• Estimated probability (0.1664) is very close to theoretical probability (1/6 ≈ 0.1667)
• Demonstrates how Monte Carlo method can estimate probabilities

Summary

In this lab, we learned how to define a random experiment using Monte Carlo simulation in C. We created a simple coin flip simulation to demonstrate the key concepts. First, we seeded the random number generator and simulated coin flips, counting the number of successful outcomes (heads). Then, we estimated the probability of heads by dividing the number of successful outcomes by the total number of trials. The output showed the estimated probability was close to the expected value of 0.5 for a fair coin flip.