プリパラで手持ちのプリチケからレア度ごと（N、R、SR、PR）の出現確率を推定する

このエントリは2014 1stライブから2014-2015 3rdライブまでのプリパラのシステムを前提として書いています。 というか2014-2015 3rdが稼働中に公開するつもりだったのが間に合いませんでした。 2015/4/2に稼働した2015 1stライブからはドリチケとPRよりも珍しいMR（ミステリーレア）が導入され少しシステムが変わっていますが、これらには対応しておりません。そういうものだと思って下さい。

このエントリでは理系学部レベルの数学を聞いたことがあることを前提としており、そういった箇所に関しては説明を省いている場合があります。 一方、学部では習わなそうな箇所やマニアックな箇所については説明を入れているつもりです。

また、

重要箇所

重要箇所についてはこのような表記を使っています。 markdownの引用表記です。 この表記を追うだけでなんとなくわかるように書いたつもりです。

はじめに

みなさんプリパラしていますか？ 私は沢山プリパラしているし、プリパラをしているとどんどん気持ちよくなってきます。

前回はコーデのレア度の出現確率を設定し、 シミュレーションを行うことでプリチケをコンプリートするために必要なプレイ回数の平均値を求めました。

mozo.hatenablog.com

しかし前回はほぼカンでレア度の出現確率を設定していました。

皆さんはPRやSRがどれくらいの確率で出現するのか知りたくありませんか？

今回はもう少し「マトモ」な方法でレア度ごとの出現確率を推定してみようと思います。

まず言葉の定義をしましょう。

• 「コーデ」とはお洋服のことを指します。プリパラでは1プレイで2つのコーデがランダムにプリズムストーンというお洋服屋さんに出現します。そして好きなほうを選んでプリチケにして入手できるというシステムになっています。

• 「プリチケ」とは入手したコーデであり、物理的なカードを指します。

• 「レア度」とは文字通りの意味です。各コーデに１つずつ割り振られています。現在はN（ノーマル）、R（レア）、SR（スーパーレア）、PR（プリパラレア）の4種類があり、右にいくほどレア度が高く出現しにくくなっていきます。しかし、各レア度の出現確率については明言されていません。

プリパラ問題？

手持ちのプリチケの情報を使って、レア度ごとの出現確率を推定せよ

今回の目標は各レア度の出現確率 \( P(Z) \)を推定することです。例えば、Nが50％、Rが40％、SRが7％、PRが3％といったような感じです。 実際にはレア度ごとにコーデは複数ありますし、同じレア度であってもコーデごとに出現確率が違うかもしれません。 本当は各コーデの出現確率を求めたいところですが、データ数が少ないこともあり、レア度に限定して推定します。これは暗に同じレア度のコーデは等しい確率で出現するという仮定を置いていることになります（たぶん）。

コーデが出現し、プリチケになるまでの過程を以下のようにモデル化します。 今後、「コーデ」を「レア度」と同じ意味で使う場合があります。 大文字が変数、小文字が具体的な値を表しています。

図１：プリパラ問題のモデル

最初に2つのレア度\( z_0 \)、\( z_1 \)の組み合わせが\( P(Z_0,Z_1) \)の確率でプリズムストーンに出現します。 ここで、2つコーデの出現確率は同一の確率変数\( Z \)に独立に従うという仮定を置いています。 即ち、\( P(Z_0|Z_1)=Z_0, P(Z_1|Z_0)=Z_1 \)であり、\( P(Z_0,Z_1)=P(Z_0)P(Z_1) \)となります。 PRが同時に2つでることだけは例外的に絶対にありえないとか、そういうことは考えないことにします。

この仮定をおかないとすると、プリズムストーンに出現するのは4種類のレア度の2つの重複在り組み合わせですから、4*4-1=15個のパラメータ（-1は確率の和は1であるという制約のため）を求める必要がでてきます。 しかし、この仮定をおくことによって、4-1=3個のパラメータで、プリズムストーンに出現するレア度の組み合わせを表現できます。 この仮定によって問題を解くのが容易になります。

今回の目的は\( P(Z) \)を推定することです。 あるプリチケを作った時、一緒に出たコーデのレア度は何だったのか私は記録していません。 つまり、\( P(Z) \)を頻度で直接推定することはできません。 このように\( Z \)は直接わからない、観測できない確率変数なので、一般的な用語を使うと隠れ変数であると言えます。 後に出てくる\( Y \)の影に潜んでいるので、潜在変数とも呼ばれます。

次に、プリズムストーンに出現した2つのレア度からプリチケにするレア度を条件付き確率\( P(Y|Z_0,Z_1) \)によって選択します。 ここでは、選択基準はレア度のみだという仮定を置いています。 実際には持っているかどうかやデザインが好みかどうかといった要素にも影響されるかと思いますが、それらは全て無視します。それらがこの確率に含まれていると考えてもよいでしょう。 また、プリパラには出現した２つのコーデをプリチケにする以外の選択肢（着ているコーデのスターランクアップ等）もありますが、今回は考慮しません。 \( P(Y|Z_0,Z_1) \)は、自分のこれまでの経験からカンで設定します。結局あやふやで適当じゃないか、と言われれば全くその通りです。 しかし、\( P(Z) \)を推定するために、\( P(Y|Z_0,Z_1) \)という分布を用いることは、 前回のように「Nと一緒にR出現した気がする回数」なんかを使うよりは理論的なバックグラウンドがあるような気にならないでしょうか？

