コンビネーションとか二項分布とか三分探索とか

まえがき

ちょっとハースストーンで遊んでいたら検証したいことが出てきたので、尤度を計算したよ。

そのときのコードです。何やってたかしりたかったらこっちも見てね ヨグ＝サロンは果たして自分のヒーローの呪文を多く選ぶのか？ - ふらふらHearthstone

コンビネーションの計算とか二項分布の計算とか三分探索とかやってます。

コンビネーションのdoubleバージョンに関してはある程度大きくなると誤差がでてきます。誤差保証はしていないので、注意してください。

使いたかったらご自由にコピペしてどうぞ。

以下スパゲッティーなコード

#include <iostream>
#include <vector>
#include <map>
#include <algorithm>
#include <cassert>
#include <functional>
#include <iomanip>
using namespace std;
class combination
{
    map<pair<int,int>,long long> memo;
    map<pair<int,int>,double> memo_f;
public:
    long long combi(int n, int x) {
        assert(n>=0);
        assert(x>=0&&x<=n);
        if(n==x)return 1;
        else if(x==0) return 1;
        else{
            auto f=memo.find(make_pair(n,x));
            if(f!=memo.end()){
                return (*f).second;
            }else{
                auto res= combi(n-1, x)+combi(n-1, x-1);
                assert(res>=0);
                memo[make_pair(n,x)]=res;
                return res;
            }
        }
    }
    double combi_f(int n, int x) {
        assert(n>=0);
        assert(x>=0&&x<=n);
        if(n==x)return 1;
        else if(x==0) return 1;
        else{
            auto f=memo_f.find(make_pair(n,x));
            if(f!=memo_f.end()){
                return (*f).second;
            }else{
                auto res= combi_f(n-1, x)+combi_f(n-1, x-1);
                assert(res>=0);
                memo_f[make_pair(n,x)]=res;
                return res;
            }
        }
    }
};
long long combi(long long n,long long x){
    static combination c;
    return c.combi(n,x);
}
double combi_f(long long n,long long x){
    static combination c;
    return c.combi_f(n,x);
}
// 確率pのベルヌーイ試行をn回したときにx回出る確率
double nikou(int x,int n,double p){
    return combi_f(n,x)*pow(p,x)*pow(1-p,n-x);
}
struct Yogdata{int x;int n;int myhero;int all;};
double likelihood(vector<Yogdata> data,double mass){
    double ans=1;
    for(auto item:data){
        ans*=nikou(item.x,item.n,item.myhero*mass/(item.myhero*mass+item.all-item.myhero));
    }
    return ans;
}

template<class In,class Out>
void test(Out (*f)(In)){
    In x;
    while(true){
        cin>>x;
        cout<<f(x)<<endl;
    }
}
template<class In1,class In2,class Out>
void test(Out (*f)(In1,In2)){
    In1 x;
    In2 y;
    while(true){
        cin>>x>>y;
        cout<<f(x,y)<<endl;
    }
}
template<class In1,class In2,class In3,class Out>
void test(Out (*f)(In1,In2,In3)){
    In1 x;
    In2 y;
    In3 z;
    while(true){
        cin>>x>>y>>z;
        cout<<f(x,y,z)<<endl;
    }
}
template<class In1,class In2,class Out1,class Out2>
void test(Out1 (*f1)(In1,In2),Out2 (*f2)(In1,In2)){
    In1 x;
    In2 y;
    while(true){
        cin>>x>>y;
        cout<<setprecision(20);
        cout<<f1(x,y)<<"\t"<<f2(x,y)<<endl;
    }
}
double sanbutan(double bn,double ed,double f(double)){
    if(abs(bn-ed)<0.0000001)return bn;
    else{
        double le=(bn*2+ed)/3,ri=(bn+ed*2)/3;
        if(f(le)>=f(ri))return sanbutan(bn,ri,f);
        else return sanbutan(le,ed,f);
    }
}

vector<Yogdata> data;
double wrap(double mass){
    return likelihood(data,mass);
}
void run(){
    int N;
    cin>>N;
    data=vector<Yogdata>(N);
    for(int i=0;i<N;i++){
        int a,b,c,d;
        cin>>a>>b>>c>>d;
        data[i]={a,b,c,d};
    }
    for(int i=0;i<5;i++){
        double mass=i/1.0;
        cout<<mass<<"\t"<<likelihood(data,mass)<<endl;
    }
    cout<<sanbutan(0.0,5.0,wrap)<<endl;
}
int main() {
    run();
}