モンテカルロ法でeを求める

モンテカルロ法でπを求めるのは,よく見かける()。eは [0,1) 一様分布の乱数の和が1を超える個数の期待値として求められる(Twitterで教えていただいた証明例:一様分布の和の平均到達時間 ―eを近似する確率シミュレーション―18th Putnam 1958 Problem A3)。

Pythonで書いてみよう:

import random

def ne():
    s = 0
    i = 0
    while s < 1:
        s += random.random()
        i += 1
    return i


n = 1000000
t = 0
for _ in range(n):
    t += ne()
print(t / n)

ほぼ 2.718 となる。


M1 Mac mini では n = 100000000 で 29.6s だった(arm64 ネイティブの Python 3.8.8)。

せっかくなので Numba を使ってみる。

import random
from numba import jit

@jit(nopython=True)
def ne():
    s = 0
    i = 0
    while s < 1:
        s += random.random()
        i += 1
    return i

def main():
    n = 100000000
    t = 0
    for _ in range(n):
        t += ne()
    print(t / n)

%time main()

8.65s になった。

メインループもJIT化してみる:

@jit(nopython=True)
def main():

2.8s になった。

一つの関数にしてJIT化する:

@jit(nopython=True)
def main():
    def ne():
        s = 0
        i = 0
        while s < 1:
            s += random.random()
            i += 1
        return i
    n = 100000000
    t = 0
    for _ in range(n):
        t += ne()
    print(t / n)

%time main()

2.23s になった。


Last modified: