モンテカルロ法でπを求めるのは,よく見かける(例)。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: