March 8, 2018

Simple sampling with Box-Muller transforms pt 1

I'm a technically uneducated idiot, so talking about math is a bit above my pay grade, but Box-Muller transform is a straightforward, computationally simple way to generate a distribution of points.

This is useful in many cases where you need to generate samples for simulations, so I think every programmer should have familiarity with it.

Lets take the example code from the wikipedia and convert it to standard C. This program will just print 60 x/y samples.

First, checkout the cool visualization of the transform on uniform inputs from the wikipedia! If you right click and view the image directly, you can hover on the crosses to see points there.

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

static const double epilson = DBL_MIN;
static const double pi2 = 2.0 * M_PI;

struct sample {
    double z0;
    double z1;
};

struct sample box_muller_sample() {
    double u1;
    double u2;
    struct sample ret;

    do {
        u1 = rand() * (1.0 / RAND_MAX);
        u2 = rand() * (1.0 / RAND_MAX);
    } while (u1 <= epilson);

    ret.z0 = sqrt(-2.0 * log(u1)) * cos(pi2 * u2);
    ret.z1 = sqrt(-2.0 * log(u1)) * sin(pi2 * u2);
    /* no sigma or mu */
    return ret;
}

int main(int argc, char **argv) {
    time_t curtime;
    srand(time(&amp;curtime));

    puts("Box-Muller samples:");
    for (int line = 0; line < 10; line++) {
        for (int i = 0; i < 6; i++) {
            struct sample bm = box_muller_sample();
            printf("%0.3f:%0.3f\t", bm.z0, bm.z1);
        }
        printf("\n");
    }
    return 0;
}

a sample output looks like this:

Box-Muller samples:
-0.559:0.786    -0.798:-0.640   -0.365:-0.140   1.417:-0.197    0.210:1.734     0.738:1.026
-2.334:0.079    -0.973:-0.634   -0.674:1.263    0.765:1.705     -0.560:0.860    -0.576:-0.044
-1.242:-1.705   1.705:-0.471    0.063:-0.079    0.732:-1.173    -0.242:-0.756   -1.798:0.853
-1.088:1.138    2.096:1.439     0.556:0.168     -0.105:0.132    -0.397:-0.728   0.090:0.109
-0.850:-0.034   0.930:0.108     -1.105:0.016    -0.592:-1.795   1.073:-0.227    -1.198:-0.517
-0.901:-1.978   -0.524:-0.090   0.821:-1.595    0.783:-0.051    0.580:-0.458    0.416:1.115
0.888:0.337     -0.074:-0.526   1.667:0.233     -1.339:0.075    0.563:-0.379    0.236:-0.229
0.222:0.472     -0.364:0.061    0.413:1.185     0.265:0.083     -1.120:1.978    0.939:0.807
-0.344:1.586    0.162:-0.274    -0.309:-0.098   -0.728:-0.441   -0.358:-0.310   1.324:-0.816
0.187:-0.160    0.035:0.083     1.386:-0.823    0.117:-1.897    0.475:-0.387    0.295:-0.382

In a future post I'll cover an even more efficient algorithm to generate uniform sampling, which is appropriate when a very large amount of random sampling numbers is needed. It's called Ziggurat


c sampling math


Previous post
Creating multiboot ELF kernels with FASM I've been meaning to write about this one for a while, since I hacked together some simple OSes last summer! Multiboot is magical. Multiboot is a
Next post
MSP430 assembly part 1: Blinking LED hello world Lets blink an LED on a MSP430 launchpad the old fashioned way! Assembly! Don't worry, MSP430 has a nice, small RISC instruction set. It's only true