If you always pay to minimize the weight of your wallet, what is the distribution of your wallet?
Let's check it by simulation.
――10,000 yen is infinite. ――Buy the product 10,000 times and see the distribution. (Ignore the first 100 times) --The price of the product is 100 yen to 9999 yen, [Benford's Law](https://ja.wikipedia.org/wiki/%E3%83%99%E3%83%B3%E3%83%95 % E3% 82% A9% E3% 83% BC% E3% 83% 89% E3% 81% AE% E6% B3% 95% E5% 89% 87). --Use Walker's Alias Method to generate random numbers (rand_from_prob). --Minimize weight after payment is calculated by Combinatorial Optimization (http://qiita.com/Tsutomu-KKE@github/items/bfbf4c185ed7004b5721).
First, let's make a weight (wgt) of 100 yen to 9999 yen according to Benford's law.
python
import numpy as np, matplotlib.pyplot as plt
from math import log
#Benford's law
wgt = np.array([log((i+1)/i, 10000) for i in range(100, 10000)])
wgt /= sum(wgt)
plt.plot(wgt)
plt.xlabel('Price')
Defines a change that returns the number after minimizing the weight.
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() #Mathematical model
x = [LpVariable('x%d'%i, lowBound=0, cat=LpInteger)
for i in range(len(money_val))] #Quantity after payment
m += lpDot(money_wgt, x) #Objective function(Weight after payment)
m += lpDot(money_val, x) == price #Amount after payment
m.solve()
return [int(value(i)) for i in x]
Let's simulate. Let's try the distribution of 1000 yen bills.
python
price = 0 #Current money
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)) #Distribution of 1000 yen bills
It's an equal probability. The same was true for other coins and 5000 yen bills.
It seems to be equal probability. The distribution of the total amount also has an equal probability of 0 yen to 9999 yen.
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
that's all
Recommended Posts