最終的に、手元のマイチケがあるレア度である確率を\( P(Y) \)とおきます。これは手元のマイチケを集計した、レア度の頻度として定義できます。

今回の目標は\( P(Y|Z_0,Z_1) \)と\( P(Y) \)がわかっていて、これらから\( P(Z) \)を推定することです。

プリパラ問題（PP: Pripara Problem）

離散確率変数\( Z = \{ N,R,SR,PR \} \)と \( Y = \{ N,R,SR,PR \} \)が存在し、\( Z_0 \)と\( Z_1 \)がそれぞれ独立に\( Z \)に従うとする。 \( P(Y) \)と\( P(Y|Z_0,Z_1) \)が既知である時、\( P(Z) \)を推定せよ

余談になりますが、上の図はグラフィカルモデルの記法を使うと以下のように描画できます。 白抜きのノードは観測できない変数、影付きのノードは観測されている変数を表します。

図２：プリパラ問題をグラフィカルモデルの表記を用いて示した図

一般化プリパラ問題

一度プリパラの状況を離れて、もう少し一般的なシチュエーションを考えてみましょう。 一般的な形でこの問題を書くと以下のようになります。

一般化プリパラ問題（GPP: Generalized Pripara Problem）

離散確率変数\( Z = \{ z_0, \cdots , z_{n-1} \} \)と \( Y = \{ y_0, \cdots ,y_{m - 1} \} \)が存在し、\( P(Y|Z) \)と\( P(Y) \)が既知であるとき、\( P(Z) \)を推定せよ

この一般的な形式では先ほどまでの\( Z_0 \)と\( Z_1 \)のペアをまとめて\( Z \)と呼んでいるので注意して下さい。 ググってもこの問題が見つからなかったので名前をつけてみました。

条件付き確率といえば、そう、ベイズですね。 ベイズ統計学ではベイズの定理を用いて、事後分布\( P(Z|Y) \)を事前分布\( P(Z) \)と尤度\( P(Y|Z) \)の積で推定することがよく行われます。

\( \begin{align*} P(Z|Y) = P(Z)\frac{P(Y|Z)}{P(Y)} = P(Z)\frac{P(Y|Z)}{\sum_{Z}P(Y|Z)} \propto P(Z)P(Y|Z) \end{align*} \)

しかし、今回は\( P(Z) \)が知りたいのでこの有名な式を使うことはできなさそうです。

仕方がないので、わかっている確率の関係をまとめてみましょう。

\( P(Y)=\sum_{Z}P(Y,Z)=\sum_{Z}P(Z)P(Y|Z) \)

\( P(Y) \)を結合確率\( P(Y,Z) \)の\( Z \)についての和で表し、さらに結合確率を確率と条件付き確率の積に分解しました。 いたって普通の式変形をしているだけです。ここまでくるとなんだか解けそうな気がしてきます。

次に\( Z \)についての和を展開してみましょう。

\( P(Y)=P(Z=z_0)P(Y|Z=z_0)+ \cdots + P(Z=z_{n-1})P(Y|Z=z_{n-1}) \)

\( Z \)に具体的な値を入れて和を展開しました。\( P(Z=z_i) \)は\( P(Z) \)において\( Z=z_i \)である時の具体的な確率値を表します。 さらに今度は\( Y \)に具体的な値を入れてみましょう。

