プリズム☆ジャンプ11
プリジャンでよいと思った本の感想をネタバレありで書きます。
どきどき！ま夜中のメイキングドラマ 秋おこつばさ/ococo.
- でびえんというかあろまちゃんがメイキングドラマを考えてネコ姐さんと出会う回。頭身が低いかわいい絵柄でギャグを交えつつも、気弱なあろまちゃん←無条件の信頼を寄せるみかんちゃんがしっかり描かれていてエモい百合になってしまっている。ノリがアニメそのものだし、起承転結がしっかりしててめちゃ面白い。ファルルオチ、ソラミとすれ違うエピローグもいい。1期の頃にプリパラに出入りしていたであろうあろまげどんのね、こういうお話がみたかったんですよ……このままOVA化してほしいマジで。ゲストの持ち物検査の話も優しい人ばかりのパプリカ学園って感じでいい。アメリカン・フットボールのシーンが狂った時のプリパラって感じで最高（一箇所だけ・がなくアメリカンフットボールになってるのに気づいた）。
あてのないはなし horogoru/ヌマル
- プリパラ外シオみれデート回。この人の描くシオみれは素晴らしい（既刊含め三冊全部持ってます）……お話にもやっとした感じの絵柄？がマッチしてると思う。相手に似てる人形を買うシーンやばくないですか？百合か？は？ってなってたら手を繋いでRealize!で絶叫してしまった。全体を通して物理的な声と心の声が別の吹き出しに書かれてて、二人の言動とその時の気持ちがわかるようになっています。どこまで素直になっていいのか？相手がどう感じるだろう？と悩んだり同じ気持ちだった時に喜んだり……そして最後はちゃんと素直な気持ちをちゃんとお互いに伝えるんですね。こういうのがエモい百合なんだよな（ろくろ）。ありがとうございます。
元気で明るい曲 もけお/ビショビショスマイル
- イラスト本。イラスト本はあまり買わないんですけど、まず表紙がまさに元気で明るい曲という感じで思わず手にとってしまった。ごちゃごちゃしてるけど統一感があるし何より色使いがすきだし絵がめっちゃ上手い。中野高円寺阿佐ヶ谷って感じ？なのかしら。ユリカ様（！？）とあろまちゃんの絵（これが買う決め手になった）、お腹が空いたみかんちゃんがあろまちゃんを××してしまう漫画、ミニあろま達が神々しい巨大みかん像を建立する漫画が特にすき。違反チケットでノンブルふってるのもいい……RLのみんなもいる！中盤にちょろっと描かれているのはマイキャラちゃん？
10年後の一番星,コピ本 あかりん殿下/変色だんご。
- アイドルを辞めたみれぃとアイドルを続けたらぁらのお話。クラシカルならぁみれを感じさせる設定でありながら、最後はみれぃが（勝手に自己完結して）前に進んでらぁらが諦めるという……一歩を踏み出せないらぁら、難聴系主人公みたいに勝手に納得するみれぃ……これですよ、これ！コピ本のWEB再録②もやばくないですか？白背景エモ百合ポエムという概念が質量を持って襲い掛かってきてませんか？は？
芸カ9とGLF16の本の感想も書きたいけどなかなか手がつけられない…… というか同人誌の感想を書くなんていう残酷なことをしてもいいのだろうか？という不安から恐らく作者に見つからないであろうブログに書いている。
未来向きの今
呪いの言葉かよ
プリパラのアイテム画像をBoVWで分類する
私はプリパラのことがとても好きだしプリパラのことを考えている時だけ”生きている”という感じがします。 皆さんもきっとそうだろうと思います。
しかしパソ君にとってはどうでしょうか？ パソ君は我々の生活になくてはならないものですが、パソ君はプリパラのことを”理解”しているでしょうか？ このパソ君の無機質な表情からはとてもそうは思えません。
パソ君にはプリパラを理解してもらう必要があります。 これはプリパラが現実を侵食するために絶対に必要なことです。
なんと幸いなことにプリパラの公式サイトにはプリパラのアイテム画像がしかもメタ情報付きでアップロードされています。
パソ君を調教するのにこいつを使わない手はないでしょう。
ということでパソ君にプリパラのアイテム画像を｛アゲアゲアイテム、ヘアアクセ、トップス、ボトムズ、ワンピース、シューズ｝の6種類に分類させることで、パソ君のプリパラに対する理解を深めてもらいました。 （識別率の絶対値よりも、プリパラの画像を認識する際の傾向を調査するのが本当にやりたいことだったりします。また、後で述べますが重大な実験上の誤りを見逃しています。）
今回採用したのはBag of visual words (BoVW)という方法です。
画像は生のRGBであれば各ピクセルに対し8*3=24bitの情報があり、画像全体の情報はそれのピクセル数倍になります。 冗長な情報が多く含まれているので生の画素値をそのまま使うのは最悪のアプローチで（あると普通は思いま）す。（CNNはそこを上手くやってしまうわけですが、それは置いておいて）
そこで画像を認識する際によく使われ（ていた）るのが局所特徴量というやつです。 局所特徴量の抽出はまずdetectorで特徴点を検出し、次に各特徴点に対してdescriptorで特徴ベクトルを算出するという2つのステップからなります。 局所特徴量はある画像に対して、特徴点と特徴ベクトルのペアの集合として定義されます。
ここで第一に問題になるのが、画像によって抽出される特徴点の数が違うということです。 画像ごとに局所特徴量の次元数が変わるので、ここままではよく使われている機械学習の手法をそのまま使うことができません。
第二の問題は、異なる画像に似たような特徴ベクトルが多く含まれている場合があるということです。 例えば、他の画像と全く違う局所特徴量が1つ抽出されたとしても、他の画像とよく似ている局所特徴量が10個抽出されたのでは、画像をパソ君に区別させるのは難しいです。
（私の理解では）以上2つの問題を解決したのが、BoVWというやつです。
まず、学習データから局所特徴量を抽出します。 そして、特徴ベクトルに対しクラスタリングを行いK個のクラスタに分類します。 最後に各画像の各特徴ベクトルに対し、クラスタ重心との最近傍探索による投票を行い、K次元のヒストグラムを作成します。 このK次元のヒストグラムを画像の特徴量とするのです。
第一の問題は全ての画像が同じK次元の特徴量を持つので解決しますし、第二の問題も多くの画像に含まれる似た特徴ベクトルは同じクラスタに潰れるので解決します。
解説はググれば色々でてきます。
ノウハウに関しては以下の記事が大変参考になりました。
何をしたのか
BoVWでプリパラのアイテム画像を｛アゲアゲアイテム、ヘアアクセ、トップス、ボトムス、ワンピース、シューズ｝の6種類に分類した。 局所特徴量の種類とクラスタ数を変えて識別率を比較した。 データは公式サイトにあるラベル付き画像データを使った。
この実験で明らかにしたかったこと
- どれくらい性能でるのか？
- どの局所特徴量が性能が良いのか？
- クラスタリングの次元数を上げるとどれくらい性能は上がるのか？
実験環境
ソフト：Windows、Python 3.4.1、OpenCV 3.0RC（特徴抽出からヒストグラム作成まで）、scikit-learn（SVMによる識別）
データ：2015/07/25時点にプリパラ公式サイトにアップロードされていたアイテム画像（ アゲアゲアイテム：127、 ヘアアクセ：247、 トップス：328、 ボトムス：288、 ワンピース：296、 シューズ：473、 合計：1759 ）
画像サイズ：｛256×256｝（元サイズは｛512×512｝）
局所特徴量：｛KAZE, AKAZE, ORB, BRISK｝
クラスタリングの次元数（K）：｛128, 256, 512, 1024, 2048, 4096, 8192｝
識別器：SVM（RBFカーネル、C=1、Γ=8、1対1分類器による6クラス識別）
評価方法：5 fold cross validationを行う。5つのfoldの識別率の平均値で評価する。
表1：実験に使った局所特徴量
| 特徴量名 | dense sampling方法 | 種類 | 距離尺度 | OpenCVのコンストラクタ(特徴量名_create())のパラメータ | 引用 |
|---|---|---|---|---|---|
| KAZE | 特徴点検出の閾値を0にする | 実数 | ユークリッド距離 | threshold=0.0 | KAZE Features. Pablo F. Alcantarilla, Adrien Bartoli and Andrew J. Davison. In European Conference on Computer Vision (ECCV), Fiorenze, Italy, October 2012. |
| AKAZE | 特徴点検出の閾値を0にする | バイナリ | ハミング距離 | threshold=0.0 | Fast Explicit Diffusion for Accelerated Features in Nonlinear Scale Spaces. Pablo F. Alcantarilla, Jesús Nuevo and Adrien Bartoli. In British Machine Vision Conference (BMVC), Bristol, UK, September 2013. |
| BRISK | 特徴点検出を使わず、背景以外から4点置きにサンプリング | バイナリ | ハミング距離 | thresh=0.0 | Stefan Leutenegger, Margarita Chli and Roland Siegwart: BRISK: Binary Robust Invariant Scalable Keypoints. ICCV 2011: 2548-2555. |
| ORB | 特徴点検出を使わず、背景以外から4点置きにサンプリング | バイナリ | ハミング距離 | edgeThreshold=11, patchSize=11 | Ethan Rublee, Vincent Rabaud, Kurt Konolige, Gary R. Bradski: ORB: An efficient alternative to SIFT or SURF. ICCV 2011: 2564-2571. |
上記記事のように、detectorによる特徴点検出を使わず密に等間隔にサンプリングした点に対して特徴ベクトルの算出を行うdense samplingを行いたかったのですが、KAZEとAKAZEではそれができませんでした。（OpenCVのドキュメントにもできないとはっきり書いてありました） そこでKAZEとAKAZEについては特徴点検出の閾値を0にし可能な限り多く特徴点を算出することでdense samplingを代用しました。 こんな感じです。（特徴ベクトルを算出した点を黄緑で示しています）
図1：AKAZEで検出された特徴点の例
ORBとBRISKについてはdetectorを使わず背景を除く4ピクセル置きに特徴ベクトルを算出しました。dense samplingを行ってくれる関数がOpenCV3.0からなくなってしまっていたので、この処理は自分で実装しました。本来はスケールも変えたほうがいいはずですが、自分で書くのは大変だったのでやめました。（C++のオレオレ実装はネットで見つけましたが今回はPythonなので使えませんでした。） こんな感じです。（特徴ベクトルを算出した点を黄緑で示しています）
図2：密にとった特徴点の例
特徴点が多い＝特徴ベクトルが多いほど、特徴ベクトルのクラスタリングの信頼度が高まるはずです。（こんな表現をしたら確率統計マンに殺害されそうだけど上手い言い方がわからない） より正しくクラスタリングの重心が求まれば、より高い識別率が期待できるでしょう。 画像をみてわかるように今回使っている特徴点の数はORBやBRISKのほうが遥かに多いです。 仮にKAZEやAKAZEのdescriptorの性能が良かったとしても、特徴ベクトルの数が全く違うのでORBやBRISKのほうが性能が良いではないかと予想できます。
実験結果
予想通り、BRISKとORBのほうがKAZEとAKAZEより識別率が良いという結果になりました。 BRISK＞ORB＞AKAZE＞KAZEという順序になりました。
表2：識別率とクラスタ数の関係
| クラスタ数 | KAZE | AKAZE | BRISK | ORB |
|---|---|---|---|---|
| 128 | 0.701 | 0.739 | 0.842 | 0.814 |
| 256 | 0.720 | 0.753 | 0.865 | 0.843 |
| 512 | 0.766 | 0.793 | 0.901 | 0.868 |
| 1024 | 0.822 | 0.839 | 0.914 | 0.883 |
| 2048 | 0.835 | 0.850 | 0.925 | 0.898 |
| 4096 | 0.866 | 0.860 | 0.937 | 0.906 |
| 8192 | 0.861 | 0.876 | 0.936 | 0.916 |
図3：識別率とクラスタ数の関係
考察
この実験から言えるのはdescriptorの性能はKAZEよりAKAZEのほうが良く、またORBよりBRISKのほうが良いということです。 特徴点の取り方が全然違うので識別率の絶対値が高いからといって一概にORBやBRISKのほうがKAZEやAKAZEより良いとはいえません。
手を動かす前はAKAZEが一番新しいのでAKAZEが一番性能がいいのではないかと思っていましたが、そうはなりませんでした。
KAZEとAKAZEの特徴は特徴点検出にあると私は考えています。 SHIFTやSURF等（今回は比較していませんが）のdetectorはDoG（Difference of Gaussian）ですが、これはエッジまでもぼかしてしまい微細な特徴点が検出できないという問題がありました。 ORBとBRISKのdetectorはFASTベースで、これは速いのですが誤検出が多い印象があります。 一方、KAZEとAKAZEではエッジを選択的に残すnon-linear diffusionを使い、微細な特徴点も比較的正確に検出することができます。（のようなことが論文に書いてあります、たぶん）
今回はdetectorを活用しきれなかったのでAKAZEの性能が伸びなかったのではないかと考えています。
不正行為
この実験では最悪なズルをしていて、特徴抽出を行う際の学習データにテストデータが含まれています。 全データから抽出した全特徴ベクトルに対してK-meansをかけ重心を求め、その重心とマッチングを行うことで全データのヒストグラムを作成しています。 （認識時はcross validationをしています。あるfold以外のヒストグラムを使って識別器を学習し、そのfoldでテストしています）
正しくは各foldに対し、それ以外のfoldに含まれる特徴ベクトルに対してクラスタリングを行いクラスタ重心を求め、そのクラスタ重心に対してヒストグラムを作成すべきです。 本当は、学習時にはテスト時に入ってくるデータの特徴ベクトルなんて知らないはずですから。
全部終わる頃になって、これはダメダメダッメーズルズルズッル－だよと気づきましたが、ただでさえクラスタリングにめっちゃ時間がかかったし、foldごとにクラスタリングするとなると約fold数倍の時間がかかるはずなのでめんどくさくなってやめました。 これは本当はやってはいけないことですが、これは実験レポートでも論文でもないし、なんとなく傾向がわかったのでまあよしとしましょう。
何が言いたいかというと、識別率の絶対値にあまり意味はありません。識別率は全体的に本来の値より高くなっていると思われます。
この実験で明らかになったこと
- どれくらい性能でるのか？→90%くらい。しかしズルしてるので本当はもっと低いはず。
- どの局所特徴量が性能が良いのか？→今回比較した中ではBRISKが一番よかった。descriptorの性能もあるが密に特徴ベクトルを算出したことの影響も大きいはず。
- クラスタリングの次元数を上げるとどれくらい性能は上がるのか？→次元数2倍で識別率は数%上がる
結論
パソ君の調教にある程度成功した。パソ君がプリパラの素晴らしさを理解し始めたのは大変良いことだ。
その他
Kを大きくするとクラスタリングはめちゃくちゃ時間がかかりました。K=8192でcore i7のコアを全部使って半日くらいかかったと思います。
また、ユークリッド距離しかサポートしていないK-meansライブラリを使ってバイナリ特徴量をクラスタリングする際はちょっと注意しないといけないのですが、それについては気が向いたら書きます。
余談ですが、プリパラ公式サイトのHTMLはニンゲンの手のぬくもりが感じられるHTMLで、本来”トップス”が入るべきところに”ラブリー”が入っていたりして全自動でスクレイピングできませんでした。 どうも世界を支配するプリパラの”システム”が作っているのではないようです。
プリパラで手持ちのプリチケからレア度ごと（N、R、SR、PR）の出現確率を推定する
このエントリは2014 1stライブから2014-2015 3rdライブまでのプリパラのシステムを前提として書いています。 というか2014-2015 3rdが稼働中に公開するつもりだったのが間に合いませんでした。 2015/4/2に稼働した2015 1stライブからはドリチケとPRよりも珍しいMR（ミステリーレア）が導入され少しシステムが変わっていますが、これらには対応しておりません。そういうものだと思って下さい。
このエントリでは理系学部レベルの数学を聞いたことがあることを前提としており、そういった箇所に関しては説明を省いている場合があります。 一方、学部では習わなそうな箇所やマニアックな箇所については説明を入れているつもりです。
また、
重要箇所
重要箇所についてはこのような表記を使っています。 markdownの引用表記です。 この表記を追うだけでなんとなくわかるように書いたつもりです。
はじめに
みなさんプリパラしていますか？ 私は沢山プリパラしているし、プリパラをしているとどんどん気持ちよくなってきます。
前回はコーデのレア度の出現確率を設定し、 シミュレーションを行うことでプリチケをコンプリートするために必要なプレイ回数の平均値を求めました。
しかし前回はほぼカンでレア度の出現確率を設定していました。
皆さんは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()
プリパラのコーデをコンプするための期待値
プリパラにハマっています。課金が前提のゲームに熱中するのは初めてです。
12月末にプリパラ筐体デビューを果たしたのですが、ファイルやらグミやら付録目当ての関連書籍やらを含めると既に万単位でお金をつぎ込んでいる気がします。
射幸心を煽るというやつなのでしょうか。コーデを全身揃えたくなって次から次へとは百円玉を筐体に入れてしまいます。その時の自分は何も考えていないのでとても幸せです。
このまま無計画に百円玉を入れていてはいつか破産するかもしれません。プリパラのために破産するなら本望です。しかし、どれくらい百円玉を入れればコーデが揃うのか見通しを立てておいたほうが借金の計画も立てやすいでしょう。
というわけでプリパラのコーデを揃えるには何枚百円玉を入れればいいのか計算してみることにしました。
http://members3.jcom.home.ne.jp/zakii/enumeration/30_capsuletoy.htm
http://doryokujin.hatenablog.jp/entry/2012/05/09/034209
「ガチャ コンプ 期待値」などでググるとこういう記事がヒットします。こういう記事によるとこれはCoupon Collector's Problem と呼ばれている問題であり、ガチャをコンプする期待値はガチャの種類をn個とすると
だそうです。これで終わりです。ぱちぱちぱち。
ここで東堂シオンが「一件落着！」と野太い声で叫び、すかさずレズの南委員長から「私の計算では云々」ツッコミが入ります。
上の式はガチャがn種類で、その全てが等確率で出現する場合の期待値です。
プリパラはシチュエーションがだいぶ違います。
プリパラ概論
プリパラのシステムについて簡単に説明しましょう。
プリチケをプレイすると貰えるキラキラした可愛らしい代物、プリチケはパキるために上部と下部に分かれています。（「トモチケをパキる」は本当に美しい日本語だと思います。）
プリチケの部品は大きく分けてマイチケ（下部）、トモチケ（上部）、ガチャチケ（上部）の3つがあります。
マイチケは百円玉を入れると必ず貰えます。
しかし、トモチケはライブをした時だけ、ガチャチケは逆にライブをしなかった時にしか貰えません。
さらにコレクションの対象となるのは
マイチケ（下部）に印刷されるコーデ
ガチャチケ（上部）に印刷されるヘアアクセ
ガチャチケ（上部）に印刷されるアゲアゲアイテム
の3つです。ややこしいですね。ガチャチケにはさらに二つの種類があるのです。
コーデとヘアアクセは月替りで、アゲアゲアイテムは数ヶ月のシーズンごとにアップデートされるようです。
トモチケは自分が印刷されているものであり、アイテムが貰えるわけではないのでコレクションの対象にはならないでしょう。
ライブをする、しないの二つの場合にわけて何が貰えるのか整理します。
1.ライブをする
ライブ後にプリズムストーンに行くことができトモチケとマイチケがもらえる。
トモチケは自分が印刷されているもので、友達と交換して遊ぶ。
プリズムストーンに行くとコーデ（トップス、ワンピース、ボトムス、シューズ）がランダムで二種類表示され、そこから好きな片方を選んでマイチケにできる。
コーデにはレア度がありN、R、SR、PRの順で出にくくなる。
2.ライブをしない
ライブ後にさらに百円玉を入れる、もしくはライブをせずに直接プリズムストーンに行くと、ガチャチケとマイチケがもらえる。
マイチケに関してはライブ時と同じ。
ガチャチケは 所謂普通のガチャで、アゲアゲアイテムかヘアアクセが貰える。
アゲアゲアイテムとヘアアクセにはレア度がなく、変わりにレベルがある。レベルが高いと出にくいのか、それとは別にアイテム別に出にくさが設定されているのか不明。
ややこしいので後者に絞って考えます。レオナの気持ちになってめがにぃに貢ぎまくることを想定します。
問題設定
想定するシチュエーションを整理すると
・百円玉を一枚入れるとマイチケとガチャチケが一つずつ貰える
・マイチケはランダムに出現した二つのコーデから選択する自由がある。出現率はレア度によって変わる
・ガチャチケは選択の自由がない。出現率がどうやって決まるのか不明
・全身コーデを揃えるにはヘアアクセをガチャチケとして、それ以外をマイチケとして入手する必要がある
となります。
冒頭で紹介したシンプルな設定でも結構ややこしいので、この場合の期待値をきちんと数学を使って求めるのは難しそうです。
というか不明点が多々あります。
そこで不明点を妄想で補ったシミュレーションを行って、コンプまでの平均プレイ回数を計算してみます。
コンプするのは現在行われている2015 3rdライブ 2月のアイテムのうち、見た目に関わるコーデとヘアアクセを対象とします。
また、ヘアアクセのレベルは問わないことにします。
アゲアゲアイテムはシーズン替りだし、見た目に反映されないので血眼になって集める必要はないんじゃないかと思います。
シミュレーションの方針
シミュレーションの方向性は２つあると思います。
・コーデが2つでてきて片方を選びプリチケ化するというプロセスをそのままシミュレーションする
・2つから片方を選ぶプロセスを無視し、最終的にプリチケ化するところだけをシミュレーションする
今回は素直に前者を採用しました。
前者を採用するとなると、2つのコーデから片方を選択する方針を決める必要があります。
選択方針としては、「持っていないほうを選ぶ。両方とも持っていない場合にはレア度が高いほうを選ぶ。両方とも持っておらずさらにレア度が等しい場合は、ランダムに片方を選ぶ」という方針を採用します。
実際には、持っていないPRとSRが出てきてSRの全身コーデを揃えたいからSRを選ぶということも考えられるでしょうが、今回はシンプルな方針を採用しました。
わかっていること
http://pripara.jp/item/pdf/2015_3rd02_itemlist.pdf
http://pripara.jp/item/2015_3rd_01.html
プリパラの公式サイトをみれば
・2015 3rdライブ 2月のレア度別のコーデの総数
N：6種類
R：6種類
SR ：7種類
PR：6種類
・2015 3rdライブ のアゲアゲアイテムの総数：9種
・2015 3rdライブ 2月のヘアアクセの総数：12種
がわかります
わからないこと
・レア度とコーデの出現率の関係
仕方がないので手持ちのプリチケから最尤推定します。最尤推定という言葉がかっこいいので使いますがただ数えてるだけです。
N：24枚（スターランク差分：1、重複：2を含む）
R：19枚（スターランク差分：1）
SR ：10枚（スターランク差分：1）
PR：7枚（スターランク差分：1）
月次ライブで獲得した手持ちのマイチケを数えると以上でした。
サイリウムコーデやパラダイスコーデ等、貰えるコーデが確定しているライブで手に入れたプリチケは無視しています。
しかし、コーデは同時に２つでてきます。実際にプリズムストーンに出現したコーデ数は手元にあるマイチケの倍であるはずです。
今回のシミュレーションを行う上では、マイチケにする確率ではなく、プリズムストーンにコーデが出現する確率が必要になります。
このままではプリズムストーンにコーデが出現する確率は推定できません。
そこで、基本的にレア度が高いほうを選んでいるはずという仮定とあいまいな記憶にもとづいて、マイチケにしなかったコーデの枚数を補完します。
| レア度 | 手持ちのマイチケ枚数 | プリチケ化したNと同時に出現した気がする回数 | プリチケ化したRと同時に出現した気がする回数 | プリチケ化したSRと同時に出現した気がする回数 | プリチケ化したPRと同時に出現した気がする回数 | プリズムストーンに出現した想定回数（合計） | `推定出現確率 |
|---|---|---|---|---|---|---|---|
| N | 24 | 24 | 10 | 6 | 3 | 67 | 67/120 |
| R | 19 | 0 | 9 | 3 | 3 | 34 | 34/120 |
| SR | 10 | 0 | 0 | 1 | 1 | 12 | 12/120 |
| PR | 7 | 0 | 0 | 0 | 0 | 7 | 7/120 |
水増しした頻度に基づきレア度ごとの出現確率を推定しています。
さらにシミュレーションを行う上で、同じレア度のコーデは同じ確率で出ることにします。これはそれほどまずい仮定ではないと思います。
私は出現したコーデが両方とも持っていた場合に、スターランクアップあるいは手持ちと重複するコーデを選んでいます。
このへんも考慮したほうが良さそうですが、今回は無視します。
・ガチャチケの出現確率
全然情報がないので、一様分布を仮定します。つまり、どのアゲアゲアイテムもガチャチケも同じ確率で出るということにします。レア度が高いコーデのヘアアクセが出にくいような気がするのですが、そういうところを推測でやるのはやめます。
ここまで情報を集めるとコードが書けます。
長いので最後に載せます。
シミュレーション結果
コーデコンプ、ヘアアクセコンプ、コーデとヘアアクセコンプ、PRコーデ（ヘアアクセ含む）の片方をコンプ、PRの両方をコンプ
について、100000回試行してコンプまでの平均プレイ回数と標準偏差を出しました。
PRコーデについて計算したのはやっぱりPRコーデを揃えたいからです。今月のPRはドロシーとレオナのノクターンスカイアイドルコーデです。とても可愛いのでとてもよいと思います。
| コンプ対象 | コンプまでの平均プレイ回数 | 標準偏差 |
|---|---|---|
| コーデとヘアアクセ | 136.8 | 55.62 |
| コーデ | 134.7 | 56.93 |
| ヘアアクセ | 65.0 | 24.91 |
| PRの片方 | 64.3 | 33.34 |
| PR両方 | 126.0 | 60.65 |
標準偏差がデカいですね。
百円玉が140枚くらいあれば2月のコーデが揃えられるみたいです。もっとあくどいかと思いましたが、こんなもんなんですね。しかし、女児がコンプを目指すのは厳しいと思います。
最近ドロシーのノクターンスカイアイドルコーデが揃ったのですが、この結果を見て自分は運がよかったほうなのかなと思いました。
この結果はわからないことを推測して出したので大体あってる保証すらもありません、一応。
コード
最後にコードです。
pythonですが、無駄にpython3です。
#! coding:UTF-8 import random from itertools import accumulate class coordInfo: def __init__(self, number, rate): self.number = number self.rate = rate #コーデはレア度によって重みをつける def createCoordDistribution(coords): z = 0.0 dist = [] #レア度が同じコーデは全て同じ確率で出る for coord in coords: dist += [coord.rate] * coord.number z += coord.rate * coord.number #正規化 dist = [d/z for d in dist] #累積分布の計算 accum = list(accumulate(dist)) return dist, accum #ガチャは一様分布 def createGachaDistribution(ageage, hairAcce): dist = [1] * (ageage + hairAcce) z = sum(dist) dist = [d/z for d in dist] accum = list(accumulate(dist)) return dist, accum #離散累積分布から(0,1)の一様分布を使ってサンプリングする def sampling(accum): r = random.random() for i, a in enumerate(accum): if r < a: return i def isCompleted(gotten): return all(gotten) #後半acceNum個をヘアアクセに、残りの前半をアゲアゲアイテムに割り当てる def isCompletedHairAcce(gotten, acceNum): return isCompleted(gotten[-acceNum:]) #コーデは出現率の高い順から並んでいる。 #PR最後３つをドロシーのノクターンスカイアイドルコーデに、PR残りの3つをレオナに割り当てる。 #ヘアアクセも動揺。 def isCompletedPrareOne(coordGotten, gachaGotten): dorothy = isCompleted(coordGotten[-3:]) and gachaGotten[-1] reona = isCompleted(coordGotten[-6:-3]) and gachaGotten[-2] return dorothy or reona def isCompletedPrareTwo(coordGotten, gachaGotten): dorothy = isCompleted(coordGotten[-3:]) and gachaGotten[-1] reona = isCompleted(coordGotten[-6:-3]) and gachaGotten[-2] return dorothy and reona def updateCoordGotten(coordAccum, coordGotten): s1 = sampling(coordAccum) s2 = sampling(coordAccum) #違うコーデが出るまで粘る while s1 == s2: s2 = sampling(coordAccum) #持っていないほうをとる if coordGotten[s1] and not coordGotten[s2]: coordGotten[s2] = True elif not coordGotten[s1] and coordGotten[s2]: coordGotten[s1] = True elif not coordGotten[s1] and not coordGotten[s2]: #両方持っていない場合には確率が低い=レア度が高いほうをとる if coordDist[s1] < coordDist[s2]: coordGotten[s1] = True else: coordGotten[s2] = True else: #両方持っている場合は何もしない pass def updateGachaGotten(gachaAccum, gachaGotten): s = sampling(gachaAccum) gachaGotten[s] = True def calc(values): average = sum(values) / len(values) stdev = 0 for v in values: stdev += (v - average)**2 stdev = (stdev / len(values)) ** 0.5 return average, stdev def printResult(message, average, stdev=None): if stdev is None: stdev = average[1] average = average[0] print("{0:<23}".format(message), end="") print("{0:.1f}".format(average), "{0:.2f}".format(stdev)) if __name__ == '__main__': #手元にあるプリチケ枚数 + Nと同時に出た気がする回数 + Rと～ + SRと～ + PRと～ pRareNum = 7 + 0 + 0 + 0 + 0 #スターランク差分1、重複0 sRareNum = 10 + 0 + 0 + 1 + 1 #スターランク差分1、重複0 rRareNum = 19 + 0 + 9 + 3 + 3 #スターランク差分1、重複0 nRareNum = 24 + 24 + 10 + 6 + 3 #スターランク差分1、重複2 pRare = coordInfo(6, pRareNum) sRare = coordInfo(7, sRareNum) rRare = coordInfo(6, rRareNum) nRare = coordInfo(6, nRareNum) ageage = 9 hairAcce = 12 expRun = 100000 coordinates = [nRare, rRare, sRare, pRare] coordDist, coordAccum = createCoordDistribution(coordinates) gachaDist, gachaAccum = createGachaDistribution(ageage, hairAcce) coordComp, gachaComp, allComp, pRareOneComp, pRareTwoComp = [], [], [], [], [] for i in range(expRun): run, coordRun, gachaRun, pRareOneRun, pRareTwoRun = [0] * 5 coordFlag, gachaFlag, pRareOneFlag, pRareTwoFlag = [True] * 4 coordGotten = [False] * len(coordAccum) gachaGotten = [False] * len(gachaAccum) while not isCompleted(coordGotten) or not isCompletedHairAcce(gachaGotten, hairAcce): updateCoordGotten(coordAccum, coordGotten) updateGachaGotten(gachaAccum, gachaGotten) run += 1 if coordFlag and isCompleted(coordGotten): coordRun = run coordFlag = False if gachaFlag and isCompletedHairAcce(gachaGotten, hairAcce): gachaRun = run gachaFlag = False #PRのどちらか片方が揃ったら if pRareOneFlag and isCompletedPrareOne(coordGotten, gachaGotten): pRareOneRun = run pRareOneFlag = False #PRのどちらか両方が揃ったら if pRareTwoFlag and isCompletedPrareTwo(coordGotten, gachaGotten): pRareTwoRun = run pRareTwoFlag = False allComp.append(run) coordComp.append(coordRun) gachaComp.append(gachaRun) pRareOneComp.append(pRareOneRun) pRareTwoComp.append(pRareTwoRun) printResult("all items", calc(allComp)) printResult("all coodinates", calc(coordComp)) printResult("all hair accessories", calc(gachaComp)) print() printResult("one PR", calc(pRareOneComp)) printResult("two PRs", calc(pRareTwoComp))
