概要
競技統計学みたいなのの試作（8問）しました。
競技プログラミング→拡張のモチベーション
競技プログラミングは、広範なアルゴリズムの勉強方法として有用であるだけでなく、生産的で現実的な問題を解くことができます。参加者には、大学の先生をはじめとして、信じられないほど能力の高い人が大勢います。そのような、現実世界で生きているだけでは会うことすらできない人たちと、Twitterやコンテスト基盤を介して活発的にコミュニケーションする文化もあります。また、レーティングによって能力が定量化されるため、人財の発掘による経済的メリットも望めます。
これを統計・ロボット・その他に拡張することは世界にとって有益です。今回、古典的な頻度主義的な統計学の問題を数問作成・模範解答を作って、統計解析がジャッジできそうかどうかを議論できるといいなと思っています。
競技学習基盤の拡張
競技学習基盤はどこまで拡張できるでしょうか。アルゴリズムでは、離散アルゴリズム(Topcoder SRMなど)・最適化アルゴリズム(Topcoder Marathon)・予測アルゴリズム(Kaggle)・セキュリティ(SECCON)などがすでにあります。ロボット系では、DARPA Robotics Challengeがシミュレータでの予選を通して、勝ち進んだ人が実際のフィールドで戦います（ただし参入障壁が非常に高く、極めて限定的）。
このような競技学習基盤は、学習だけではなくKaggleのように生産的なケースも多いです。例えば、原発ロボットのコンテストを開催すれば、（現状個人では極めて難しい）ロボット開発スキルを学ぶことができ、それによって社会貢献がなされ、レーティングにより就活にも役立ちます！すごい！
一般的な基盤の制作は困難だと思いますが、とにかく拡張しようという気概を共有したいと思いました。最近勉強しているデータ解析・統計・ベイズ推定・ニューラルネットワークによる最適化などについて、試していきたいと思います。その中でも、実装が軽い頻度主義的統計検定について、「競技プログラミング的問題の試作」を行ってみました。
個人的には頻度主義的検定は好きじゃないので、ベイジアンの人たちは許してください。
あんまり1秒間に100万回統計検定したいみたいなことはないと思うけど、MCMCみたいな話になってくると速度も大事になってくるのかなと思います。
リンク
以下のコードは僕のgithubに置きました。確率密度ライブラリはここのものを使いました。
問題
A - イカサマサイコロ
問題文
問題
n面サイコロを何回か振ったところ、出目iが出た回数がa_i回だった。サイコロは正しく作られているかを、5%有意で検定せよ。
入力
n
a_1 ... a_n
出力
YES or NO
サンプル1
入力
6
103 95 102 97 108 95
出力
NO
解説
正しく作られていると仮定して、この分布が出る確率が5%以下かを両側カイ二乗検定。
解答
#include <bits/stdc++.h>
#include <sys/time.h>
using namespace std;
#define rep(i,n) for(long long i = 0; i < (long long)(n); i++)
#define repi(i,a,b) for(long long i = (long long)(a); i < (long long)(b); i++)
#define pb push_back
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define mt make_tuple
#define mp make_pair
template<class T1, class T2> bool chmin(T1 &a, T2 b) { return b < a && (a = b, true); }
template<class T1, class T2> bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
#define exists find_if
#define forall all_of
using ll = long long; using vll = vector<ll>; using vvll = vector<vll>; using P = pair<ll, ll>;
using ld = double; using vld = vector<ld>;
using vi = vector<int>; using vvi = vector<vi>; vll conv(vi& v) { vll r(v.size()); rep(i, v.size()) r[i] = v[i]; return r; }
using Pos = complex<double>;
template <typename T, typename U> ostream &operator<<(ostream &o, const pair<T, U> &v) { o << "(" << v.first << ", " << v.second << ")"; return o; }
template<size_t...> struct seq{}; template<size_t N, size_t... Is> struct gen_seq : gen_seq<N-1, N-1, Is...>{}; template<size_t... Is> struct gen_seq<0, Is...> : seq<Is...>{};
template<class Ch, class Tr, class Tuple, size_t... Is>
void print_tuple(basic_ostream<Ch,Tr>& os, Tuple const& t, seq<Is...>){ using s = int[]; (void)s{0, (void(os << (Is == 0? "" : ", ") << get<Is>(t)), 0)...}; }
template<class Ch, class Tr, class... Args>
auto operator<<(basic_ostream<Ch, Tr>& os, tuple<Args...> const& t) -> basic_ostream<Ch, Tr>& { os << "("; print_tuple(os, t, gen_seq<sizeof...(Args)>()); return os << ")"; }
ostream &operator<<(ostream &o, const vvll &v) { rep(i, v.size()) { rep(j, v[i].size()) o << v[i][j] << " "; o << endl; } return o; }
template <typename T> ostream &operator<<(ostream &o, const vector<T> &v) { o << '['; rep(i, v.size()) o << v[i] << (i != v.size()-1 ? ", " : ""); o << "]"; return o; }
template <typename T> ostream &operator<<(ostream &o, const set<T> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U> ostream &operator<<(ostream &o, const map<T, U> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U, typename V> ostream &operator<<(ostream &o, const unordered_map<T, U, V> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it; o << "]"; return o; }
vector<int> range(const int x, const int y) { vector<int> v(y - x + 1); iota(v.begin(), v.end(), x); return v; }
template <typename T> istream& operator>>(istream& i, vector<T>& o) { rep(j, o.size()) i >> o[j]; return i;}
string bits_to_string(ll input, ll n=64) { string s; rep(i, n) s += '0' + !!(input & (1ll << i)); reverse(all(s)); return s; }
template <typename T> unordered_map<T, ll> counter(vector<T> vec){unordered_map<T, ll> ret; for (auto&& x : vec) ret[x]++; return ret;};
string substr(string s, P x) {return s.substr(x.fi, x.se - x.fi); }
struct ci : public iterator<forward_iterator_tag, ll> { ll n; ci(const ll n) : n(n) { } bool operator==(const ci& x) { return n == x.n; } bool operator!=(const ci& x) { return !(*this == x); } ci &operator++() { n++; return *this; } ll operator*() const { return n; } };
size_t random_seed; namespace std { using argument_type = P; template<> struct hash<argument_type> { size_t operator()(argument_type const& x) const { size_t seed = random_seed; seed ^= hash<ll>{}(x.fi); seed ^= (hash<ll>{}(x.se) << 1); return seed; } }; }; // hash for various class
namespace myhash{ const int Bsizes[]={3,9,13,17,21,25,29,33,37,41,45,49,53,57,61,65,69,73,77,81}; const int xor_nums[]={0x100007d1,0x5ff049c9,0x14560859,0x07087fef,0x3e277d49,0x4dba1f17,0x709c5988,0x05904258,0x1aa71872,0x238819b3,0x7b002bb7,0x1cf91302,0x0012290a,0x1083576b,0x76473e49,0x3d86295b,0x20536814,0x08634f4d,0x115405e8,0x0e6359f2}; const int hash_key=xor_nums[rand()%20]; const int mod_key=xor_nums[rand()%20]; template <typename T> struct myhash{ std::size_t operator()(const T& val) const { return (hash<T>{}(val)%mod_key)^hash_key; } }; };
template <typename T> class uset:public std::unordered_set<T,myhash::myhash<T>> { using SET=std::unordered_set<T,myhash::myhash<T>>; public: uset():SET(){SET::rehash(myhash::Bsizes[rand()%20]);} };
template <typename T,typename U> class umap:public std::unordered_map<T,U,myhash::myhash<T>> { public: using MAP=std::unordered_map<T,U,myhash::myhash<T>>; umap():MAP(){MAP::rehash(myhash::Bsizes[rand()%20]);} };
struct timeval start; double sec() { struct timeval tv; gettimeofday(&tv, NULL); return (tv.tv_sec - start.tv_sec) + (tv.tv_usec - start.tv_usec) * 1e-6; }
struct init_{init_(){ gettimeofday(&start, NULL); ios::sync_with_stdio(false); cin.tie(0); srand((unsigned int)time(NULL)); random_seed = RAND_MAX / 2 + rand() / 2; }} init__;
static const double EPS = 1e-14;
static const long long INF = 1e18;
static const long long mo = 1e9+7;
#define ldout fixed << setprecision(40)
// 正規分布の累積確率
// A.m.マーリの方法
double qnorm(double u)
{
static double a[9] = { 1.24818987e-4, -1.075204047e-3, 5.198775019e-3,
-0.019198292004, 0.059054035642,-0.151968751364,
0.319152932694,-0.5319230073, 0.797884560593};
static double b[15] = { -4.5255659e-5, 1.5252929e-4, -1.9538132e-5,
-6.76904986e-4, 1.390604284e-3,-7.9462082e-4,
-2.034254874e-3, 6.549791214e-3,-0.010557625006,
0.011630447319,-9.279453341e-3, 5.353579108e-3,
-2.141268741e-3, 5.35310549e-4, 0.999936657524};
double w, y, z;
int i;
if(u == 0.) return 0.5;
y = u / 2.;
if(y < -6.) return 0.;
if(y > 6.) return 1.;
if(y < 0.) y = - y;
if(y < 1.) {
w = y * y;
z = a[0];
for(i = 1; i < 9; i++) z = z * w + a[i];
z *= (y * 2.);
} else {
y -= 2.;
z = b[0];
for(i = 1; i < 15; i++) z = z * y + b[i];
}
if(u < 0.) return (1. - z) / 2.;
return (1. + z) / 2.;
}
// 正規分布のパーセント点
// 戸田の近似式
double pnorm(double qn)
{
static double b[11] = { 1.570796288, 0.03706987906, -0.8364353589e-3,
-0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
-0.104527497e-5, 0.8360937017e-7,-0.3231081277e-8,
0.3657763036e-10,0.6936233982e-12};
double w1, w3;
int i;
if(qn < 0. || 1. < qn) {
fprintf(stderr, "Error : qn <= 0 or qn >= 1 in pnorm()!\n");
return 0.;
}
if(qn == 0.5) return 0.;
w1 = qn;
if(qn > 0.5) w1 = 1. - w1;
w3 = -log(4. * w1 * (1. - w1));
w1 = b[0];
for(i = 1; i < 11; i++) w1 += (b[i] * pow(w3, (double)i));
if(qn > 0.5) return sqrt(w1 * w3);
return -sqrt(w1 * w3);
}// 自由度nのカイ2乗の上側確率
double qchi(double x2, int n)
{
double w, pw, x, qc;
int i, i1;
if(n < 1) {
fprintf(stderr,"Error : 自由度 < 1 in qchi()!\n");
return 0.;
}
if(x2 <= 0.) return 1.;
if(x2 > 400.) return 0.;
if(n > 10) {
w = 2. / (9. * (double)n);
pw = pow(x2 / (double)n, 1. / 3.);
return 1. - qnorm((pw - 1. + w) / sqrt(w));
}
w = exp(-x2 / 2.);
if(n == 2) return w;
x = sqrt(x2);
if(n == 1) return 2. * (1. - qnorm(x));
if((n % 2) == 0) {
i1 = 2;
qc = w;
} else {
i1 = 1;
qc = 2. * (1. - qnorm(x));
w *= (0.797884560750774 / x);
}
for(i = i1; i <= n - 2; i += 2) {
w *= (x2 / (double)i);
qc += w;
}
return qc;
}
// 自由度nのカイ2乗のパーセント点
double pchi(double qc, int n)
{
double c1, c2, gam, x, w, wx;
int i, j;
if(qc <= 0. || qc >= 1. || n < 1) {
fprintf(stderr,"Error : Illigal parameter in pchi()!\n");
return 0.;
}
if(n == 1) {
w = pnorm(qc / 2.);
return w * w;
}
if(n == 2) return -2. * log(qc);
x = -pnorm(qc);
if(n > 10) {
w = x * x;
wx = sqrt(2. * (double)n);
c1 = (w - 7.) * x / 9. / wx;
c2 = ((3. * w + 7.) * w - 16.) * 2. / 405. / (double)n;
wx = (double)n + wx * x + 2. * (w - 1.) / 3. + c1 - c2;
if(wx < 0.) return 0.;
return wx;
}
w = 2. / 9. / (double)n;
w = 1. - w + x * sqrt(w);
wx = (double)n * w * w * w;
if((n % 2) == 0) gam = 1.;
else gam = 1.772453850905516;
j = (n + 1) / 2 - 1;
w = (double)n / 2.;
for(i = 1; i <= j; i++) gam *= (w - (double)i);
x = wx / 2.;
c1 = pow(x, w - 1.);
c2 = exp(-x) / 2.;
return wx + (qchi(wx, n) - qc) * gam / c1 / c2;
}
double mean(vector<double> a) {
double ret = 0;
rep(i, a.size()) {
ret += a[i];
}
return ret / a.size();
}
double stdev(vector<double> a) {
double ret = 0;
double m = mean(a);
rep(i, a.size()) {
ret += (a[i] - m) * (a[i] - m);
}
return sqrt(ret / (a.size() - 1));
}
int main(void) {
ll n; cin >> n;
vector<double> a(n); cin >> a;
double m = accumulate(all(a), 0ll);
rep(i, n) a[i] -= (m / n);
rep(i, n) a[i] *= a[i];
rep(i, n) a[i] /= (m / n);
double chisq = accumulate(all(a), 0.0);
if (chisq > pchi(0.05, n-1)) {
cout << "YES" << endl;
} else {
cout << "NO" << endl;
}
return 0;
}
B - はむこダイエット
問題文
問題
はむこダイエット法を被験者n人に行った。すると、行う前後で体重が、a_iからb_iに変化した。はむこダイエット法は、5%有意で効果があるか？
入力
n
a_1 ... a_n
b_1 ... b_n
出力
YES or NO
サンプル1
入力
8
54.3 51.9 55.2 55.3 53.8 52.1 53.6 50.2
52.6 48.7 56.7 52.5 52.1 53.3 51.4 49.2
出力
YES
解説
対応ある片側t検定でホイ
解答
#include <bits/stdc++.h>
#include <sys/time.h>
using namespace std;
#define rep(i,n) for(long long i = 0; i < (long long)(n); i++)
#define repi(i,a,b) for(long long i = (long long)(a); i < (long long)(b); i++)
#define pb push_back
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define mt make_tuple
#define mp make_pair
template<class T1, class T2> bool chmin(T1 &a, T2 b) { return b < a && (a = b, true); }
template<class T1, class T2> bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
#define exists find_if
#define forall all_of
using ll = long long; using vll = vector<ll>; using vvll = vector<vll>; using P = pair<ll, ll>;
using ld = double; using vld = vector<ld>;
using vi = vector<int>; using vvi = vector<vi>; vll conv(vi& v) { vll r(v.size()); rep(i, v.size()) r[i] = v[i]; return r; }
using Pos = complex<double>;
template <typename T, typename U> ostream &operator<<(ostream &o, const pair<T, U> &v) { o << "(" << v.first << ", " << v.second << ")"; return o; }
template<size_t...> struct seq{}; template<size_t N, size_t... Is> struct gen_seq : gen_seq<N-1, N-1, Is...>{}; template<size_t... Is> struct gen_seq<0, Is...> : seq<Is...>{};
template<class Ch, class Tr, class Tuple, size_t... Is>
void print_tuple(basic_ostream<Ch,Tr>& os, Tuple const& t, seq<Is...>){ using s = int[]; (void)s{0, (void(os << (Is == 0? "" : ", ") << get<Is>(t)), 0)...}; }
template<class Ch, class Tr, class... Args>
auto operator<<(basic_ostream<Ch, Tr>& os, tuple<Args...> const& t) -> basic_ostream<Ch, Tr>& { os << "("; print_tuple(os, t, gen_seq<sizeof...(Args)>()); return os << ")"; }
ostream &operator<<(ostream &o, const vvll &v) { rep(i, v.size()) { rep(j, v[i].size()) o << v[i][j] << " "; o << endl; } return o; }
template <typename T> ostream &operator<<(ostream &o, const vector<T> &v) { o << '['; rep(i, v.size()) o << v[i] << (i != v.size()-1 ? ", " : ""); o << "]"; return o; }
template <typename T> ostream &operator<<(ostream &o, const set<T> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U> ostream &operator<<(ostream &o, const map<T, U> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U, typename V> ostream &operator<<(ostream &o, const unordered_map<T, U, V> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it; o << "]"; return o; }
vector<int> range(const int x, const int y) { vector<int> v(y - x + 1); iota(v.begin(), v.end(), x); return v; }
template <typename T> istream& operator>>(istream& i, vector<T>& o) { rep(j, o.size()) i >> o[j]; return i;}
string bits_to_string(ll input, ll n=64) { string s; rep(i, n) s += '0' + !!(input & (1ll << i)); reverse(all(s)); return s; }
template <typename T> unordered_map<T, ll> counter(vector<T> vec){unordered_map<T, ll> ret; for (auto&& x : vec) ret[x]++; return ret;};
string substr(string s, P x) {return s.substr(x.fi, x.se - x.fi); }
struct ci : public iterator<forward_iterator_tag, ll> { ll n; ci(const ll n) : n(n) { } bool operator==(const ci& x) { return n == x.n; } bool operator!=(const ci& x) { return !(*this == x); } ci &operator++() { n++; return *this; } ll operator*() const { return n; } };
size_t random_seed; namespace std { using argument_type = P; template<> struct hash<argument_type> { size_t operator()(argument_type const& x) const { size_t seed = random_seed; seed ^= hash<ll>{}(x.fi); seed ^= (hash<ll>{}(x.se) << 1); return seed; } }; }; // hash for various class
namespace myhash{ const int Bsizes[]={3,9,13,17,21,25,29,33,37,41,45,49,53,57,61,65,69,73,77,81}; const int xor_nums[]={0x100007d1,0x5ff049c9,0x14560859,0x07087fef,0x3e277d49,0x4dba1f17,0x709c5988,0x05904258,0x1aa71872,0x238819b3,0x7b002bb7,0x1cf91302,0x0012290a,0x1083576b,0x76473e49,0x3d86295b,0x20536814,0x08634f4d,0x115405e8,0x0e6359f2}; const int hash_key=xor_nums[rand()%20]; const int mod_key=xor_nums[rand()%20]; template <typename T> struct myhash{ std::size_t operator()(const T& val) const { return (hash<T>{}(val)%mod_key)^hash_key; } }; };
template <typename T> class uset:public std::unordered_set<T,myhash::myhash<T>> { using SET=std::unordered_set<T,myhash::myhash<T>>; public: uset():SET(){SET::rehash(myhash::Bsizes[rand()%20]);} };
template <typename T,typename U> class umap:public std::unordered_map<T,U,myhash::myhash<T>> { public: using MAP=std::unordered_map<T,U,myhash::myhash<T>>; umap():MAP(){MAP::rehash(myhash::Bsizes[rand()%20]);} };
struct timeval start; double sec() { struct timeval tv; gettimeofday(&tv, NULL); return (tv.tv_sec - start.tv_sec) + (tv.tv_usec - start.tv_usec) * 1e-6; }
struct init_{init_(){ gettimeofday(&start, NULL); ios::sync_with_stdio(false); cin.tie(0); srand((unsigned int)time(NULL)); random_seed = RAND_MAX / 2 + rand() / 2; }} init__;
static const double EPS = 1e-14;
static const long long INF = 1e18;
static const long long mo = 1e9+7;
#define ldout fixed << setprecision(40)
// 正規分布のパーセント点
// 戸田の近似式
double pnorm(double qn)
{
static double b[11] = { 1.570796288, 0.03706987906, -0.8364353589e-3,
-0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
-0.104527497e-5, 0.8360937017e-7,-0.3231081277e-8,
0.3657763036e-10,0.6936233982e-12};
double w1, w3;
int i;
if(qn < 0. || 1. < qn) {
fprintf(stderr, "Error : qn <= 0 or qn >= 1 in pnorm()!\n");
return 0.;
}
if(qn == 0.5) return 0.;
w1 = qn;
if(qn > 0.5) w1 = 1. - w1;
w3 = -log(4. * w1 * (1. - w1));
w1 = b[0];
for(i = 1; i < 11; i++) w1 += (b[i] * pow(w3, (double)i));
if(qn > 0.5) return sqrt(w1 * w3);
return -sqrt(w1 * w3);
}
// 自由度nのt分布の両側確率
void qtsub(double *q, int n, double w2, double w3, double t2)
{
int i, j;
double w;
j = (n - 2) / 2;
for(i = 1; i <= j; i++) {
w = 2. * (double)i - w3;
*q -= (w2 *= (t2 * w / (1. + w)));
}
}
double qt(double t, int n)
{
double t1, t2, w1, w2, wq;
if(n < 1)
{
fprintf(stderr,"Error : n < 1 in qt()!\n");
return 0.;
}
w1 = 0.636619772284456;
if(t < 0.) t = - t;
t1 = t / sqrt((double)n);
t2 = 1. / (1. + t1 * t1);
if((n % 2) != 0) {
wq = 1. - w1 * atan(t1);
if(n != 1) {
wq -= (w2 = w1 * t1 * t2);
if (n != 3) qtsub(&wq, n, w2, 0., t2);
}
if(wq > 0.) return wq;
return 0.;
}
wq = 1. - (w2 = t1 * sqrt(t2));
if(n != 2) qtsub(&wq, n, w2, 1., t2);
if(wq > 0.) return wq;
return 0.;
}
// 自由度nのt分布のパーセント点
double ptsub(double q, int n)
{
double eps, qe, s, w;
if(n == 1 && 0.001 <= q && q < 0.01) eps = 1.e-4;
else if (n == 2 && q < 0.0001) eps = 1.e-4;
else if (n == 1 && q < 0.001) eps = 1.e-2;
else eps = 1.e-5;
s = 10000.;
w = 0.;
for(;;) {
w += s;
if(s <= eps) return w;
if((qe = qt(w,n) - q) == 0.) return w;
if(qe < 0.) {
w -= s;
s /= 10.;
}
}
}
double pt(double q, int n)
{
double f, f1, f2, f3, f4, f5, u, u2, w, w0, w1, w2, w3, w4;
if(q < 1.e-5 || q > 1. || n < 1) {
fprintf(stderr,"Error : Illigal parameter in pt()!\n");
return 0.;
}
if(n <= 5) return ptsub(q, n);
if(q <= 5.e-3 && n <= 13) return ptsub(q, n);
f1 = 4. * (f = (double)n);
f5 = (f4 = (f3 = (f2 = f * f) * f) * f) * f;
f2 *= 96.;
f3 *= 384.;
f4 *= 92160.;
f5 *= 368640.;
u = pnorm(1. - q / 2.);
w0 = (u2 = u * u) * u;
w1 = w0 * u2;
w2 = w1 * u2;
w3 = w2 * u2;
w4 = w3 * u2;
w = ((w0 + u) / f1);
w += ((5. * w1 + 16. * w0 + 3. * u) / f2);
w += ((3. * w2 + 19. * w1 + 17. * w0 - 15. * u) / f3);
w += ((79. * w3 + 776. * w2 + 1482. * w1 - 1920. * w0 - 945. * u) / f4);
w += ((27. * w4 + 339. * w3 + 930. * w2 - 1782. * w1 - 765. * w0
+ 17955. * u) / f5);
return u + w;
}
double mean(vector<double> a) {
double ret = 0;
rep(i, a.size()) {
ret += a[i];
}
return ret / a.size();
}
double stdev(vector<double> a) {
double ret = 0;
double m = mean(a);
rep(i, a.size()) {
ret += (a[i] - m) * (a[i] - m);
}
return sqrt(ret / (a.size() - 1));
}
int main(void) {
ll n; cin >> n;
vector<double> a(n), b(n);
cin >> a >> b;
// 後から見た前の体重が有意に大きいか？を調べたい
rep(i, n) a[i] -= b[i];
// 対応あるt検定、片側
double t = mean(a) / stdev(a) * sqrt(n);
cout << t << endl;
if (t > pt(0.05*2, n-1)) {
cout << "YES" << endl;
} else {
cout << "NO" << endl;
}
return 0;
}
C - 男女差
問題文
問題
男子n人女子m人の期末テストを受けた。男子iのテスト結果がa_i, 女子jのテスト結果がb_jである。この時、男女間で平均値に5%有意な差はあるか？
入力
n m
a_1 ... a_n
b_1 ... a_m
出力
YES or NO
サンプル1
入力
15 10
65 69 69 66 66 51 63 46 63 73 64 61 72 59 74
54 68 68 55 55 55 56 55 52 56
出力
YES
サンプル2
入力
10 11
15.3 14.9 14.5 14.4 14.0 13.9 14.1 14.7 15.3 14.6
13.9 14.2 14.1 14.3 14.1 13.7 14.7 13.9 14.1 13.8 14.3
出力
YES
解説
対応のないt検定、等分散仮定なし。
なので、F検定で等分散仮定をチェックしてから、t検定を選ぶ。
解答
#include <bits/stdc++.h>
#include <sys/time.h>
using namespace std;
#define rep(i,n) for(long long i = 0; i < (long long)(n); i++)
#define repi(i,a,b) for(long long i = (long long)(a); i < (long long)(b); i++)
#define pb push_back
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define mt make_tuple
#define mp make_pair
template<class T1, class T2> bool chmin(T1 &a, T2 b) { return b < a && (a = b, true); }
template<class T1, class T2> bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
#define exists find_if
#define forall all_of
using ll = long long; using vll = vector<ll>; using vvll = vector<vll>; using P = pair<ll, ll>;
using ld = double; using vld = vector<ld>;
using vi = vector<int>; using vvi = vector<vi>; vll conv(vi& v) { vll r(v.size()); rep(i, v.size()) r[i] = v[i]; return r; }
using Pos = complex<double>;
template <typename T, typename U> ostream &operator<<(ostream &o, const pair<T, U> &v) { o << "(" << v.first << ", " << v.second << ")"; return o; }
template<size_t...> struct seq{}; template<size_t N, size_t... Is> struct gen_seq : gen_seq<N-1, N-1, Is...>{}; template<size_t... Is> struct gen_seq<0, Is...> : seq<Is...>{};
template<class Ch, class Tr, class Tuple, size_t... Is>
void print_tuple(basic_ostream<Ch,Tr>& os, Tuple const& t, seq<Is...>){ using s = int[]; (void)s{0, (void(os << (Is == 0? "" : ", ") << get<Is>(t)), 0)...}; }
template<class Ch, class Tr, class... Args>
auto operator<<(basic_ostream<Ch, Tr>& os, tuple<Args...> const& t) -> basic_ostream<Ch, Tr>& { os << "("; print_tuple(os, t, gen_seq<sizeof...(Args)>()); return os << ")"; }
ostream &operator<<(ostream &o, const vvll &v) { rep(i, v.size()) { rep(j, v[i].size()) o << v[i][j] << " "; o << endl; } return o; }
template <typename T> ostream &operator<<(ostream &o, const vector<T> &v) { o << '['; rep(i, v.size()) o << v[i] << (i != v.size()-1 ? ", " : ""); o << "]"; return o; }
template <typename T> ostream &operator<<(ostream &o, const set<T> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U> ostream &operator<<(ostream &o, const map<T, U> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U, typename V> ostream &operator<<(ostream &o, const unordered_map<T, U, V> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it; o << "]"; return o; }
vector<int> range(const int x, const int y) { vector<int> v(y - x + 1); iota(v.begin(), v.end(), x); return v; }
template <typename T> istream& operator>>(istream& i, vector<T>& o) { rep(j, o.size()) i >> o[j]; return i;}
string bits_to_string(ll input, ll n=64) { string s; rep(i, n) s += '0' + !!(input & (1ll << i)); reverse(all(s)); return s; }
template <typename T> unordered_map<T, ll> counter(vector<T> vec){unordered_map<T, ll> ret; for (auto&& x : vec) ret[x]++; return ret;};
string substr(string s, P x) {return s.substr(x.fi, x.se - x.fi); }
struct ci : public iterator<forward_iterator_tag, ll> { ll n; ci(const ll n) : n(n) { } bool operator==(const ci& x) { return n == x.n; } bool operator!=(const ci& x) { return !(*this == x); } ci &operator++() { n++; return *this; } ll operator*() const { return n; } };
size_t random_seed; namespace std { using argument_type = P; template<> struct hash<argument_type> { size_t operator()(argument_type const& x) const { size_t seed = random_seed; seed ^= hash<ll>{}(x.fi); seed ^= (hash<ll>{}(x.se) << 1); return seed; } }; }; // hash for various class
namespace myhash{ const int Bsizes[]={3,9,13,17,21,25,29,33,37,41,45,49,53,57,61,65,69,73,77,81}; const int xor_nums[]={0x100007d1,0x5ff049c9,0x14560859,0x07087fef,0x3e277d49,0x4dba1f17,0x709c5988,0x05904258,0x1aa71872,0x238819b3,0x7b002bb7,0x1cf91302,0x0012290a,0x1083576b,0x76473e49,0x3d86295b,0x20536814,0x08634f4d,0x115405e8,0x0e6359f2}; const int hash_key=xor_nums[rand()%20]; const int mod_key=xor_nums[rand()%20]; template <typename T> struct myhash{ std::size_t operator()(const T& val) const { return (hash<T>{}(val)%mod_key)^hash_key; } }; };
template <typename T> class uset:public std::unordered_set<T,myhash::myhash<T>> { using SET=std::unordered_set<T,myhash::myhash<T>>; public: uset():SET(){SET::rehash(myhash::Bsizes[rand()%20]);} };
template <typename T,typename U> class umap:public std::unordered_map<T,U,myhash::myhash<T>> { public: using MAP=std::unordered_map<T,U,myhash::myhash<T>>; umap():MAP(){MAP::rehash(myhash::Bsizes[rand()%20]);} };
struct timeval start; double sec() { struct timeval tv; gettimeofday(&tv, NULL); return (tv.tv_sec - start.tv_sec) + (tv.tv_usec - start.tv_usec) * 1e-6; }
struct init_{init_(){ gettimeofday(&start, NULL); ios::sync_with_stdio(false); cin.tie(0); srand((unsigned int)time(NULL)); random_seed = RAND_MAX / 2 + rand() / 2; }} init__;
static const double EPS = 1e-14;
static const long long INF = 1e18;
static const long long mo = 1e9+7;
#define ldout fixed << setprecision(40)
// 正規分布のパーセント点
// 戸田の近似式
double pnorm(double qn)
{
static double b[11] = { 1.570796288, 0.03706987906, -0.8364353589e-3,
-0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
-0.104527497e-5, 0.8360937017e-7,-0.3231081277e-8,
0.3657763036e-10,0.6936233982e-12};
double w1, w3;
int i;
if(qn < 0. || 1. < qn) {
fprintf(stderr, "Error : qn <= 0 or qn >= 1 in pnorm()!\n");
return 0.;
}
if(qn == 0.5) return 0.;
w1 = qn;
if(qn > 0.5) w1 = 1. - w1;
w3 = -log(4. * w1 * (1. - w1));
w1 = b[0];
for(i = 1; i < 11; i++) w1 += (b[i] * pow(w3, (double)i));
if(qn > 0.5) return sqrt(w1 * w3);
return -sqrt(w1 * w3);
}
// 自由度nのt分布の両側確率
void qtsub(double *q, int n, double w2, double w3, double t2)
{
int i, j;
double w;
j = (n - 2) / 2;
for(i = 1; i <= j; i++) {
w = 2. * (double)i - w3;
*q -= (w2 *= (t2 * w / (1. + w)));
}
}
double qt(double t, int n)
{
double t1, t2, w1, w2, wq;
if(n < 1)
{
fprintf(stderr,"Error : n < 1 in qt()!\n");
return 0.;
}
w1 = 0.636619772284456;
if(t < 0.) t = - t;
t1 = t / sqrt((double)n);
t2 = 1. / (1. + t1 * t1);
if((n % 2) != 0) {
wq = 1. - w1 * atan(t1);
if(n != 1) {
wq -= (w2 = w1 * t1 * t2);
if (n != 3) qtsub(&wq, n, w2, 0., t2);
}
if(wq > 0.) return wq;
return 0.;
}
wq = 1. - (w2 = t1 * sqrt(t2));
if(n != 2) qtsub(&wq, n, w2, 1., t2);
if(wq > 0.) return wq;
return 0.;
}
// 自由度nのt分布のパーセント点
double ptsub(double q, int n)
{
double eps, qe, s, w;
if(n == 1 && 0.001 <= q && q < 0.01) eps = 1.e-4;
else if (n == 2 && q < 0.0001) eps = 1.e-4;
else if (n == 1 && q < 0.001) eps = 1.e-2;
else eps = 1.e-5;
s = 10000.;
w = 0.;
for(;;) {
w += s;
if(s <= eps) return w;
if((qe = qt(w,n) - q) == 0.) return w;
if(qe < 0.) {
w -= s;
s /= 10.;
}
}
}
double pt(double q, int n)
{
double f, f1, f2, f3, f4, f5, u, u2, w, w0, w1, w2, w3, w4;
if(q < 1.e-5 || q > 1. || n < 1) {
fprintf(stderr,"Error : Illigal parameter in pt()!\n");
return 0.;
}
if(n <= 5) return ptsub(q, n);
if(q <= 5.e-3 && n <= 13) return ptsub(q, n);
f1 = 4. * (f = (double)n);
f5 = (f4 = (f3 = (f2 = f * f) * f) * f) * f;
f2 *= 96.;
f3 *= 384.;
f4 *= 92160.;
f5 *= 368640.;
u = pnorm(1. - q / 2.);
w0 = (u2 = u * u) * u;
w1 = w0 * u2;
w2 = w1 * u2;
w3 = w2 * u2;
w4 = w3 * u2;
w = ((w0 + u) / f1);
w += ((5. * w1 + 16. * w0 + 3. * u) / f2);
w += ((3. * w2 + 19. * w1 + 17. * w0 - 15. * u) / f3);
w += ((79. * w3 + 776. * w2 + 1482. * w1 - 1920. * w0 - 945. * u) / f4);
w += ((27. * w4 + 339. * w3 + 930. * w2 - 1782. * w1 - 765. * w0
+ 17955. * u) / f5);
return u + w;
}
// 自由度n1, n2のF分布の上側確率
double qf(double f, int n1, int n2)
{
double d, p, w, y, z;
int j, k, m, mn, n;
if(f < 0. || n1 < 1 || n2 < 1) {
fprintf(stderr,"Error : Illegal parameter in qf()!\n");
return 0.;
}
w = f * (double)n1 / (double)n2;
z = 1. / (1. + w);
m = 2 * (n1 / 2) - n1 + 2;
n = 2 * (n2 / 2) - n2 + 2;
if((mn = 4 - 2 * (n1 % 2) - (n2 % 2)) == 1) {
y = 0.318309886142228;
p = sqrt(w);
d = y * z / p;
p = 2. * y * atan(p);
} else if(mn == 2) {
p = sqrt(z * w);
d = p * z / w / 2.;
} else if(mn == 3) {
p = sqrt(z);
d = z * p / 2.;
p = 1. - p;
} else {
d = z * z;
p = z * w;
}
if(n1 <= 2 && n2 <= 2) return 1. - p;
y = 2. * w / z;
if((k = n + 2) <= n2) {
for( ; k <= n2; k += 2) {
d *= ((1. + (double)m / (double)(k - 2)) * z);
if(m != 1) p = (p + w) * z;
else p += (d * y / (double)(k - 1));
}
}
k = m + 2;
if(k > n1) return 1. - p;
y = z * w;
z = 2. / z;
j = n2 - 2;
for( ; k <= n1; k += 2) {
d *= (y * (double)(k + j) / (double)(k - 2));
p -= (z * d / (double)(k + j));
}
return 1. - p;
}
double pfsub(double x, double y, double z) {
return (sqrt(z) - y) / x / 2.;
}
// 自由度n1, n2のF分布のパーセント点
double pf(double q, int n1, int n2)
{
double a, b, c, d, eps, fw, qe, qn, s, u, u2, w1, w2, w3, w4;
if(q < 0. || q > 1. || n1 < 1 || n2 < 1) {
fprintf(stderr,"Error : Illegal parameter in pf()!\n");
return 0.;
}
if(n1 <= 240 || n2 <= 240) {
eps = 1.e-5;
if(n2 == 1) eps = 1.e-4;
s = 1000.;
fw = 0.;
for(;;) {
fw += s;
if(s <= eps) return fw;
if((qe = qf(fw,n1,n2) - q) == 0.) return fw;
if(qe < 0.) {
fw -= s;
s /= 10.;
}
}
}
eps = 1.e-6;
qn = q;
if(q < 0.5) qn = 1. - q;
u = pnorm(qn);
w1 = 2. / (double)n1 / 9.;
w2 = 2. / (double)n2 / 9.;
w3 = 1. - w1;
w4 = 1. - w2;
u2 = u * u;
a = w4 * w4 - u2 * w2;
b = -2. * w3 * w4;
c = w3 * w3 - u2 * w1;
d = b * b - 4 * a * c;
if(d < 0.)
fw = pfsub(a, b, 0.);
else {
if(fabs(a) > eps)
fw = pfsub(a, b, d);
else {
if(fabs(b) > eps) return -c / b;
fw = pfsub(a, b, 0.);
}
}
return fw * fw * fw;
}
// http://www5.airnet.ne.jp/tomy/cpro/sslib11.htm
double mean(vector<double> a) {
double ret = 0;
rep(i, a.size()) {
ret += a[i];
}
return ret / a.size();
}
double stdev(vector<double> a) {
double ret = 0;
double m = mean(a);
rep(i, a.size()) {
ret += (a[i] - m) * (a[i] - m);
}
return sqrt(ret / (a.size() - 1));
}
int main(void) {
int m, n; cin >> m >> n;
vector<double> a(m), b(n); cin >> a >> b;
// (m, a, sa)のほうが分散が多くswap
double sa = stdev(a), sb = stdev(b);
if (sa < sb) swap(sa, sb), swap(n, m), swap(a, b);
double ua = sa * sa, ub = sb * sb;
double f = ua / ub;
double t = 0, v = 0;
if (f < pf(0.05, m-1, n-1)) {
// 等分散仮定OK
// 等分散のt検定, 両側
double u = ((m-1)*ua+(n-1)*ub) / (n + m - 2);
t = abs(mean(a) - mean(b)) / sqrt(u * (1.0 / m + 1.0 / n));
v = n + m - 2;
} else {
// 等分散仮定not OK
// ウェルチのt検定, 両側
t = abs(mean(a) - mean(b)) / sqrt(ua/m+ub/n);
v = pow(ua/m+ub/n, 2) / (ua*ua/m/m/(m-1)+ub*ub/n/n/(n-1));
}
if (t > pt(0.05, v)){
cout << "YES" << endl;
} else {
cout << "NO" << endl;
}
return 0;
}
D - 正確なビタミンC
問題文
問題
ある錠剤は1粒あたり平均m_target [mg]のビタミンCを含むものと定められている．
錠剤の標本n個を無作為抽出したとき，ビタミンCの含有量は1粒当り平均m[mg]，標本単純標準偏差s [mg]であった．
この標本を取り出した母集団の錠剤は規格を満たすか？有意水準1%で検定せよ。
入力
m_target n m s
出力
検定により、規格を満たしそうならばYES、満たさなそうならばNOを出力せよ。
サンプル1
入力
250 20 249.5 0.5
出力
NO
解説
t分布
解答
#include <bits/stdc++.h>
#include <sys/time.h>
using namespace std;
#define rep(i,n) for(long long i = 0; i < (long long)(n); i++)
#define repi(i,a,b) for(long long i = (long long)(a); i < (long long)(b); i++)
#define pb push_back
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define mt make_tuple
#define mp make_pair
template<class T1, class T2> bool chmin(T1 &a, T2 b) { return b < a && (a = b, true); }
template<class T1, class T2> bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
#define exists find_if
#define forall all_of
using ll = long long; using vll = vector<ll>; using vvll = vector<vll>; using P = pair<ll, ll>;
using ld = double; using vld = vector<ld>;
using vi = vector<int>; using vvi = vector<vi>; vll conv(vi& v) { vll r(v.size()); rep(i, v.size()) r[i] = v[i]; return r; }
using Pos = complex<double>;
template <typename T, typename U> ostream &operator<<(ostream &o, const pair<T, U> &v) { o << "(" << v.first << ", " << v.second << ")"; return o; }
template<size_t...> struct seq{}; template<size_t N, size_t... Is> struct gen_seq : gen_seq<N-1, N-1, Is...>{}; template<size_t... Is> struct gen_seq<0, Is...> : seq<Is...>{};
template<class Ch, class Tr, class Tuple, size_t... Is>
void print_tuple(basic_ostream<Ch,Tr>& os, Tuple const& t, seq<Is...>){ using s = int[]; (void)s{0, (void(os << (Is == 0? "" : ", ") << get<Is>(t)), 0)...}; }
template<class Ch, class Tr, class... Args>
auto operator<<(basic_ostream<Ch, Tr>& os, tuple<Args...> const& t) -> basic_ostream<Ch, Tr>& { os << "("; print_tuple(os, t, gen_seq<sizeof...(Args)>()); return os << ")"; }
ostream &operator<<(ostream &o, const vvll &v) { rep(i, v.size()) { rep(j, v[i].size()) o << v[i][j] << " "; o << endl; } return o; }
template <typename T> ostream &operator<<(ostream &o, const vector<T> &v) { o << '['; rep(i, v.size()) o << v[i] << (i != v.size()-1 ? ", " : ""); o << "]"; return o; }
template <typename T> ostream &operator<<(ostream &o, const set<T> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U> ostream &operator<<(ostream &o, const map<T, U> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U, typename V> ostream &operator<<(ostream &o, const unordered_map<T, U, V> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it; o << "]"; return o; }
vector<int> range(const int x, const int y) { vector<int> v(y - x + 1); iota(v.begin(), v.end(), x); return v; }
template <typename T> istream& operator>>(istream& i, vector<T>& o) { rep(j, o.size()) i >> o[j]; return i;}
string bits_to_string(ll input, ll n=64) { string s; rep(i, n) s += '0' + !!(input & (1ll << i)); reverse(all(s)); return s; }
template <typename T> unordered_map<T, ll> counter(vector<T> vec){unordered_map<T, ll> ret; for (auto&& x : vec) ret[x]++; return ret;};
string substr(string s, P x) {return s.substr(x.fi, x.se - x.fi); }
struct ci : public iterator<forward_iterator_tag, ll> { ll n; ci(const ll n) : n(n) { } bool operator==(const ci& x) { return n == x.n; } bool operator!=(const ci& x) { return !(*this == x); } ci &operator++() { n++; return *this; } ll operator*() const { return n; } };
size_t random_seed; namespace std { using argument_type = P; template<> struct hash<argument_type> { size_t operator()(argument_type const& x) const { size_t seed = random_seed; seed ^= hash<ll>{}(x.fi); seed ^= (hash<ll>{}(x.se) << 1); return seed; } }; }; // hash for various class
namespace myhash{ const int Bsizes[]={3,9,13,17,21,25,29,33,37,41,45,49,53,57,61,65,69,73,77,81}; const int xor_nums[]={0x100007d1,0x5ff049c9,0x14560859,0x07087fef,0x3e277d49,0x4dba1f17,0x709c5988,0x05904258,0x1aa71872,0x238819b3,0x7b002bb7,0x1cf91302,0x0012290a,0x1083576b,0x76473e49,0x3d86295b,0x20536814,0x08634f4d,0x115405e8,0x0e6359f2}; const int hash_key=xor_nums[rand()%20]; const int mod_key=xor_nums[rand()%20]; template <typename T> struct myhash{ std::size_t operator()(const T& val) const { return (hash<T>{}(val)%mod_key)^hash_key; } }; };
template <typename T> class uset:public std::unordered_set<T,myhash::myhash<T>> { using SET=std::unordered_set<T,myhash::myhash<T>>; public: uset():SET(){SET::rehash(myhash::Bsizes[rand()%20]);} };
template <typename T,typename U> class umap:public std::unordered_map<T,U,myhash::myhash<T>> { public: using MAP=std::unordered_map<T,U,myhash::myhash<T>>; umap():MAP(){MAP::rehash(myhash::Bsizes[rand()%20]);} };
struct timeval start; double sec() { struct timeval tv; gettimeofday(&tv, NULL); return (tv.tv_sec - start.tv_sec) + (tv.tv_usec - start.tv_usec) * 1e-6; }
struct init_{init_(){ gettimeofday(&start, NULL); ios::sync_with_stdio(false); cin.tie(0); srand((unsigned int)time(NULL)); random_seed = RAND_MAX / 2 + rand() / 2; }} init__;
static const double EPS = 1e-14;
static const long long INF = 1e18;
static const long long mo = 1e9+7;
#define ldout fixed << setprecision(40)
// 正規分布のパーセント点
// 戸田の近似式
double pnorm(double qn)
{
static double b[11] = { 1.570796288, 0.03706987906, -0.8364353589e-3,
-0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
-0.104527497e-5, 0.8360937017e-7,-0.3231081277e-8,
0.3657763036e-10,0.6936233982e-12};
double w1, w3;
int i;
if(qn < 0. || 1. < qn) {
fprintf(stderr, "Error : qn <= 0 or qn >= 1 in pnorm()!\n");
return 0.;
}
if(qn == 0.5) return 0.;
w1 = qn;
if(qn > 0.5) w1 = 1. - w1;
w3 = -log(4. * w1 * (1. - w1));
w1 = b[0];
for(i = 1; i < 11; i++) w1 += (b[i] * pow(w3, (double)i));
if(qn > 0.5) return sqrt(w1 * w3);
return -sqrt(w1 * w3);
}
// 自由度nのt分布の両側確率
void qtsub(double *q, int n, double w2, double w3, double t2)
{
int i, j;
double w;
j = (n - 2) / 2;
for(i = 1; i <= j; i++) {
w = 2. * (double)i - w3;
*q -= (w2 *= (t2 * w / (1. + w)));
}
}
double qt(double t, int n)
{
double t1, t2, w1, w2, wq;
if(n < 1)
{
fprintf(stderr,"Error : n < 1 in qt()!\n");
return 0.;
}
w1 = 0.636619772284456;
if(t < 0.) t = - t;
t1 = t / sqrt((double)n);
t2 = 1. / (1. + t1 * t1);
if((n % 2) != 0) {
wq = 1. - w1 * atan(t1);
if(n != 1) {
wq -= (w2 = w1 * t1 * t2);
if (n != 3) qtsub(&wq, n, w2, 0., t2);
}
if(wq > 0.) return wq;
return 0.;
}
wq = 1. - (w2 = t1 * sqrt(t2));
if(n != 2) qtsub(&wq, n, w2, 1., t2);
if(wq > 0.) return wq;
return 0.;
}
// 自由度nのt分布のパーセント点
double ptsub(double q, int n)
{
double eps, qe, s, w;
if(n == 1 && 0.001 <= q && q < 0.01) eps = 1.e-4;
else if (n == 2 && q < 0.0001) eps = 1.e-4;
else if (n == 1 && q < 0.001) eps = 1.e-2;
else eps = 1.e-5;
s = 10000.;
w = 0.;
for(;;) {
w += s;
if(s <= eps) return w;
if((qe = qt(w,n) - q) == 0.) return w;
if(qe < 0.) {
w -= s;
s /= 10.;
}
}
}
double pt(double q, int n)
{
double f, f1, f2, f3, f4, f5, u, u2, w, w0, w1, w2, w3, w4;
if(q < 1.e-5 || q > 1. || n < 1) {
fprintf(stderr,"Error : Illigal parameter in pt()!\n");
return 0.;
}
if(n <= 5) return ptsub(q, n);
if(q <= 5.e-3 && n <= 13) return ptsub(q, n);
f1 = 4. * (f = (double)n);
f5 = (f4 = (f3 = (f2 = f * f) * f) * f) * f;
f2 *= 96.;
f3 *= 384.;
f4 *= 92160.;
f5 *= 368640.;
u = pnorm(1. - q / 2.);
w0 = (u2 = u * u) * u;
w1 = w0 * u2;
w2 = w1 * u2;
w3 = w2 * u2;
w4 = w3 * u2;
w = ((w0 + u) / f1);
w += ((5. * w1 + 16. * w0 + 3. * u) / f2);
w += ((3. * w2 + 19. * w1 + 17. * w0 - 15. * u) / f3);
w += ((79. * w3 + 776. * w2 + 1482. * w1 - 1920. * w0 - 945. * u) / f4);
w += ((27. * w4 + 339. * w3 + 930. * w2 - 1782. * w1 - 765. * w0
+ 17955. * u) / f5);
return u + w;
}
int main(void) {
double m_target, n, m, s; cin >> m_target >> n >> m >> s;
double t = (m - m_target) / s * sqrt(n-1);
if (abs(t) < abs(pt(0.01, n-1))) {
cout << "YES" << endl;
} else {
cout << "NO" << endl;
}
cout << abs(t) << " " << abs(pt(0.01, n-1)) << endl;
return 0;
}
E - ゴルフボール
問題文
問題
ある会社の製品であるゴルフボールからn個抽出して重さを測定したところ，標本平均m[g]，標本単純標準偏差s[g]であった．この製品の母平均のうち、この条件で検定して5%有意となる範囲を出力せよ。
入力
n m s
n <= 30
0 <= m <= 50
0 < s <= 1
出力
l r
lは題意の範囲の下限値、rは題意の範囲の上限値
サンプル1
入力
12 45.8 1
出力
45.13 46.46
解説
t分布の信頼区間。
解答
#include <bits/stdc++.h>
#include <sys/time.h>
using namespace std;
#define rep(i,n) for(long long i = 0; i < (long long)(n); i++)
#define repi(i,a,b) for(long long i = (long long)(a); i < (long long)(b); i++)
#define pb push_back
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define mt make_tuple
#define mp make_pair
template<class T1, class T2> bool chmin(T1 &a, T2 b) { return b < a && (a = b, true); }
template<class T1, class T2> bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
#define exists find_if
#define forall all_of
using ll = long long; using vll = vector<ll>; using vvll = vector<vll>; using P = pair<ll, ll>;
using ld = double; using vld = vector<ld>;
using vi = vector<int>; using vvi = vector<vi>; vll conv(vi& v) { vll r(v.size()); rep(i, v.size()) r[i] = v[i]; return r; }
using Pos = complex<double>;
template <typename T, typename U> ostream &operator<<(ostream &o, const pair<T, U> &v) { o << "(" << v.first << ", " << v.second << ")"; return o; }
template<size_t...> struct seq{}; template<size_t N, size_t... Is> struct gen_seq : gen_seq<N-1, N-1, Is...>{}; template<size_t... Is> struct gen_seq<0, Is...> : seq<Is...>{};
template<class Ch, class Tr, class Tuple, size_t... Is>
void print_tuple(basic_ostream<Ch,Tr>& os, Tuple const& t, seq<Is...>){ using s = int[]; (void)s{0, (void(os << (Is == 0? "" : ", ") << get<Is>(t)), 0)...}; }
template<class Ch, class Tr, class... Args>
auto operator<<(basic_ostream<Ch, Tr>& os, tuple<Args...> const& t) -> basic_ostream<Ch, Tr>& { os << "("; print_tuple(os, t, gen_seq<sizeof...(Args)>()); return os << ")"; }
ostream &operator<<(ostream &o, const vvll &v) { rep(i, v.size()) { rep(j, v[i].size()) o << v[i][j] << " "; o << endl; } return o; }
template <typename T> ostream &operator<<(ostream &o, const vector<T> &v) { o << '['; rep(i, v.size()) o << v[i] << (i != v.size()-1 ? ", " : ""); o << "]"; return o; }
template <typename T> ostream &operator<<(ostream &o, const set<T> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U> ostream &operator<<(ostream &o, const map<T, U> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U, typename V> ostream &operator<<(ostream &o, const unordered_map<T, U, V> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it; o << "]"; return o; }
vector<int> range(const int x, const int y) { vector<int> v(y - x + 1); iota(v.begin(), v.end(), x); return v; }
template <typename T> istream& operator>>(istream& i, vector<T>& o) { rep(j, o.size()) i >> o[j]; return i;}
string bits_to_string(ll input, ll n=64) { string s; rep(i, n) s += '0' + !!(input & (1ll << i)); reverse(all(s)); return s; }
template <typename T> unordered_map<T, ll> counter(vector<T> vec){unordered_map<T, ll> ret; for (auto&& x : vec) ret[x]++; return ret;};
string substr(string s, P x) {return s.substr(x.fi, x.se - x.fi); }
struct ci : public iterator<forward_iterator_tag, ll> { ll n; ci(const ll n) : n(n) { } bool operator==(const ci& x) { return n == x.n; } bool operator!=(const ci& x) { return !(*this == x); } ci &operator++() { n++; return *this; } ll operator*() const { return n; } };
size_t random_seed; namespace std { using argument_type = P; template<> struct hash<argument_type> { size_t operator()(argument_type const& x) const { size_t seed = random_seed; seed ^= hash<ll>{}(x.fi); seed ^= (hash<ll>{}(x.se) << 1); return seed; } }; }; // hash for various class
namespace myhash{ const int Bsizes[]={3,9,13,17,21,25,29,33,37,41,45,49,53,57,61,65,69,73,77,81}; const int xor_nums[]={0x100007d1,0x5ff049c9,0x14560859,0x07087fef,0x3e277d49,0x4dba1f17,0x709c5988,0x05904258,0x1aa71872,0x238819b3,0x7b002bb7,0x1cf91302,0x0012290a,0x1083576b,0x76473e49,0x3d86295b,0x20536814,0x08634f4d,0x115405e8,0x0e6359f2}; const int hash_key=xor_nums[rand()%20]; const int mod_key=xor_nums[rand()%20]; template <typename T> struct myhash{ std::size_t operator()(const T& val) const { return (hash<T>{}(val)%mod_key)^hash_key; } }; };
template <typename T> class uset:public std::unordered_set<T,myhash::myhash<T>> { using SET=std::unordered_set<T,myhash::myhash<T>>; public: uset():SET(){SET::rehash(myhash::Bsizes[rand()%20]);} };
template <typename T,typename U> class umap:public std::unordered_map<T,U,myhash::myhash<T>> { public: using MAP=std::unordered_map<T,U,myhash::myhash<T>>; umap():MAP(){MAP::rehash(myhash::Bsizes[rand()%20]);} };
struct timeval start; double sec() { struct timeval tv; gettimeofday(&tv, NULL); return (tv.tv_sec - start.tv_sec) + (tv.tv_usec - start.tv_usec) * 1e-6; }
struct init_{init_(){ gettimeofday(&start, NULL); ios::sync_with_stdio(false); cin.tie(0); srand((unsigned int)time(NULL)); random_seed = RAND_MAX / 2 + rand() / 2; }} init__;
static const double EPS = 1e-14;
static const long long INF = 1e18;
static const long long mo = 1e9+7;
#define ldout fixed << setprecision(40)
// 正規分布のパーセント点
// 戸田の近似式
double pnorm(double qn)
{
static double b[11] = { 1.570796288, 0.03706987906, -0.8364353589e-3,
-0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
-0.104527497e-5, 0.8360937017e-7,-0.3231081277e-8,
0.3657763036e-10,0.6936233982e-12};
double w1, w3;
int i;
if(qn < 0. || 1. < qn) {
fprintf(stderr, "Error : qn <= 0 or qn >= 1 in pnorm()!\n");
return 0.;
}
if(qn == 0.5) return 0.;
w1 = qn;
if(qn > 0.5) w1 = 1. - w1;
w3 = -log(4. * w1 * (1. - w1));
w1 = b[0];
for(i = 1; i < 11; i++) w1 += (b[i] * pow(w3, (double)i));
if(qn > 0.5) return sqrt(w1 * w3);
return -sqrt(w1 * w3);
}
// 自由度nのt分布の両側確率
void qtsub(double *q, int n, double w2, double w3, double t2)
{
int i, j;
double w;
j = (n - 2) / 2;
for(i = 1; i <= j; i++) {
w = 2. * (double)i - w3;
*q -= (w2 *= (t2 * w / (1. + w)));
}
}
double qt(double t, int n)
{
double t1, t2, w1, w2, wq;
if(n < 1)
{
fprintf(stderr,"Error : n < 1 in qt()!\n");
return 0.;
}
w1 = 0.636619772284456;
if(t < 0.) t = - t;
t1 = t / sqrt((double)n);
t2 = 1. / (1. + t1 * t1);
if((n % 2) != 0) {
wq = 1. - w1 * atan(t1);
if(n != 1) {
wq -= (w2 = w1 * t1 * t2);
if (n != 3) qtsub(&wq, n, w2, 0., t2);
}
if(wq > 0.) return wq;
return 0.;
}
wq = 1. - (w2 = t1 * sqrt(t2));
if(n != 2) qtsub(&wq, n, w2, 1., t2);
if(wq > 0.) return wq;
return 0.;
}
// 自由度nのt分布のパーセント点
double ptsub(double q, int n)
{
double eps, qe, s, w;
if(n == 1 && 0.001 <= q && q < 0.01) eps = 1.e-4;
else if (n == 2 && q < 0.0001) eps = 1.e-4;
else if (n == 1 && q < 0.001) eps = 1.e-2;
else eps = 1.e-5;
s = 10000.;
w = 0.;
for(;;) {
w += s;
if(s <= eps) return w;
if((qe = qt(w,n) - q) == 0.) return w;
if(qe < 0.) {
w -= s;
s /= 10.;
}
}
}
double pt(double q, int n)
{
double f, f1, f2, f3, f4, f5, u, u2, w, w0, w1, w2, w3, w4;
if(q < 1.e-5 || q > 1. || n < 1) {
fprintf(stderr,"Error : Illigal parameter in pt()!\n");
return 0.;
}
if(n <= 5) return ptsub(q, n);
if(q <= 5.e-3 && n <= 13) return ptsub(q, n);
f1 = 4. * (f = (double)n);
f5 = (f4 = (f3 = (f2 = f * f) * f) * f) * f;
f2 *= 96.;
f3 *= 384.;
f4 *= 92160.;
f5 *= 368640.;
u = pnorm(1. - q / 2.);
w0 = (u2 = u * u) * u;
w1 = w0 * u2;
w2 = w1 * u2;
w3 = w2 * u2;
w4 = w3 * u2;
w = ((w0 + u) / f1);
w += ((5. * w1 + 16. * w0 + 3. * u) / f2);
w += ((3. * w2 + 19. * w1 + 17. * w0 - 15. * u) / f3);
w += ((79. * w3 + 776. * w2 + 1482. * w1 - 1920. * w0 - 945. * u) / f4);
w += ((27. * w4 + 339. * w3 + 930. * w2 - 1782. * w1 - 765. * w0
+ 17955. * u) / f5);
return u + w;
}
int main(void) {
double n, m, s; cin >> n >> m >> s;
double t = pt(0.05, n-1);
double delta = t / sqrt(n-1) * s;
cout << m - delta << " " << m + delta << endl;
return 0;
}
F - サボり教師
問題文
問題
多数の答案の中からn枚を無作為抽出して採点した結果、点数はa_iであった。試験の点数は正規分布する。この試験で母平均の信頼度95%の信頼区間を求めよ．
入力
n
a_1 ... a_n
出力
l r
サンプル1
入力
9
38 44 70 53 50 34 51 49 58
出力
41.4 57.9
解説
正規分布で信頼区間
解答
#include <bits/stdc++.h>
#include <sys/time.h>
using namespace std;
#define rep(i,n) for(long long i = 0; i < (long long)(n); i++)
#define repi(i,a,b) for(long long i = (long long)(a); i < (long long)(b); i++)
#define pb push_back
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define mt make_tuple
#define mp make_pair
template<class T1, class T2> bool chmin(T1 &a, T2 b) { return b < a && (a = b, true); }
template<class T1, class T2> bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
#define exists find_if
#define forall all_of
using ll = long long; using vll = vector<ll>; using vvll = vector<vll>; using P = pair<ll, ll>;
using ld = double; using vld = vector<ld>;
using vi = vector<int>; using vvi = vector<vi>; vll conv(vi& v) { vll r(v.size()); rep(i, v.size()) r[i] = v[i]; return r; }
using Pos = complex<double>;
template <typename T, typename U> ostream &operator<<(ostream &o, const pair<T, U> &v) { o << "(" << v.first << ", " << v.second << ")"; return o; }
template<size_t...> struct seq{}; template<size_t N, size_t... Is> struct gen_seq : gen_seq<N-1, N-1, Is...>{}; template<size_t... Is> struct gen_seq<0, Is...> : seq<Is...>{};
template<class Ch, class Tr, class Tuple, size_t... Is>
void print_tuple(basic_ostream<Ch,Tr>& os, Tuple const& t, seq<Is...>){ using s = int[]; (void)s{0, (void(os << (Is == 0? "" : ", ") << get<Is>(t)), 0)...}; }
template<class Ch, class Tr, class... Args>
auto operator<<(basic_ostream<Ch, Tr>& os, tuple<Args...> const& t) -> basic_ostream<Ch, Tr>& { os << "("; print_tuple(os, t, gen_seq<sizeof...(Args)>()); return os << ")"; }
ostream &operator<<(ostream &o, const vvll &v) { rep(i, v.size()) { rep(j, v[i].size()) o << v[i][j] << " "; o << endl; } return o; }
template <typename T> ostream &operator<<(ostream &o, const vector<T> &v) { o << '['; rep(i, v.size()) o << v[i] << (i != v.size()-1 ? ", " : ""); o << "]"; return o; }
template <typename T> ostream &operator<<(ostream &o, const set<T> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U> ostream &operator<<(ostream &o, const map<T, U> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U, typename V> ostream &operator<<(ostream &o, const unordered_map<T, U, V> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it; o << "]"; return o; }
vector<int> range(const int x, const int y) { vector<int> v(y - x + 1); iota(v.begin(), v.end(), x); return v; }
template <typename T> istream& operator>>(istream& i, vector<T>& o) { rep(j, o.size()) i >> o[j]; return i;}
string bits_to_string(ll input, ll n=64) { string s; rep(i, n) s += '0' + !!(input & (1ll << i)); reverse(all(s)); return s; }
template <typename T> unordered_map<T, ll> counter(vector<T> vec){unordered_map<T, ll> ret; for (auto&& x : vec) ret[x]++; return ret;};
string substr(string s, P x) {return s.substr(x.fi, x.se - x.fi); }
struct ci : public iterator<forward_iterator_tag, ll> { ll n; ci(const ll n) : n(n) { } bool operator==(const ci& x) { return n == x.n; } bool operator!=(const ci& x) { return !(*this == x); } ci &operator++() { n++; return *this; } ll operator*() const { return n; } };
size_t random_seed; namespace std { using argument_type = P; template<> struct hash<argument_type> { size_t operator()(argument_type const& x) const { size_t seed = random_seed; seed ^= hash<ll>{}(x.fi); seed ^= (hash<ll>{}(x.se) << 1); return seed; } }; }; // hash for various class
namespace myhash{ const int Bsizes[]={3,9,13,17,21,25,29,33,37,41,45,49,53,57,61,65,69,73,77,81}; const int xor_nums[]={0x100007d1,0x5ff049c9,0x14560859,0x07087fef,0x3e277d49,0x4dba1f17,0x709c5988,0x05904258,0x1aa71872,0x238819b3,0x7b002bb7,0x1cf91302,0x0012290a,0x1083576b,0x76473e49,0x3d86295b,0x20536814,0x08634f4d,0x115405e8,0x0e6359f2}; const int hash_key=xor_nums[rand()%20]; const int mod_key=xor_nums[rand()%20]; template <typename T> struct myhash{ std::size_t operator()(const T& val) const { return (hash<T>{}(val)%mod_key)^hash_key; } }; };
template <typename T> class uset:public std::unordered_set<T,myhash::myhash<T>> { using SET=std::unordered_set<T,myhash::myhash<T>>; public: uset():SET(){SET::rehash(myhash::Bsizes[rand()%20]);} };
template <typename T,typename U> class umap:public std::unordered_map<T,U,myhash::myhash<T>> { public: using MAP=std::unordered_map<T,U,myhash::myhash<T>>; umap():MAP(){MAP::rehash(myhash::Bsizes[rand()%20]);} };
struct timeval start; double sec() { struct timeval tv; gettimeofday(&tv, NULL); return (tv.tv_sec - start.tv_sec) + (tv.tv_usec - start.tv_usec) * 1e-6; }
struct init_{init_(){ gettimeofday(&start, NULL); ios::sync_with_stdio(false); cin.tie(0); srand((unsigned int)time(NULL)); random_seed = RAND_MAX / 2 + rand() / 2; }} init__;
static const double EPS = 1e-14;
static const long long INF = 1e18;
static const long long mo = 1e9+7;
#define ldout fixed << setprecision(40)
// 正規分布のパーセント点
// 戸田の近似式
double pnorm(double qn)
{
static double b[11] = { 1.570796288, 0.03706987906, -0.8364353589e-3,
-0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
-0.104527497e-5, 0.8360937017e-7,-0.3231081277e-8,
0.3657763036e-10,0.6936233982e-12};
double w1, w3;
int i;
if(qn < 0. || 1. < qn) {
fprintf(stderr, "Error : qn <= 0 or qn >= 1 in pnorm()!\n");
return 0.;
}
if(qn == 0.5) return 0.;
w1 = qn;
if(qn > 0.5) w1 = 1. - w1;
w3 = -log(4. * w1 * (1. - w1));
w1 = b[0];
for(i = 1; i < 11; i++) w1 += (b[i] * pow(w3, (double)i));
if(qn > 0.5) return sqrt(w1 * w3);
return -sqrt(w1 * w3);
}
// 自由度nのt分布の両側確率
void qtsub(double *q, int n, double w2, double w3, double t2)
{
int i, j;
double w;
j = (n - 2) / 2;
for(i = 1; i <= j; i++) {
w = 2. * (double)i - w3;
*q -= (w2 *= (t2 * w / (1. + w)));
}
}
double qt(double t, int n)
{
double t1, t2, w1, w2, wq;
if(n < 1)
{
fprintf(stderr,"Error : n < 1 in qt()!\n");
return 0.;
}
w1 = 0.636619772284456;
if(t < 0.) t = - t;
t1 = t / sqrt((double)n);
t2 = 1. / (1. + t1 * t1);
if((n % 2) != 0) {
wq = 1. - w1 * atan(t1);
if(n != 1) {
wq -= (w2 = w1 * t1 * t2);
if (n != 3) qtsub(&wq, n, w2, 0., t2);
}
if(wq > 0.) return wq;
return 0.;
}
wq = 1. - (w2 = t1 * sqrt(t2));
if(n != 2) qtsub(&wq, n, w2, 1., t2);
if(wq > 0.) return wq;
return 0.;
}
// 自由度nのt分布のパーセント点
double ptsub(double q, int n)
{
double eps, qe, s, w;
if(n == 1 && 0.001 <= q && q < 0.01) eps = 1.e-4;
else if (n == 2 && q < 0.0001) eps = 1.e-4;
else if (n == 1 && q < 0.001) eps = 1.e-2;
else eps = 1.e-5;
s = 10000.;
w = 0.;
for(;;) {
w += s;
if(s <= eps) return w;
if((qe = qt(w,n) - q) == 0.) return w;
if(qe < 0.) {
w -= s;
s /= 10.;
}
}
}
double pt(double q, int n)
{
double f, f1, f2, f3, f4, f5, u, u2, w, w0, w1, w2, w3, w4;
if(q < 1.e-5 || q > 1. || n < 1) {
fprintf(stderr,"Error : Illigal parameter in pt()!\n");
return 0.;
}
if(n <= 5) return ptsub(q, n);
if(q <= 5.e-3 && n <= 13) return ptsub(q, n);
f1 = 4. * (f = (double)n);
f5 = (f4 = (f3 = (f2 = f * f) * f) * f) * f;
f2 *= 96.;
f3 *= 384.;
f4 *= 92160.;
f5 *= 368640.;
u = pnorm(1. - q / 2.);
w0 = (u2 = u * u) * u;
w1 = w0 * u2;
w2 = w1 * u2;
w3 = w2 * u2;
w4 = w3 * u2;
w = ((w0 + u) / f1);
w += ((5. * w1 + 16. * w0 + 3. * u) / f2);
w += ((3. * w2 + 19. * w1 + 17. * w0 - 15. * u) / f3);
w += ((79. * w3 + 776. * w2 + 1482. * w1 - 1920. * w0 - 945. * u) / f4);
w += ((27. * w4 + 339. * w3 + 930. * w2 - 1782. * w1 - 765. * w0
+ 17955. * u) / f5);
return u + w;
}
int main(void) {
int n; cin >> n;
vector<double> a(n);
rep(i, n) cin >> a[i];
double m = accumulate(all(a), 0.0) / n;
// 単純標準偏差
double s = 0;
rep(i, n)
s += pow(a[i] - m, 2);
s /= n;
s = sqrt(s);
double delta = s / sqrt(n-1) * pt(0.05, 1.0*n-1);
cout << m - delta << " " << m + delta << endl;
return 0;
}
G - 電話の多い日
問題文
問題
ある会社では、1時間あたり平均n回電話がかかってくる。1日の業務時間8時間あたりm回以上電話がかかってくる確率を求めよ。
入力
n m
出力
確率を1行で。
サンプル1
入力
2 30
出力
YES
解説
ランダムに起きると仮定してポワソン分布を数値積分
解答
#include <bits/stdc++.h>
#include <sys/time.h>
using namespace std;
#define rep(i,n) for(long long i = 0; i < (long long)(n); i++)
#define repi(i,a,b) for(long long i = (long long)(a); i < (long long)(b); i++)
#define pb push_back
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define mt make_tuple
#define mp make_pair
template<class T1, class T2> bool chmin(T1 &a, T2 b) { return b < a && (a = b, true); }
template<class T1, class T2> bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
#define exists find_if
#define forall all_of
using ll = long long; using vll = vector<ll>; using vvll = vector<vll>; using P = pair<ll, ll>;
using ld = double; using vld = vector<ld>;
using vi = vector<int>; using vvi = vector<vi>; vll conv(vi& v) { vll r(v.size()); rep(i, v.size()) r[i] = v[i]; return r; }
using Pos = complex<double>;
template <typename T, typename U> ostream &operator<<(ostream &o, const pair<T, U> &v) { o << "(" << v.first << ", " << v.second << ")"; return o; }
template<size_t...> struct seq{}; template<size_t N, size_t... Is> struct gen_seq : gen_seq<N-1, N-1, Is...>{}; template<size_t... Is> struct gen_seq<0, Is...> : seq<Is...>{};
template<class Ch, class Tr, class Tuple, size_t... Is>
void print_tuple(basic_ostream<Ch,Tr>& os, Tuple const& t, seq<Is...>){ using s = int[]; (void)s{0, (void(os << (Is == 0? "" : ", ") << get<Is>(t)), 0)...}; }
template<class Ch, class Tr, class... Args>
auto operator<<(basic_ostream<Ch, Tr>& os, tuple<Args...> const& t) -> basic_ostream<Ch, Tr>& { os << "("; print_tuple(os, t, gen_seq<sizeof...(Args)>()); return os << ")"; }
ostream &operator<<(ostream &o, const vvll &v) { rep(i, v.size()) { rep(j, v[i].size()) o << v[i][j] << " "; o << endl; } return o; }
template <typename T> ostream &operator<<(ostream &o, const vector<T> &v) { o << '['; rep(i, v.size()) o << v[i] << (i != v.size()-1 ? ", " : ""); o << "]"; return o; }
template <typename T> ostream &operator<<(ostream &o, const set<T> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U> ostream &operator<<(ostream &o, const map<T, U> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]"; return o; }
template <typename T, typename U, typename V> ostream &operator<<(ostream &o, const unordered_map<T, U, V> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it; o << "]"; return o; }
vector<int> range(const int x, const int y) { vector<int> v(y - x + 1); iota(v.begin(), v.end(), x); return v; }
template <typename T> istream& operator>>(istream& i, vector<T>& o) { rep(j, o.size()) i >> o[j]; return i;}
string bits_to_string(ll input, ll n=64) { string s; rep(i, n) s += '0' + !!(input & (1ll << i)); reverse(all(s)); return s; }
template <typename T> unordered_map<T, ll> counter(vector<T> vec){unordered_map<T, ll> ret; for (auto&& x : vec) ret[x]++; return ret;};
string substr(string s, P x) {return s.substr(x.fi, x.se - x.fi); }
struct ci : public iterator<forward_iterator_tag, ll> { ll n; ci(const ll n) : n(n) { } bool operator==(const ci& x) { return n == x.n; } bool operator!=(const ci& x) { return !(*this == x); } ci &operator++() { n++; return *this; } ll operator*() const { return n; } };
size_t random_seed; namespace std { using argument_type = P; template<> struct hash<argument_type> { size_t operator()(argument_type const& x) const { size_t seed = random_seed; seed ^= hash<ll>{}(x.fi); seed ^= (hash<ll>{}(x.se) << 1); return seed; } }; }; // hash for various class
namespace myhash{ const int Bsizes[]={3,9,13,17,21,25,29,33,37,41,45,49,53,57,61,65,69,73,77,81}; const int xor_nums[]={0x100007d1,0x5ff049c9,0x14560859,0x07087fef,0x3e277d49,0x4dba1f17,0x709c5988,0x05904258,0x1aa71872,0x238819b3,0x7b002bb7,0x1cf91302,0x0012290a,0x1083576b,0x76473e49,0x3d86295b,0x20536814,0x08634f4d,0x115405e8,0x0e6359f2}; const int hash_key=xor_nums[rand()%20]; const int mod_key=xor_nums[rand()%20]; template <typename T> struct myhash{ std::size_t operator()(const T& val) const { return (hash<T>{}(val)%mod_key)^hash_key; } }; };
template <typename T> class uset:public std::unordered_set<T,myhash::myhash<T>> { using SET=std::unordered_set<T,myhash::myhash<T>>; public: uset():SET(){SET::rehash(myhash::Bsizes[rand()%20]);} };
template <typename T,typename U> class umap:public std::unordered_map<T,U,myhash::myhash<T>> { public: using MAP=std::unordered_map<T,U,myhash::myhash<T>>; umap():MAP(){MAP::rehash(myhash::Bsizes[rand()%20]);} };
struct timeval start; double sec() { struct timeval tv; gettimeofday(&tv, NULL); return (tv.tv_sec - start.tv_sec) + (tv.tv_usec - start.tv_usec) * 1e-6; }
struct init_{init_(){ gettimeofday(&start, NULL); ios::sync_with_stdio(false); cin.tie(0); srand((unsigned int)time(NULL)); random_seed = RAND_MAX / 2 + rand() / 2; }} init__;
static const double EPS = 1e-14;
static const long long INF = 1e18;
static const long long mo = 1e9+7;
#define ldout fixed << setprecision(40)
// tgammaを使わないとだめ。あとlgammaなる関数がある
double poisson(double k, double lamda) {
return pow(lamda, k) * exp(-lamda) / tgamma(k+1);
}
int main(void) {
ll n, m; cin >> n >> m;
n *= 8;
double ret = 1;
rep(i, m) {
ret -= poisson(i, n);
}
cout << ret << endl;
return 0;
}