#include <iostream>
#include <cstdlib>
#include <map>
#include <tuple>
#include <iomanip>
#include <cmath>
#include <algorithm>
#include <cassert>
using namespace std;
class Hypothesis {
double weight;
map<tuple<int, int>, double> memo;
public:
Hypothesis(double weight) :
weight(weight)
{
}
double open_chance(int opened) {
int closed = 36 - opened;
double result = (double)closed / (weight * opened + closed);
assert(0 <= result && result <= 1);
return result;
}
double likelihood(int total, int opened) {
if (total == 0 && opened == 0) return 1;
if (total == 0 || opened == 0) return 0;
auto hit = memo.find(make_tuple(total, opened));
if (hit != memo.end()) return hit->second;
auto result =
likelihood(total - 1, opened) * (1 - open_chance(opened)) +
likelihood(total - 1, opened - 1) * open_chance(opened - 1);
memo[make_tuple(total, opened)] = result;
assert(0 <= result && result <= 1);
return result;
}
double log_likelihood(int total, int opened) {
return log(likelihood(total, opened));
}
};
int main()
{
struct input_type {
int opened;
int total;
};
input_type inputs[] = {
input_type{ 25, 53 },
input_type{ 28, 72 },
input_type{ 25, 72 },
input_type{ 29, 72 },
input_type{ 27, 69 }
};
int N = 50;
for (int i = -N; i <= N; ++i) {
double wight = pow(10, (double)i / N);
Hypothesis hypothesis(wight);
double sum = 0;
for (auto input : inputs) {
sum += hypothesis.log_likelihood(input.total, input.opened);
}
int barlength = (int)round(20 + sum);
if (barlength < 0) barlength = 0;
if (barlength > 20) barlength = 20;
string bar = string(barlength, ' ') + string(20 - barlength, '*') + "|";
cout << setw(12) << wight << "\t" << setw(12) << sum << "\t# " << bar << endl;
}
}
I2luY2x1ZGUgPGlvc3RyZWFtPgojaW5jbHVkZSA8Y3N0ZGxpYj4KI2luY2x1ZGUgPG1hcD4KI2luY2x1ZGUgPHR1cGxlPgojaW5jbHVkZSA8aW9tYW5pcD4KI2luY2x1ZGUgPGNtYXRoPgojaW5jbHVkZSA8YWxnb3JpdGhtPgojaW5jbHVkZSA8Y2Fzc2VydD4KCgp1c2luZyBuYW1lc3BhY2Ugc3RkOwoKY2xhc3MgSHlwb3RoZXNpcyB7Cglkb3VibGUgd2VpZ2h0OwoJbWFwPHR1cGxlPGludCwgaW50PiwgZG91YmxlPiBtZW1vOwoKcHVibGljOgoJSHlwb3RoZXNpcyhkb3VibGUgd2VpZ2h0KSA6CgkJd2VpZ2h0KHdlaWdodCkKCXsKCgl9CgoJZG91YmxlIG9wZW5fY2hhbmNlKGludCBvcGVuZWQpIHsKCQlpbnQgY2xvc2VkID0gMzYgLSBvcGVuZWQ7CgkJZG91YmxlIHJlc3VsdCA9IChkb3VibGUpY2xvc2VkIC8gKHdlaWdodCAqIG9wZW5lZCArIGNsb3NlZCk7CgkJYXNzZXJ0KDAgPD0gcmVzdWx0ICYmIHJlc3VsdCA8PSAxKTsKCQlyZXR1cm4gcmVzdWx0OwoJfQoKCWRvdWJsZSBsaWtlbGlob29kKGludCB0b3RhbCwgaW50IG9wZW5lZCkgewoJCWlmICh0b3RhbCA9PSAwICYmIG9wZW5lZCA9PSAwKSByZXR1cm4gMTsKCQlpZiAodG90YWwgPT0gMCB8fCBvcGVuZWQgPT0gMCkgcmV0dXJuIDA7CgoJCWF1dG8gaGl0ID0gbWVtby5maW5kKG1ha2VfdHVwbGUodG90YWwsIG9wZW5lZCkpOwoJCWlmIChoaXQgIT0gbWVtby5lbmQoKSkgcmV0dXJuIGhpdC0+c2Vjb25kOwoKCQlhdXRvIHJlc3VsdCA9CgkJCWxpa2VsaWhvb2QodG90YWwgLSAxLCBvcGVuZWQpICogKDEgLSBvcGVuX2NoYW5jZShvcGVuZWQpKSArCgkJCWxpa2VsaWhvb2QodG90YWwgLSAxLCBvcGVuZWQgLSAxKSAqIG9wZW5fY2hhbmNlKG9wZW5lZCAtIDEpOwoKCQltZW1vW21ha2VfdHVwbGUodG90YWwsIG9wZW5lZCldID0gcmVzdWx0OwoJCWFzc2VydCgwIDw9IHJlc3VsdCAmJiByZXN1bHQgPD0gMSk7CgkJcmV0dXJuIHJlc3VsdDsKCX0KCglkb3VibGUgbG9nX2xpa2VsaWhvb2QoaW50IHRvdGFsLCBpbnQgb3BlbmVkKSB7CgkJcmV0dXJuIGxvZyhsaWtlbGlob29kKHRvdGFsLCBvcGVuZWQpKTsKCX0KfTsKCmludCBtYWluKCkKewoJc3RydWN0IGlucHV0X3R5cGUgewoJCWludCBvcGVuZWQ7CgkJaW50IHRvdGFsOwoJfTsKCWlucHV0X3R5cGUgaW5wdXRzW10gPSB7CgkJaW5wdXRfdHlwZXsgMjUsIDUzIH0sCgkJaW5wdXRfdHlwZXsgMjgsIDcyIH0sCgkJaW5wdXRfdHlwZXsgMjUsIDcyIH0sCgkJaW5wdXRfdHlwZXsgMjksIDcyIH0sCgkJaW5wdXRfdHlwZXsgMjcsIDY5IH0KCX07CgoKCWludCBOID0gNTA7CgoJZm9yIChpbnQgaSA9IC1OOyBpIDw9IE47ICsraSkgewoJCWRvdWJsZSB3aWdodCA9IHBvdygxMCwgKGRvdWJsZSlpIC8gTik7CgkJSHlwb3RoZXNpcyBoeXBvdGhlc2lzKHdpZ2h0KTsKCQlkb3VibGUgc3VtID0gMDsKCQlmb3IgKGF1dG8gaW5wdXQgOiBpbnB1dHMpIHsKCQkJc3VtICs9IGh5cG90aGVzaXMubG9nX2xpa2VsaWhvb2QoaW5wdXQudG90YWwsIGlucHV0Lm9wZW5lZCk7CgkJfQoKCQlpbnQgYmFybGVuZ3RoID0gKGludClyb3VuZCgyMCArIHN1bSk7CgkJaWYgKGJhcmxlbmd0aCA8IDApIGJhcmxlbmd0aCA9IDA7CgkJaWYgKGJhcmxlbmd0aCA+IDIwKSBiYXJsZW5ndGggPSAyMDsKCgkJc3RyaW5nIGJhciA9IHN0cmluZyhiYXJsZW5ndGgsICcgJykgKyBzdHJpbmcoMjAgLSBiYXJsZW5ndGgsICcqJykgKyAifCI7CgoJCWNvdXQgPDwgc2V0dygxMikgPDwgd2lnaHQgPDwgIlx0IiA8PCBzZXR3KDEyKSA8PCBzdW0gPDwgIlx0IyAiIDw8IGJhciA8PCBlbmRsOwoJfQoKfQoK