Python
シミュレーション
数学
最適化
組合せ最適化
More than 3 years have passed since last update.


お札と小銭の分布

もし、必ず財布の重量が最小になるように支払っていたら、財布の中身はどのような分布になるでしょうか?


シミュレーション

シミュレーションで確かめてみましょう。


  • 1万円は無限にあるものとします。

  • 1万回、商品を買って分布をみます。(最初の100回は無視します)

  • 商品の金額は、100円~9999円とし、ベンフォードの法則に従う分布とします。

  • 乱数の発生(rand_from_prob)は、Walker's Alias Methodを用います。

  • 支払い後に重量を最小化するのは、組合せ最適化で計算します。


Pythonで計算

まずは、ベンフォードの法則で100円~9999円の重み(wgt)を作ってみましょう。


python

import numpy as np, matplotlib.pyplot as plt

from math import log
# ベンフォードの法則
wgt = np.array([log((i+1)/i, 10000) for i in range(100, 10000)])
wgt /= sum(wgt)
plt.plot(wgt)
plt.xlabel('Price')

image

重量を最小化した後の個数を返すchangeを定義します。


python

from pulp import *

money_val = (1, 5, 10, 50, 100, 500, 1000, 5000)
money_wgt = (1.0, 3.7, 4.5, 4.0, 4.8, 7.0, 1.0, 1.0)
def change(price):
m = LpProblem() # 数理モデル
x = [LpVariable('x%d'%i, lowBound=0, cat=LpInteger)
for i in range(len(money_val))] # 支払い後の個数
m += lpDot(money_wgt, x) #目的関数(支払い後の重量)
m += lpDot(money_val, x) == price # 支払い後の金額
m.solve()
return [int(value(i)) for i in x]

シミュレーションしてみましょう。試しに1000円札の分布をみてみます。


python

price = 0 # 現在の所持金

warm, nrun = 100, 10000
res = []
for i, p in enumerate(rand_from_prob(wgt, warm+nrun)):
price -= p
if price < 0:
price += 10000
if price:
res.append(change(price))
a = np.array(res[-nrun:])
plt.hist(a[:,6], bins=5, range=(0, 5)) # 1000円札の分布

image

等確率ですね。他の硬貨や5000円札も同様でした。


結論

等確率になるようです。総額の分布も0円から9999円の等確率になります。


python

import pandas as pd

from itertools import product
r2, r5 = range(2), range(5)
ptn = [np.dot(money_val, n) for nn in
product(r5, r2, r5, r2, r5, r2, r5, r2)]
plt.hist(ptn)
print(pd.DataFrame(ptn).describe())
>>>
0
count 10000.00000
mean 4999.50000
std 2886.89568
min 0.00000
25% 2499.75000
50% 4999.50000
75% 7499.25000
max 9999.00000

image

以上