Hypothesis Testing

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Compute standard normal CDF using erfc (numerically stable)
5
// Phi(z) = 0.5 * erfc(-z / sqrt(2))
6
static inline double normal_cdf(double z) {
7
    return 0.5 * erfc(-z / M_SQRT2); // M_SQRT2 = sqrt(2)
8
}
9
10
struct ZTestResult {
11
    double mean;
12
    double z;
13
    double p_two_sided;
14
};
15
16
ZTestResult z_test_mean_known_sigma(const vector<double>& x, double mu0, double sigma) {
17
    if (x.empty()) throw invalid_argument("Sample must be non-empty");
18
    if (!(sigma > 0)) throw invalid_argument("Sigma must be positive");
19
20
    // Compute sample mean in O(n)
21
    double sum = accumulate(x.begin(), x.end(), 0.0);
22
    double mean = sum / x.size();
23
24
    // Standard error and z-statistic
25
    double se = sigma / sqrt((double)x.size());
26
    double z = (mean - mu0) / se;
27
28
    // Two-sided p-value: p = 2*(1 - Phi(|z|)) = erfc(|z|/sqrt(2))
29
    double p_two = erfc(fabs(z) / M_SQRT2);
30
31
    return {mean, z, p_two};
32
}
33
34
int main() {
35
    // Example: test H0: mu = 10 with known sigma = 2
36
    vector<double> data = {9.5, 10.4, 9.8, 10.1, 10.7, 9.9, 10.3, 9.6};
37
    double mu0 = 10.0;
38
    double sigma = 2.0;
39
40
    auto res = z_test_mean_known_sigma(data, mu0, sigma);
41
42
    cout.setf(ios::fixed); cout << setprecision(6);
43
    cout << "Sample mean = " << res.mean << "\n";
44
    cout << "z-statistic = " << res.z << "\n";
45
    cout << "Two-sided p-value = " << res.p_two_sided << "\n";
46
47
    double alpha = 0.05;
48
    cout << (res.p_two_sided < alpha ? "Reject H0 at alpha=0.05" : "Fail to reject H0 at alpha=0.05") << "\n";
49
    return 0;
50
}
51

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Log of binomial PMF using lgamma for stability
5
static inline long double log_binom_pmf(int n, int k, long double p) {
6
    if (k < 0 || k > n) return -numeric_limits<long double>::infinity();
7
    if (p == 0.0L) return (k == 0) ? 0.0L : -numeric_limits<long double>::infinity();
8
    if (p == 1.0L) return (k == n) ? 0.0L : -numeric_limits<long double>::infinity();
9
    long double lnC = lgammal(n + 1.0L) - lgammal(k + 1.0L) - lgammal((n - k) + 1.0L);
10
    return lnC + k * logl(p) + (n - k) * logl(1.0L - p);
11
}
12
13
struct BinomialTestResult {
14
    int n, k;
15
    long double p0;
16
    long double p_value_two_sided;
17
};
18
19
BinomialTestResult exact_binomial_test_two_sided(int n, int k_obs, long double p0) {
20
    if (n <= 0) throw invalid_argument("n must be positive");
21
    if (k_obs < 0 || k_obs > n) throw invalid_argument("k must be in [0,n]");
22
    if (!(p0 > 0.0L && p0 < 1.0L)) throw invalid_argument("p0 must be in (0,1)");
23
24
    long double logp_obs = log_binom_pmf(n, k_obs, p0);
25
26
    // Conservative two-sided p-value: sum probabilities with pmf <= pmf_observed
27
    long double p_sum = 0.0L;
28
    for (int k = 0; k <= n; ++k) {
29
        long double lp = log_binom_pmf(n, k, p0);
30
        if (lp <= logp_obs + 1e-18L) { // tiny tolerance
31
            p_sum += expl(lp);
32
        }
33
    }
34
35
    // Clamp due to numerical rounding
36
    p_sum = min<long double>(1.0L, max<long double>(0.0L, p_sum));
37
    return {n, k_obs, p0, p_sum};
38
}
39
40
int main() {
41
    // Example: 20 trials, 7 successes, test H0: p = 0.25
42
    int n = 20, k = 7; long double p0 = 0.25L;
43
    auto res = exact_binomial_test_two_sided(n, k, p0);
44
45
    cout.setf(ios::fixed); cout << setprecision(10);
46
    cout << "n = " << res.n << ", k = " << res.k << ", p0 = " << (double)res.p0 << "\n";
47
    cout << "Two-sided exact p-value = " << (double)res.p_value_two_sided << "\n";
48
49
    double alpha = 0.05;
50
    cout << (res.p_value_two_sided < alpha ? "Reject H0" : "Fail to reject H0") << " at alpha=0.05\n";
51
    return 0;
52
}
53

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct PermTestResult {
5
    double observed_abs_diff;
6
    double p_value_two_sided;
7
};
8
9
static double mean_of(const vector<double>& v) {
10
    return accumulate(v.begin(), v.end(), 0.0) / max<size_t>(1, v.size());
11
}
12
13
PermTestResult permutation_test_diff_means(const vector<double>& A, const vector<double>& B, int B_iters, uint64_t seed = 1234567ULL) {
14
    if (A.empty() || B.empty()) throw invalid_argument("Both groups must be non-empty");
15
    if (B_iters <= 0) throw invalid_argument("Number of permutations must be positive");
16
17
    // Observed statistic: absolute difference in means
18
    double obs = fabs(mean_of(A) - mean_of(B));
19
20
    // Combine data
21
    vector<double> all; all.reserve(A.size() + B.size());
22
    all.insert(all.end(), A.begin(), A.end());
23
    all.insert(all.end(), B.begin(), B.end());
24
25
    size_t n1 = A.size(); size_t n2 = B.size();
26
27
    // Random generator
28
    mt19937_64 rng(seed);
29
30
    int extreme = 0;
31
    vector<double> tmp = all; // will be shuffled in-place
32
33
    for (int it = 0; it < B_iters; ++it) {
34
        shuffle(tmp.begin(), tmp.end(), rng);
35
36
        // Compute means of the first n1 and remaining n2 elements
37
        double sum1 = 0.0; for (size_t i = 0; i < n1; ++i) sum1 += tmp[i];
38
        double sum2 = 0.0; for (size_t i = n1; i < n1 + n2; ++i) sum2 += tmp[i];
39
        double diff = fabs((sum1 / n1) - (sum2 / n2));
40
41
        if (diff >= obs - 1e-15) ++extreme; // count as or more extreme
42
    }
43
44
    // Monte Carlo p-value with +1 correction for stability
45
    double p = (extreme + 1.0) / (B_iters + 1.0);
46
    return {obs, p};
47
}
48
49
int main() {
50
    vector<double> A = {10.1, 9.9, 10.3, 10.0, 9.8};
51
    vector<double> B = {10.7, 10.5, 10.6, 10.4};
52
53
    int permutations = 20000; // increase for more precision
54
    auto res = permutation_test_diff_means(A, B, permutations, 42);
55
56
    cout.setf(ios::fixed); cout << setprecision(6);
57
    cout << "Observed |mean(A) - mean(B)| = " << res.observed_abs_diff << "\n";
58
    cout << "Two-sided permutation p-value (Monte Carlo) = " << res.p_value_two_sided << "\n";
59
60
    double alpha = 0.05;
61
    cout << (res.p_value_two_sided < alpha ? "Reject H0 (equal means)" : "Fail to reject H0 (equal means)") << " at alpha=0.05\n";
62
    return 0;
63
}
64