user2770287
user2770287

Reputation: 153

Generating numbers from normal distribution in Python

I'm trying to test the speed of generating numbers from normal distribution by using Box–Muller transform against Marsaglia polar method. It is said that Marsaglia polar method is suppose to be faster than Box–Muller transform because it does not need to compute sin and cos. However, when I code this in Python, this is not true. Can someone verify this or explain to me why this is happening?

def marsaglia_polar():
    while True:
        x = (random.random() * 2) - 1
        y = (random.random() * 2) - 1
        s = x * x + y * y
        if s < 1:
            t = math.sqrt((-2) * math.log(s)/s)
            return x * t, y * t

def box_muller():
    u1 = random.random()
    u2 = random.random()

    t = math.sqrt((-2) * math.log(u1))
    v = 2 * math.pi * u2

    return t * math.cos(v), t * math.sin(v)

Upvotes: 5

Views: 1978

Answers (1)

Joshua Grigonis
Joshua Grigonis

Reputation: 768

For "fun", I wrote it up in go. The box_muller function is faster there as well. Also, it's about 10 times faster than the python version.

package main

import (
    "fmt"
    "math"
    "math/rand"
    "time"
)

func main() {
    rand.Seed(time.Now().UnixNano())
    now := time.Now()
    for i := 0; i < 1000000; i++ {
        marsaglia_polar()
    }
    fmt.Println("marsaglia_polar duration = ", time.Since(now))
    now = time.Now()
    for i := 0; i < 1000000; i++ {
        box_muller()
    }
    fmt.Println("box_muller duration      = ", time.Since(now))
}

func marsaglia_polar() (float64, float64) {
    for {
        x := random() * 2 - 1;
        y := random() * 2 - 1;
        s := x * x + y * y;
        if s < 1 {
            t := math.Sqrt((-2) * math.Log(s)/s);
            return x * t, y * t
        }
    }
}

func box_muller() (float64, float64) {
    u1 := random()
    u2 := random()
    t := math.Sqrt((-2) * math.Log(u1))
    v := 2 * math.Pi * u2
    return t * math.Cos(v), t * math.Sin(v)
}

func random() float64 {
    return rand.Float64()
}

Output:

marsaglia_polar duration =  104.308126ms
box_muller duration      =  88.365933ms

Upvotes: 3

Related Questions