問題
クラスの委員が5人だけ決まっていません（体育が2人、図書が1人、保健が2人）。
時間もないので抽選することになりました。
なるべく生徒の希望に応じた確率で決めてください。
希望
- 生徒1：体育か図書を希望。第1希望が体育
- 生徒2：体育か保健を希望。第1希望が体育
- 生徒3：保健か体育を希望。第1希望が保健
- 生徒4：図書か保健を希望。第1希望が図書
- 生徒5：保健か体育を希望。第1希望が保健
ここでは、第1希望の確率を2/3、第2希望の方を1/3としましょう。
また、前提として、委員ごとの確率の和は委員の人数に一致することとします。
考え方
この割当問題を数理最適化でモデル化します。
まず、確率から何らかの重みを算出します。
目的関数は、その重みを上限とする一様乱数を係数とします。
目的関数を最大化することで、確率で割り当たるようにしてみましょう。
重みと確率
ここから少し寄り道をします。
2つの重みを2と3としましょう。
この重みから2つの一様乱数（範囲は0〜2と0〜3）aとbを生成します。
aとbを2次元で表すと、取りうる範囲は2 x 3で面積は6です。
aが大きくなる面積は2で、bが大きくなる面積は4です。
この2と4は、次のように関数func2()で計算できます。
In
import numpy as np
def func2(a, b):
p = a * a / 2
yield p
q = p + a * (b - a)
yield q
weight = np.array([2, 3])
print(*func2(*weight))
Out
2.0 4.0
実際に60万組の乱数を発生して確認してみましょう。
In
n = 600_000
rng = np.random.default_rng(0)
rnd = rng.random((n, weight.size))
print(np.unique((rnd * weight).argmax(axis=1), return_counts=True))
Out
(array([0, 1]), array([199555, 400445]))
おおよそ2対4になっています。
関数func2()の逆関数r_func2()は、式変形することで次のように作成できます。
In
def r_func2(p, q):
a = (p * 2)**(1 / 2)
yield a
b = (p + q) / a
yield b
rate = np.array([2, 4])
print(*(r_func2(*rate)))
Out
2.0 3.0
weightと同じ値が逆算されました。
重みの個数が3つのときは、次のようになります。これは、3次元の体積から計算できます。
In
def func3(a, b, c):
p = 0 + 1 * (a**3 - 0**3) / 3
yield p
q = p + a * (b**2 - a**2) / 2
yield q
r = q + a * b * (c**1 - b**1) / 1
yield r
def r_func3(p, q, r):
a = (3 * (p - 0) + 0**3)**(1 / 3)
yield a
b = (2 * (q - p) / a + a**2)**(1 / 2)
yield b
c = (1 * (r - q) / a / b + b**1)**(1 / 1)
yield c
rate = np.array([1, 2, 3])
print(*r_func3(*rate))
Out
1.4422495703074083 1.8619361889584654 2.2343234267501586
逆関数が正しいことを確認しましょう。
In
weight = np.array([*r_func3(*rate)])
print(np.array([*func3(*weight)]))
Out
[1. 2. 3.]
実際に60万組の乱数を発生して確認してみましょう。指定した比率になっています。
In
n = 600_000
rnd = rng.random((n, weight.size))
print(np.unique((rnd * weight).argmax(axis=1), return_counts=True))
Out
(array([0, 1, 2]), array([100019, 199960, 300021]))
重みの個数が4つのときは、次のようになります。計算式は、なんとなく規則性がわかるかと思います。
In
def func4(a, b, c, d):
p = 0 + 1 * (a**4 - 0**4) / 4
yield p
q = p + a * (b**3 - a**3) / 3
yield q
r = q + a * b * (c**2 - b**2) / 2
yield r
s = r + a * b * c * (d**1 - c**1) / 1
yield s
def r_func4(p, q, r, s):
a = (4 * (p - 0) / 1 + 0**4)**(1 / 4)
yield a
b = (3 * (q - p) / a + a**3)**(1 / 3)
yield b
c = (2 * (r - q) / a / b + b**2)**(1 / 2)
yield c
d = (1 * (s - r) / a / b / c + c**1)**(1 / 1)
yield d
rate = np.array([1, 2, 3, 4])
weight = np.array([*r_func4(*rate)])
n = 1_000_000
rnd = rng.random((n, weight.size))
print(np.unique((rnd * weight).argmax(axis=1), return_counts=True))
Out
(array([0, 1, 2, 3]), array([100125, 200424, 299249, 400202]))
解答
各生徒の希望確率を次のようにします。
wish_probs = [
{"体育": 2 / 3, "図書": 1 / 3},
{"体育": 2 / 3, "保健": 1 / 3},
{"保健": 2 / 3, "体育": 1 / 3},
{"図書": 2 / 3, "保健": 1 / 3},
{"保健": 2 / 3, "体育": 1 / 3},
]
この確率の割合に対する重みは、次のように2.0と1.5です。
In
print(*r_func2(2, 1))
Out
2.0 1.5
ここからは、実際に問題を解いていきます。
最初に、モデル作成で使う関数などを定義します。
from dataclasses import make_dataclass
import polars as pl
import numpy as np
from mip import Model, maximize, xsum
def xdot(a, b):
"""内積"""
a = a.to_numpy()
b = b.to_numpy()
return xsum(a * b)
def int_value(v):
"""整数解"""
return v.to_numpy().astype(float).round().astype(int)
rng = np.random.default_rng(0)
次に、数理モデルを作成します。
m = Model()
df = pl.from_dicts(
{"生徒": i, "委員": committee, "確率": prob}
for i, wish_prob in enumerate(wish_probs)
for committee, prob in wish_prob.items()
)
df = df.with_columns(
係数=pl.when(pl.col("確率") > 0.5).then(2.0).otherwise(1.5),
Var=m.add_var_tensor((len(df),), "Var", var_type="B"),
Val=0,
)
col = make_dataclass("Col", df.columns)(*map(pl.col, df.columns))
for _, group in df.group_by(col.委員):
# 委員ごとの人数
m += xsum(group["Var"]) == group["確率"].sum()
for _, group in df.group_by(col.生徒):
# 生徒ごとに1つ選択
m += xsum(group["Var"]) == 1
m.verbose = 0
最後に、目的関数を変えながら1000回解いて、生徒ごとの割合を見てみましょう。
In
result = []
for _ in range(1000):
coef = df["係数"] * rng.random(len(df))
m.objective = maximize(xdot(coef, df["Var"]))
m.optimize()
df = df.with_columns(Val=int_value(df["Var"]))
result.append(df.filter(col.Val == 1).sort("生徒")["委員"].to_list())
for lst in zip(*result):
print(np.unique(lst, return_counts=True))
Out
(array(['体育', '図書'], dtype='<U2'), array([699, 301]))
(array(['体育', '保健'], dtype='<U2'), array([589, 411]))
(array(['体育', '保健'], dtype='<U2'), array([336, 664]))
(array(['保健', '図書'], dtype='<U2'), array([301, 699]))
(array(['体育', '保健'], dtype='<U2'), array([376, 624]))
生徒同士は独立ではないので、必ずしも指定した割合になるとは限りませんが、近い値にはなっているようです。
まとめ
確率で割り当てる問題を数理最適化で解いてみました。
各委員の人数は決まっているため、生徒ごとに独立に決めることはできません。
数理モデルを使うことで、なるべく確率ごとに選んだり、委員の人数を指定できました。
以上