\( {\scriptstyle \left \{\begin{align} P(Y=y_0) &= P(Z=z_0)P(Y=y_0|Z=z_0)+ \cdots + P(Z=z_{n-1})P(Y=y_0|Z=z_{n-1}) \\ \vdots \\ P(Y=y_{m - 1}) &= P(Z=z_0)P(Y=y_{m - 1}|Z=z_0)+ \cdots + P(Z=z_{n-1})P(Y=y_{m - 1}|Z=z_{n-1}) \end{align}\right. } \)

というわけで連立方程式になりました。 \(P(Z) = \{ P(Z=z_0), \cdots , P(Z=z_{n-1}) \}\)は変数ベクトルとみなせることに注意して下さい。 一般化プリパラ問題の解はこの連立方程式を解けば求められます。

とその前に、実はこの式は少し冗長です。確率は和が1になると決まっている\( \sum_{Z} P(Z)=1 \) ので、 \( P(Z=z_{j}) = 1 - \sum_{i \neq j}P(Z=z_{i}) \) という関係があります。 つまり\( Z \)が取りうる具体的な値は\( n \)個ありますが、最後の1個についての確率は他の値についての確率が決まれば一意に決まるので、パラメータ数は\( n - 1\)個で済みます。 \( Y \)についても同様で、パラメータ数は\( m - 1\)個で済みます。

そこで、先ほどの連立方程式に \( P(Z=z_{n - 1}) = 1 - \sum_{0 \leq i \leq n-2}P(Z=z_{i}) \) を代入して式を変形し、さらに一番下の冗長な式を取り除くことで式と変数の数を減らせます。 変形した式を以下に示します。パッと見ではごちゃごちゃしてまいました。

一般化プリパラ問題（GPP: Generalized Pripara Problem）の解法

一般化プリパラ問題の解\(P(Z) = \{ P(Z=z_0), \cdots , P(Z=z_{n-1}) \}\)は以下の\( m - 1\)個の式からなる\( n- 1\)元連立方程式を解くことによって求められる。 \( {\scriptstyle \left \{\begin{align} P(Y=y_0) - P(Y=y_0|Z=z_{n-1}) &= P(Z=z_0)(P(Y=y_0|Z=z_0) - P(Y=y_0|Z=z_{n-1}))+ \\ & \cdots + P(Z=z_{n - 2})(P(Y=y_0|Z=z_{ n - 2}) - P(Y=y_0|Z=z_{n-1})) \\ \vdots \\ P(Y=y_{m - 2}) - P(Y=y_{m - 2}|Z=z_{n-1}) &= P(Z=z_0)(P(Y=y_{m - 2}|Z=z_0) - P(Y=y_{m - 1}|Z=z_{n-1}))+ \\ & \cdots + P(Z=z_{n - 2})(P(Y=y_{m - 2}|Z=z_{ n - 2}) - P(Y=y_{m - 2}|Z=z_{n-1})) \end{align}\right. } \)

\( P(Z=z_{n - 1}) = 1 - \sum_{0 \leq i \leq n-2}P(Z=z_{i}) \)

ただし、\( 0 \leq P(Z=z_{i}) \leq 1 \)

ここまでくるとありふれた線形代数の話になります。 一般化プリパラ問題が解けるかどうかは状況によります。 例えば、\( n=m \)であり、係数行列が正方行列かつ逆行列を求めることができれば、解は唯一に定まります（ただし、全ての値が正かつ和が1になるという確率としての制約を満たさなければ解なし）。 しかし、それ以外の場合は、解が不定（無数の解が存在する）か、解を持ちません。

プリパラ問題

一般化プリパラ問題の知見を用いてプリパラ問題を数式に落としていきましょう。

式がかなりごちゃごちゃしますから、今後、以下の簡易的な記法でこれまでの表記を置き換えることがあります。ご了承下さい。

\( \left \{\begin{align} Z &= \{ z_0 \overset{\mathrm{def}}{=} N, z_1 \overset{\mathrm{def}}{=} R , z_2 \overset{\mathrm{def}}{=} SR, z_3 \overset{\mathrm{def}}{=} PR \} \\ Y &= \{ y_0 \overset{\mathrm{def}}{=} N, y_1 \overset{\mathrm{def}}{=} R , y_2 \overset{\mathrm{def}}{=} SR, y_3 \overset{\mathrm{def}}{=} PR \} \\ x_i & \overset{\mathrm{def}}{=} P(Z=z_i) \\ p(y_i) & \overset{\mathrm{def}}{=} P(Y=y_i) \\ p_(y_i|z_j,z_k) & \overset{\mathrm{def}}{=} P(Y=y_i|Z_0=z_j,Z_1=z_k) \\ p_{\scriptscriptstyle PR}(y_i)\ & \overset{\mathrm{def}}{=} P(Y=y_i|Z_0=z_i,Z_1=z_3) \end{align}\right. \)

まず、実際のプリパラ問題は、一般化プリパラ問題の\( Z\)を\( (Z_0,Z_1) \)に読み替えた場合に相当します。 一般化プリパラ問題の\( P(Z) \)は、プリパラ問題では\( P(Z_0,Z_1)=P(Z_0)P(Z_1) \)となります。

これは困ったことになりました。変数の積、非線形な要素が現れています。 こうなるともはや”線形”代数では扱うことができません。 つまり、逆行列を求めて終わりではありません。 不安ですが、ちょっとこのまま先に進んでみます。

次に、\( P(Y|Z_0,Z_1) \)について考えてみましょう。 図１を見て下さい。 プリズムストーンにはSRとPRが出現しています。 この場合、プリチケになりうるのはSRかPRだけで、NとRは絶対にプリチケにはなりません。 \( P(Y=N|Z_0=SR,Z_1=PR) = P(Y=R|Z_0=SR,Z_1=PR) = 0\)です。 つまり、プリズムストーンに出現していないレア度がプリチケになることは絶対にありません。

また、プリズムストーンにNが2つ出現した場合を考えて下さい。 この場合、プリチケになりうるのはNだけで、それ以外は絶対にプリチケになりません。 \( P(Y=N|Z_0=N,Z_1=N) = 1, P(Y \neq N |Z_0=N,Z_1=N) = 0\)です。 つまり、プリズムストーンにあるレア度が２つ出現した場合、そのレア度をプリチケにするしかありません。

\( P(Y|Z_0,Z_1) \)には最終的に私の主観を与えて計算しますが、上の２つの条件に関しては、誰の主観で あれ変わることがないでしょう。そういうシステムですから～～👓

以上をまとめてみましょう。

プリパラ問題において、以下が成り立つ

\( P(Z_0,Z_1)=P(Z_0)P(Z_1) \)

\( P(Y=y_{k \neq i, j}|Z_0=y_i,Z_1=y_j) = 0\)

\( P(Y=y_i|Z_0=y_i,Z_1=y_i) = 1, P(Y = y_{j \neq i } |Z_0=y_i,Z_1=y_i) = 0 \)

これを一般化プリパラ問題の式に代入すると、0によって項の数がかなり減ります。 整理すると、以下のような非線形連立方程式が得られます。

プリパラ問題（PP: Pripara Problem）の解法

プリパラ問題の解\(P(Z) = \{ x_0, x_1, x_2, x_3 \}\)は以下の3元連立非線形方程式を解くことで求められる。 \( {\scriptstyle \left\{ \begin{array}{l} f_0 = (1 - 2p_{\scriptscriptstyle PR}(y_0))x_0x_0 + 2(p(y_0 |z_0, z_1) - p_{\scriptscriptstyle PR}(y_0))x_0x_1 + 2(p(y_0 | z_0, z_2) - p_{\scriptscriptstyle PR}(y_0))x_0x_2 + 2p_{\scriptscriptstyle PR}(y_0)x_0 - p(y_0) = 0 \\ f_1 = 2(p(y_1|z_0, z_1) - p_{\scriptscriptstyle PR}(y_1))x_0x_1 + (1 - 2p_{\scriptscriptstyle PR}(y_1))x_1x_1 + 2(p(y_1 | z_1, z_2) - p_{\scriptscriptstyle PR}(y_1))x_1x_2 + 2p_{\scriptscriptstyle PR}(y_1)x_1 - p(y_1) = 0 \\ f_2 = 2(p(y_2|z_0, z_2) - p_{\scriptscriptstyle PR}(y_2))x_0x_2 + 2(p(y_2 |z_1, z_2) - p_{\scriptscriptstyle PR}(y_2))x_1x_2 + (1 - 2p_{\scriptscriptstyle PR}(y_2))x_2x_2 + 2p_{\scriptscriptstyle PR}(y_2)x_2 - p(y_2) = 0 \end{array} \right. } \) \( x_3 = 1 - (x_0 + x_1 + x_2) \)

ただし\( 0 \leq x_i \leq 1 \)かつ\( x_0 > x_1 > x_2 > x_3 \)を満たす

この先のことを考えて左辺=0という形にし、さらに式に名前をつけています。 最後の制約は、確率値は0以上1以下であり、さらにレア度が高いほど確率が低いという必須条件を表しています。 変形の際は\(p(y_i | z_j, z_k) = p(y_i|z_k, z_j)\)を用いました。2という係数がでてくるのは このためです。

これでやっとプリパラ問題を数式に落とすことができました。 しかしこれは非線形連立方程式です。理系の学部生なら誰もが習うような線形代数の初歩的な手法を使って解くことはできません。 そこで、次のチャプターでは非線形連立方程式の解法について説明します。

非線形連立方程式

非線形連立方程式は複雑すぎて一般に解析的に解くことはできません。 イテレーションを伴う数値計算によって近似解を求めるしかありません。 このチャプターでは二種類の解き方を紹介します。

逆運動学（IK: Inverse Kinematics）とレーベンバーグ・マーカート（LM：Levenberg–Marquardt）法

工学的な応用を考えた時に、非線形連立方程式がよくでてくるのが逆運動学（IK: Inverse Kinematics）という分野です。 ロボティクスにおけるロボットアームや、ゲームやアニメに出てくるキャラクターの人体3DCGモデルといった関節がある物体を考えて下さい。 ある関節iを\( \theta_i \)°曲げることに伴い、その関節の先に繋がっている部分が動き、ある座標xに移動します。θからxを計算することは簡単にできます。

三次元空間を対象とし、関節が複数ある場合を考えましょう。 \( { \bf \Theta } \)を関節パラメータベクトルとすると関節に繋がっている部位の座標\( \{x,y,z\} \)は関節の動きを表す非線形関数\( f_x, f_y, f_z \)によって

\( \left\{ \begin{array}{l} x & = f_x({\bf\Theta}) \\ y & = f_y({\bf\Theta}) \\ z & = f_z({\bf\Theta}) \end{array} \right. \)

と計算できます。 これを順運動学（FK: Forward Kinematics）と呼びます。

逆に位置\( \{x,y,z\} \)からパラメータ\( { \bf \Theta } \)を求めることをIKと呼びます。 3DCGモデルを扱っている人にとっては馴染みがある用語かもしれません。 上の式をみれば、IKが非線形連立方程式を解くことと等価であることがわかるかと思います。

IKは簡単ではありません。何故ならば解がいくつもあるからです。 憎たらしい人物の顔面に拳でパンチをすることを想定して下さい。同じ位置にパンチをするのでも、斜め方向からのフックなのか正面からのストレートなのか、いくつもの腕の関節の角度が考えられますね。 さらに、足や腰も動かすとなるとより多くの関節角度の組み合わせが考えられます。 腕を前に伸ばした結果は一意にある位置へのパンチになりますが（FK）、ある位置へのパンチをするにはやり方がいくつもある（IK）のです。 制約は\( \{x,y,z\} \)の3つですが、一般にパラメータ（関節）のほうが数が多いため、いくつもの解ができてしまうのです。

IKを日本語で詳しく説明した資料はあまりないようです。 英語ですが、以下の資料は参考になりました。 http://www.math.ucsd.edu/~sbuss/ResearchWeb/ikmethods/iksurvey.pdf

IK=非線形連立方程式=プリパラ問題を解く上で性能が良いのはDamped least squares、別名を レーベンバーグ・マーカート（LM：Levenberg–Marquardt）法と呼ばれる手法のようです。

LM法を含む最適化については日本語のわかりやすい参考書があります。 www.amazon.co.jp

LM法は2乗誤差の和を最小にするパラメータ\( \hat{ \bf \Theta } \)を推定する方法で、非線形最小2乗法の1種です。 上の連立方程式と同じ記法を無理やり使うと、LM法が最小化する二乗誤差の和は以下のように書けます。

\(J = \frac{1}{2} \sum_{i=x,y,z} ( f_{i}( { \bf \Theta }) - i )^{2} \)

最小2乗法による直線のフィッティングはやったことがある人も多いかと思いますが、 直線ではなく任意の非線形関数でフィッティングする方法の1つとしてLM法と捉えることができます。

一方でアルゴリズムの動きを理解する上では、LM法は\( J \)を目的関数とした勾配ベースの局所探索アルゴリズムの1種と捉えたほうがよいでしょう。 初期値からスタートして、局所的な傾きをみて、目的関数が減る方向にどんどんパラメータをアップデートしていくという方法です。 LM法は、解から遠い場所では一次微分（ヤコビアン、ヤコビ行列）による粗い探索（最急降下法）を行い、 解の近傍ではヤコビアンで近似した二次微分（ヘシアン、ヘッセ行列）を用いてより細かく探索（ガウス・ニュートン法）します。

進化計算（EC：Evolutionary Computation）と差分進化（DE: Differential Evolution）

非線形連立方程式を解く上で実用上最も有用な手法はLM法であることがわかりました。 しかし、LM法はあくまで局所探索です。局所的な傾きをみているだけなので、局所解にハマったら抜け出すことはできません。 アイドルたるもの、局所解にハマってしまってもよいのでしょうか？

1期最終回で実はプリパラはいくつも存在することが明らかになりました。 あるプリパラでは一等賞（局所解）でも、他のプリパラにはもっと強いアイドルがいるかもしれません。 全てのプリパラで一等賞になる＝大域最適解を探すことを目指すべきではないでしょうか？

大域最適化の手法としてポピュラーなのが進化計算（EC：Evolutionary Computation）です。 遺伝的アルゴリズム（GA: Genetic Algorithm）の仲間たち、といえばわかる人にはわかって頂けるでしょうか。

局所的な情報である勾配だけをみていては大域最適解は絶対に見つかりません。 進化計算のベースとなるアルゴリズムは基本的に勾配を使いません。 探索空間にいくつも解をばらまき、良さそうな解を組み合わせて解をアップデートしていきます。 このアップデート法が生物進化の考えにインスパイアされている場合が多いので、進化計算と呼ばれます。 進化計算という名前はインパクトがある一方、誤解を生みやすいと私は思います。 あまりにも抽象化されているからです。 アルゴリズムが自動的にヒューリスティックを見つけるので、メタヒューリスティックなアルゴリズムであるとも言われます。

LM法などの勾配ベースの手法と比較すると、進化計算は良い解を見つけることができるのですがコストパフォーマンスはよくないです。 探索空間に多くの解をばらまくので、必然的に計算量がかかります。

実用上、勾配ベースではなく進化計算を使うべきシチュエーションとしては

• 勾配に意味がない=局所解が多すぎる場合
• 計算量度外視でいい解がほしい場合
• 解析的に勾配を求められず、勾配の数値計算にコストがかかりすぎる場合

などが考えられます。

今回はLM法だけでは不安なので、比較するために進化計算を使ってみます。 上記の二番目の理由になるでしょうか。

進化計算の手法はいくつもあります。それこそトンボとか狼とか、生き物の数だけあるといってもいいでしょう。 天下り的ですが最近scipyに実装されたようなので、今回は差分進化（DE: Differential Evolution）を使ってみます。 DEは実用的な性能を持った進化計算の1つです。

DEは進化計算といいつつ進化や生き物っぽくありません。 探索空間にばらまいた解同士の差分ベクトルを使ってランダムに解をアップデートし、その中から良い解を選んでさらにアップデートしてきます。 差分ベクトルを使うことで自動的に最適な探索方向やステップ幅が決定されるそうです。

DEに関する日本語の資料はネット上に沢山あるのですが、どれもイマイチわかりやすくありません。 これでもわかりやすいとは言えませんが、図を使った説明資料として以下のページを貼っておきます。産業界でも使われているようですね。

マルチドメインソリューション：ソリューション例：Differential Evolution による最適化事例：サイバネット

今回はDEでもLM法と同様に\( J \)を目的関数として設定し、最小化します。 \( J \)を最小化するという目的は同じですが、やり方が異なるわけです。

実数値ベクトルを引数とし実数の値を返すような関数を最適化するため実用的な進化計算手法として、DEの他には粒子群最適化（PSO: Particle Swarm Optimization）やCMA-ESがあります。 進化計算には組み合わせ最適化向けの手法もあり（というか元々はこっちが先なのですが）GAや蟻コロニー最適化（ACO: Ant Colony Optimization）等があります。

非線形連立方程式の解法

非線形連立方程式の解法をまとめると以下のようになります。

非線形連立方程式の解法

\( F_i({\bf\Theta}) \)を任意の非線形関数、\( C_i \)を定数とする。

\( \left\{ \begin{array}{l} F_0({\bf\Theta}) = C_0\\ \vdots \\ F_{ n - 1}({\bf\Theta}) = C_{n - 1} \end{array} \right. \)

なる非線形連立方程式があるとき、

\( \left\{ \begin{array}{l} f_0({\bf\Theta}) = F_0({\bf\Theta}) - C_0 = 0\\ \vdots \\ f_{n - 1}({\bf\Theta}) =F_{ n - 1}({\bf\Theta}) - C_{n - 1} = 0 \end{array} \right. \)

とおく。

この非線形連立方程式の近似解は、二乗誤差の和

\(J = \frac{1}{2} \sum_{i} ( f_{i})^{2} \)

を最小にする\( \bf\Theta \)である。

\( J \)の最小化にはLM法やDEを用いることができる。

真の解では\(J = 0 \)のはずですから、\( J \)が0に近いほど真の解に近い、良い近似解が得られているということになります。

実験と考察

実際の数値を使い、プログラムを書いて、プリパラ問題の解を求めてみましょう。

私の環境におけるプリパラ問題

さて、これまでは変数を使って抽象的な話をしてきましたが、 ここからは実際の数値を具体的に考えていきましょう。

\( P(Y) \)を手持ちのプリチケから具体的に推定しましょう。 前回のブログによると、少し前の私は Nを24枚、Rを19枚、SRを10枚、PRを7枚、計60枚のプリチケを持っていたようです。 現在のプリチケの枚数は想定で数倍になっているのですが、数えるのがめんどくさいのでこの値を使いましょう。 単純にレア度の頻度を算出します。 \( P(Y=N)=\frac{24}{60}=0.4 \)、 \( P(Y=R)=\frac{19}{60}=0.3167 \)、 \( P(Y=SR)=\frac{10}{60}=0.1667 \)、 \( P(Y=PR)=\frac{7}{60}=0.1167 \)となります。

次に、\( P(Y|Z_0,Z_1) \)を具体的に以下の表のように与えます（表が綺麗につくれませんでした。markdownで綺麗な表を作るにはどうすればよいのでしょうか？）。このように離散の条件付き確率は表で表す（実装する）ことができ、条件付き確率表（CPT：Conditional Probability Table）と呼ばれます。

表１：\( P(Y|Z_0,Z_1) \)の条件付き確率表

	(N, N)	(N, R), (R, N)	(N, SR), (SR, N)	(N, PR), (PR, N)	(R, R)	(R, SR), (SR, R)	(R, PR), (PR, R)	(SR,SR)	(SR, PR), (PR, SR)	(PR, PR)
\(Y=N \)	1	0.4	0.1	0.01	0	0	0	0	0	0
\(Y=R \)	0	0.6	0	0	1	0.2	0.05	0	0	0
\(Y=SR \)	0	0	0.9	0	0	0.8	0	1	0.2	0
\(Y=PR \)	0	0	0	0.99	0	0	0.95	0	0.8	1

この表は列（縦軸）が\( (Z_0,Z_1) \)の具体的な値、行（横軸）が\( Y \)の具体的な値、あるセルが\( P(Y|Z_0,Z_1) \)の値を示しています。 左から二列目を見て下さい。 この列は\( P(Y|Z_0=N,Z_1=R) \)もしくは\( P(Y|Z_0=R,Z_1=N) \)である場合、つまりプリズムストーンにNとRが出現した時にプリチケにする確率を表しています。 この列の上から1列目をみるとNをプリチケにする確率\( P(Y=N|Z_0=N,Z_1=R) \)が0.4、2列目を見ると\( P(Y=R|Z_0=N,Z_1=R) \)が0.6であることがわかります。 レア度が高いRのほうをプリチケにする確率を高く設定しています。 3列目と4列目は\( P(Y=SR|Z_0=N,Z_1=R) =P(Y=PR|Z_0=N,Z_1=R)=0 \)ですが、NとRが出現したのにSRやPRがプリチケになることは絶対にないので0としています（これは先ほど述べました）。

このようにある列に注目すると、2つが非0であり、レア度が高いほうに大きな確率が割り振られています。 レア度の差が大きくなるほど、割り振られる確率値の差も大きくなります。例えば、NとPRが出現した時はほぼ確実にPRを取るはずなのでPRを取る確率を99%、Nを取る確率を1%に設定しています。 例外はNとNが出現した時のように\( (Z_0,Z_1) \)がどちらも同じレア度をとる場合で、この場合は選択の余地がないので、あるレア度が1にその他が0になります（これも先ほど述べましたね）。 どの列も和を取ると1になります。

以上の数値を先ほどのプリパラ問題の式に代入すると以下のようになります。

私の環境におけるプリパラ問題（PP: Pripara Problem）の解法

\( x_0 \)をNの出現確率、\( x_1 \)をRの出現確率、 \( x_3 \)をSRの出現確率、 \( x_3 \)をPRの出現確率とする。 私の手持ちのプリチケと、私のレア度を選択する感覚を用いると、 \( \{ x_0, x_1, x_2 \}\)は以下の3元連立非線形方程式の解である。

\( \left\{ \begin{array}{l} f_0 = 0.98x_0x_0 + 0.78x_0x_1 + 0.18x_0x_2 +0.02x_0 - 0.4 = 0 \\ f_1 = 1.1x_0x_1 + 0.9x_1x_1 + 0.3x_1x_2 +0.1x_1- 0.3167 = 0 \\ f_2 = 1.4x_0x_2 + 1.2x_1x_2 + 0.6x_2x_2 +0.4x_2- 0.1667 = 0 \end{array} \right. \)

以上の連立方程式を解くことは

\( J = \frac{1}{2}(f^2_0 + f^2_1 + f^2_2) \)

を最小化する\( \{ x_0, x_1, x_2 \}\)を見つけることと等価であり、 このエントリでは最小化にLM法もしくはDEを用いる。

また、\( x_3 = 1 - (x_0 + x_1 + x_2) \)である。

ただし解として意味があるのは\( 0 \leq x_i \leq 1 \)かつ\( x_0 > x_1 > x_2 > x_3 \)を満たす解だけである

あとはこれを実際にプログラミングして解いてみましょう。

私の環境におけるプリパラ問題をscipyのLM法とDEを使って解く

さて、いよいよ実際にプリパラ問題を解いてみましょう。 pythonのscipyにLM法とDEが実装されているのでそれらを使ってみましょう。 ただしDEについてはscipy 0.15.0以降でないと使えないので注意して下さい。

コードは長いので最後に載せます。 重要なところだけをみていきましょう。

ポイントとしては、まずleastsq（LM法）とdifferential_evolution（DE）をimportします。

from scipy.optimize import leastsq
from scipy.optimize import differential_evolution

LM法を使っているのはこの箇所です。二乗誤差を計算する関数ではなく、残差\( f_i \)をリストで返す関数を渡す必要があることに気をつけて下さい。calcResidualは\( [f_0, f_1, f_2] \)を返します。

#残差計算関数、初期値、残渣計算に用いるデータを渡す
pbest = leastsq(calcResidual, p0, args=(P_N, P_R, P_SR), full_output=1)

DEを使っているのはこの箇所です。こちらにはLMが内部で使っていると思われる二乗誤差の和をとる関数、つまり\( J \)を渡してあげます。calcFitnessがそれです。また、探索範囲も渡す必要があり、今回は全ての変数について[0,1]を探索範囲としています。

#評価関数、探索範囲、評価関数に渡す追加の引数を渡す
deResult = differential_evolution(calcFitness, [(0,1)]*3, args=(P_N, P_R, P_SR))

このコードではLM法とDEを10回ずつ実行しますが、LM法でもDEでも毎回ほぼ同じ結果が返ってきます。 しかも二乗誤差の和\( J \)もほぼ0です。ほぼ真の解だといってよいでしょう。

各パラメータが[0,1]の範囲では、目的関数\( J \)は唯一の局所解=大域最適解を持っていたようです。 DEを使う必要はなかったですね。 ある実行例を下に示します。

P(N) = 0.50741041685

P(R) = 0.322879972067

P(SR)= 0.106710586887

P(PR)= 0.0629990241973

error: 3.08148791102e-33

これは

ただし\( 0 \leq x_i \leq 1 \)かつ\( x_0 > x_1 > x_2 > x_3 \)を満たす

という制約も満たしています!

これでプリパラ問題を解くことができました！

私の環境におけるプリパラ問題（PP: Pripara Problem）の解

私の手持ちのプリチケと、私のレア度を選択する感覚を用いると、プリズムストーンに

Nが出現する確率は51%

Rが出現する確率は32%

SRが出現する確率は11%

PRが出現する確率は6%

である。

少しだけ南委員長に近づくことができた気がします。 ありがとうプリパラ。はやくプリパラを電極で脳に繋ぎたい。

プリパラ問題を解くコード

python3.4, scipy0.15.0を使っています。

import random
from scipy.optimize import leastsq
from scipy.optimize import differential_evolution

#手持ちのプリチケ枚数
pRareNum =  7
sRareNum = 10
rRareNum = 19
nRareNum = 24

allNum = float(pRareNum + sRareNum + rRareNum + nRareNum)

#P(Y)を手持ちのプリチケから計算
P_N  = [nRareNum / allNum]
P_R  = [rRareNum / allNum]
P_SR = [sRareNum / allNum]
P_PR = [pRareNum / allNum]

rarities = ['N', 'R', 'SR', 'PR']

#P(Y|Z0,Z1)を保持するクラス
class ConditionalProb:
    def __init__(self):
        self.cond = {}
        for instance in rarities:
            self.cond[instance] = {}
            for r0 in rarities:
                self.cond[instance][r0] = {}
                for r1 in rarities:
                    self.cond[instance][r0][r1] = 0
    
    def initialize(self):
        self.setCond('N', ('N', 'N'), 1)
        self.setCond('N', ('N', 'R'), 0.4)
        self.setCond('N', ('N', 'SR'), 0.1)
        self.setCond('N', ('N', 'PR'), 0.01)
        
        self.setCond('R', ('R', 'N'), 0.6)
        self.setCond('R', ('R', 'R'), 1)
        self.setCond('R', ('R', 'SR'), 0.2)
        self.setCond('R', ('R', 'PR'), 0.05)

        self.setCond('SR', ('SR', 'N'), 0.9)
        self.setCond('SR', ('SR', 'R'), 0.8)
        self.setCond('SR', ('SR', 'SR'), 1)
        self.setCond('SR', ('SR', 'PR'), 0.2)

        self.setCond('PR', ('PR', 'N'), 0.99)
        self.setCond('PR', ('PR', 'R'), 0.95)
        self.setCond('PR', ('PR', 'SR'), 0.8)
        self.setCond('PR', ('PR', 'PR'), 1)
        
    def setCond(self, instance, conditionPair, p):
        self.cond[instance][conditionPair[0]][conditionPair[1]] = p
        self.cond[instance][conditionPair[1]][conditionPair[0]] = p
        
    def getCond(self, instance, conditionPair):
        return self.cond[instance][conditionPair[0]][conditionPair[1]]
    
conditionalProb = ConditionalProb()
conditionalProb.initialize()

def calcResidualN(p, n):
    c = conditionalProb

    pN = c.getCond('N', ('N', 'N'))
    pR = c.getCond('N', ('N', 'R'))
    pSR = c.getCond('N', ('N', 'SR'))
    pPR = c.getCond('N', ('N', 'PR'))
    coefNN = pN - 2 * pPR
    coefNR = 2 * (pR - pPR)
    coefNSR = 2 * (pSR - pPR)
    
    residual = coefNN*p[0]*p[0] + coefNR*p[0]*p[1] + coefNSR*p[0]*p[2] \
                + 2*pPR*p[0] - n[0]

    return residual


def calcResidualR(p, r):
    c = conditionalProb

    pN = c.getCond('R', ('R', 'N'))
    pR = c.getCond('R', ('R', 'R'))
    pSR = c.getCond('R', ('R', 'SR'))
    pPR = c.getCond('R', ('R', 'PR'))
    coefRN = 2 * (pN - pPR)
    coefRR = pR - 2 * pPR
    coefRSR = 2 * (pSR - pPR)
    
    residual = coefRN*p[0]*p[1] + coefRR*p[1]*p[1] + coefRSR*p[1]*p[2] \
                + 2*pPR*p[1] - r[0]

    return residual

def calcResidualSR(p, sr):
    c = conditionalProb

    pN = c.getCond('SR', ('SR', 'N'))
    pR = c.getCond('SR', ('SR', 'R'))
    pSR = c.getCond('SR', ('SR', 'SR'))
    pPR = c.getCond('SR', ('SR', 'PR'))
    coefSRN = 2 * (pN - pPR)
    coefSRR = 2 * (pR - pPR)
    coefSRSR = pSR - 2 * pPR
    
    residual = coefSRN*p[0]*p[2] + coefSRR*p[1]*p[2] + coefSRSR*p[2]*p[2] \
                + 2*pPR*p[2] - sr[0]

    return residual

#residual（残差）、すなわち0にすべき各式の誤差を計算する。
#二乗して和をとるのはLM法が内部でやってくれる
def calcResidual(params, n, r, sr):
    residualN = calcResidualN(params, n)
    residualR = calcResidualR(params, r)
    residualSR = calcResidualSR(params, sr)

    residual = [residualN, residualR, residualSR] 

    return residual

#LM法の初期値を生成する。確率の和が1以下であるという制約を観ますように生成する。
def createInitialParams():
    pN  = random.random()
    pR  = random.uniform(0, 1 - pN, )
    pSR = random.uniform(0, 1 - pN - pR)
    return [pN, pR, pSR]

#DE用の誤差関数。LMが内部で使う誤差関数と同じ
def calcFitness(params, n, r, sr):
    fitness = 0
    residuals = calcResidual(params, n, r, sr)
    for r in residuals:
        fitness += r**2
    return fitness


if __name__ == '__main__':
    
    for i in range(1, 11):
        print("LM method" ,"run:", i)
        p0 = createInitialParams()
        print("init params ", p0)
        #残差計算関数、初期値、残渣計算に用いるデータを渡す
        pbest = leastsq(calcResidual, p0, args=(P_N, P_R, P_SR), full_output=1)
        bestparams = pbest[0]
        n, r, sr, pr = bestparams[0], bestparams[1], bestparams[2],  1-sum(bestparams)
        print("P(N) =",   n)
        print("P(R) =",   r)
        print("P(SR)=", sr)
        print("P(PR)=", pr)
        print("error: ", calcFitness([n, r, sr], P_N, P_R, P_SR))
        print()
    
    print()
    for i in range(1, 11):
        print("DE", "run:", i)
        #評価関数、探索範囲、評価関数に渡す追加の引数を渡す
        deResult = differential_evolution(calcFitness, [(0,1)]*3, args=(P_N, P_R, P_SR))
        print(deResult["x"], deResult["fun"])
        bestparams = deResult["x"]
        n, r, sr, pr = bestparams[0], bestparams[1], bestparams[2],  1-sum(bestparams)
        print("P(N) =",   n)
        print("P(R) =",   r)
        print("P(SR)=", sr)
        print("P(PR)=", pr)
        print("error: ", calcFitness([n, r, sr], P_N, P_R, P_SR))
        print()