モンテカルロ法で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: