b419. 公平的硬幣

1. Problem
2. Sample Input
3. Sample Output
4. Solution
1. 4.1. 暴力檢測
2. 4.2. 數值積分

Problem

一天，liouzhou_101 去打印店裏打印了幾張複習資料，花了 4 毛錢，他給了店主1張塊錢的鈔票，結果店主補了他兩枚硬幣，一個5毛錢，一個 1 毛錢。

在走回宿舍的途中，liouzhou_101 一直在想一個問題，5 毛錢的硬幣和 1 毛錢的硬幣哪個更公平？所謂公平，就是指將一枚硬幣拋擲 1 次，他正面朝上的概率是 12，反面朝上的概率也是 12。

結果他一會去就把1毛錢的硬幣拋了 1000 次，記下來有 531 次正面朝上。然後他突然說道：“這硬幣不公平！！！”

真的不公平嗎？他希望你計算一下出現這種情況的概率。

Sample Input

10 5 5
20 6 8
100 40 60
1000 531 531
5000 2345 2456

Sample Output

Solution

跟柳柳州互相出題目，已經來到了第 n 天，儘管我出的都是垃圾跟搬運題目，一出題被電得好愉悅，出在多的垃圾題都電不到柳柳州。

在擲硬幣找區間次數 $[a, b]$ 的機率，由於是離散的統計上就只能推組合數學，但又不是曲棍棒的模型 (Hockey Stick Pattern)，算了先手動窮舉，用 Stirling’s formula 計算組合，後來發現有幾組數據有問題。

當 n 太大的時候，機率分布就相當稀疏，直接套用數值積分，直接去積分自然分布得到公式稍微困難，利用 simpson 或者是 romberg 去解決，但這樣還會有一個問題，只有在平均值的部分機率非常高 (類似離散的 00000100000)，導致自適應辛普森積分會判斷錯誤，提醒剖半積分後測資就正確了不少。似乎合作出了一個可怕的題目，儘管如此，還是得用 mathlab 或者是 wolframalpha 來協助測資的檢驗，剩下的部分就交給柳柳州了。

關於自然分布 (normal distribution)：

骰子擲 $n$ 次，命中的機率 $p = 1/2$
期望值 $\mu = np = n/2$
標準差為 $\sigma = \sqrt{np(1-p)} = \sqrt{n/4}$。

期望值 $E[x] = \sum_{k=1}^{n} k \binom{n}{k} p^k (1-p)^{n-k}$，消階層之後提出，計算後得到 $E[x] = np$。

而標準差 $\sigma$ 先從變異數下手

定義： $Var[x] = E[(x - \mu)^2] = \sum (x - \mu)^2 P[X = x]$
拆解步驟 1： $Var[x] = \sum (x^2 - 2 \mu x + \mu^{2}) P[X = x]$
拆解步驟 2： $Var[x] = E[x^2] - 2 E[x] E[x] + E[x] E[x] = E[x^2] - E[x]^2$
拆解步驟 3： $E[x^2] = E[x(x-1)] + E[x] = n(n-1)p^2 + np$
拆解步驟 4： $Var[x] = np(1-p)$

接下來直接帶入自然分布的公式，參數 $(\mu, \sigma)$，特別小心使用辛普森積分 $[a, b]$ 時，由於問題是離散的，要改變成 $[a - 0.5, b + 0.5]$，防止 $a = b$ 時造成積分錯誤。

特別感謝學弟 firejox 的數學指導。

不久之後，firejox 說道內建函數 erf() 可以解決積分問題，詳見 wiki Gauss error function，因此這一題也可以不寫辛普森積分，套用內建函數在不到 10 行內解決此題。

暴力檢測

#include <bits/stdc++.h>
using namespace std;
const long double PI = acos(-1), e = exp(1);
const int MAXN = 100001;
double f[MAXN];
double stirling(int n) {
    if (n < MAXN)
        return f[n];
    return log2(sqrt(2*PI*n)) + n*log2(n/e);
}
double logC(int n, int m) {
    return stirling(n) - stirling(n - m) - stirling(m);
}
int main() {
    f[0] = 0;
    for (int i = 1; i < MAXN; i++)
        f[i] = log2(i) + f[i-1];
    int n, a, b;
    int cases = 0;
    while (scanf("%d %d %d", &n, &a, &b) == 3) {
        double ret = 0;
        for (int i = a; i <= b; i++)
            ret += pow(2, logC(n, i) - log2(2) * n);
        printf("%.4lf\n", ret);
        if (++cases > 100)
            break;
    }
    return 0;
}

數值積分

#include <bits/stdc++.h>
using namespace std;
const long double PI = acos(-1), e = exp(1);
const int MAXN = 1001;
double f[MAXN];
double stirling(int n) {
    if (n < MAXN)
        return f[n];
    return log2(sqrt(2*PI*n)) + n*log2(n/e);
}
double logC(int n, int m) {
    return stirling(n) - stirling(n - m) - stirling(m);
}
class SimpsonIntegral {
public:
    const double eps = 1e-6;
    const double pi = acos(-1);
    const double sqr2pi = sqrt(2 * pi);
    double mu, sigma; // function coefficient
    double f(double x) { // usually replace
        return exp(-(pow(x - mu, 2)/2.0/sigma/sigma))/ sigma / sqr2pi;
    }
    void setCoefficient(double m, double s) {
        this->mu = m, this->sigma = s;
    } 
    double simpson(double a, double b, double fa, double fb, double fab) {
        return (fa + 4 * fab + fb) * (b - a) / 6.0;
    }
    double integral(double l, double r) {
        return integral(l, r, eps);
    }
private:
    double integral(double a, double b, double c, double eps, double A, double fa, double fb, double fc) {
        double ab = (a + b)/2, bc = (b + c)/2;
        double fab = f(ab), fbc = f(bc);
        double L = simpson(a, b, fa, fb, fab), R = simpson(b, c, fb, fc, fbc);
        if (fabs(L + R - A) <= 15 * eps)
            return L + R + (L + R - A) / 15.0;
        return integral(a, ab, b, eps/2, L, fa, fab, fb) + integral(b, bc, c, eps/2, R, fb, fbc, fc);
    }
    double integral(double a, double b, double eps) {
        double ab = (a + b) /2;
        double fa = f(a), fb = f(b), fab = f(ab);
        return integral(a, ab, b, eps, simpson(a, b, fa, fb, fab), fa, fab, fb);
    }
} tool;
double binom(double n, double a, double b) {
    double a1, b1;
    double d = sqrt(n/2.0);
    a1 = 0.5 * (1.0 + erf((a - n/2.0 - 0.5) / d));
    b1 = 0.5 * (1.0 + erf((b - n/2.0 + 0.5) / d));
    return b1 - a1;
}
int main() {
    f[0] = 0;
    for (int i = 1; i < MAXN; i++)
        f[i] = log2(i) + f[i-1];
    int n, a, b;
    int cases = 0;
    while (scanf("%d %d %d", &n, &a, &b) == 3) {
        cases++;
        if (n < MAXN) {
            double ret = 0;
            for (int i = a; i <= b; i++)
                ret += pow(2, logC(n, i) - log2(2) * n);
            printf("%.4lf\n", ret);
        } else {
//			tool.setCoefficient(n/2.0, sqrt(n/4.0));
//			double ret = 0, l = a - 0.5, r = b + 0.5;
//			if (r < n/2.0 || l > n/2.0)
//				ret = tool.integral(l, r);
//			else
//				ret = tool.integral(l, n/2.0) + tool.integral(n/2.0, r);
//			printf("%.4lf\n", ret);
            printf("%.4lf\n", binom(n, a, b));
        } 
    }
    return 0;
}

Morris' Blog