最長増加部分列(LIS)の期待値についての考察(その2：本題、LISの期待値)
こんにちは、とぼです。今回の記事は、前回記事が長くなりすぎた最長増加部分列(LIS)問題の期待値に関する考察の続き...もとい本題です。
LISについては前回の記事でめっちゃ詳しく書いたので、(もしご存じでない方がいれば)(むしろプロがこの記事を見かけてくれたら、むしろさすまた投げてください。)
yadrfu.hatenablog.com では早速、LISの期待値に関する考察をします！！
本題：LISの期待値
丁寧に書いていたら、多分前書きの方が長くなったんですが、ここから本題(?)。Twitterで「整数1~nのランダムな並び替えのLISの期待値」はnに対してどのような挙動をとるのか、という話題が上がっておりまして、大変興味深かったので、少し実験をしてみました。
少し実験してみましょう。以下、1～nのLISの期待値をEL(n)と書くことにします。
（１）n=1のとき、EL(1)=1は自明。
（２）n=2のとき、数列が{1,2}ならLIS=2,{2,1}ならLIS=1なので、EL(2)=(1+2)/2=3/2。
（３）n=3のとき、数列が{1,2,3}でLIS=3,数列が{3,2,1}でLIS=1,それ以外ではLIS=2と なるのでEL(3)=(1+3+2*4)/6=2。
（４）n=4のときは...,うーん手計算では難しい(面倒)！と、いうわけでネクパ
(next_permutation。詳しくはProject Eulerの記事その３
を参照。)を使って、各組合せに対しLISを求め、（脳筋全探索で）期待値を求めてみましょう！！
以下のような脳死コードを書きました。
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <cmath>
#include <set>
#include <map>
#include <string>
#include <numeric>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <iterator>
#include <cstring>
using namespace std;
#define rep(i, a) for (int i = 0; i < a; i++)
//ユークリッドの互除法
ll gcd(ll x, ll y) {
if (x < y) swap(x, y);
if (y == 0) return x;
return gcd(y, x%y);
}
//最長増加部分列の長さの計算
int LIS(vector<int> V) {
int DP[11] = { 0 };
fill(DP, DP + V.size(), INF);
rep(i, V.size()) {
*lower_bound(DP, DP + V.size(), V[i]) = V[i];
}
return lower_bound(DP, DP + V.size(), INF) - DP;
}
int main()
{
rep(i, 10) {
vector<int> V(i+1);
iota(V.begin(), V.end(), 0);
int nEx = 0; //期待値分子(EL(n)の分子)
do {
nEx+=LIS(V); //期待値分子の加算
} while (next_permutation(V.begin(), V.end()));
int All=Kaijou(i + 1); //期待値分母=n!
int t = gcd(All, nEx); //最大公約数求めて約分
nEx /= t; All /= t; //約分
cout<<"n:"<<i+1<<", Ex:"<<nEx<<"/"<<All<<endl; //EL(n)を分数表示
}
return 0;
}
実行結果を以下に示します。
以上のようにn=1~10までのEL(n)を導出できました。これ以上は実行時間が長くてダメですね(まあ、n=12くらいまでは数分レベル待てばできるとは思うが)。
で、法則を見つけようとか、nと増加量の関係とか考えようとしましたが脳死していたのでよくわかりませんでした。(2016がいてちょっと喜ぶくらい)
そこで "" つ「オンライン整数列大辞典(OEIS)」 !!!! ""
折角分数表示したので、各EL(n)の分子{1,3,2,29,67,2261,499,7601,163673,3146141}をOEISにぶち込みます。
で、見事ヒットします。OEISマジで整数列ならなんでもあってヤバイ（語彙力）。
さらに、よく見ると...
...ん？どこかで見覚えのある分数列が...
あっっっ！！！！
まさにEL(n)ですやん...。OEISやべぇな。（語彙力)
というわけで、手元では（現状のアルゴリズムでは）計算困難なn=21までのEL(n)の値を得ることができました。
これをExcelにぶち込んで、nとEL(n)の関係を表すグラフを描いてみました。下図の青実線がnとEL(n)の関係を示すグラフで、赤点線が、それをもとに推定させた近似曲線です。近似曲線のパターンを「線形」、「対数形」、「累乗形」と試したところ一番フィットした「累乗形」を採用しています。下図に示すように
EL(n) = 0.9972n ⁰・⁶³⁴⁹ ≒ n ⁰・⁶³⁴⁹
の近似曲線を得ることができました。
しかしながら、OEISにさえFOMULAが載ってなかったので、正確なEL(n)の一般解を求めるのは困難と思われます。
まあ今回は近似曲線を求めることができたので個人的には満足です。
現場からは以上になります。
