Last active
November 2, 2025 06:32
-
-
Save rhaps0dy/7ef5f4f4cbfa8b2e795f51a847f50e87 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import math | |
| import numpy as np | |
| def logsumexp(a, b): | |
| """Stable log(exp(a) + exp(b))""" | |
| m = np.maximum(a, b) | |
| return m + np.log(np.exp(a - m) + np.exp(b - m)) | |
| def solve(N, num_points=10000): | |
| """ | |
| Find epsilon minimizing f(eps) = eps + exp(-2*N*eps^2) | |
| using grid search in log-space over log_epsilon ∈ [-10, 10]. | |
| Computes in log-space using logsumexp for stability. | |
| """ | |
| # grid of log_epsilon values | |
| log_eps_grid = np.linspace(-10.0, 10.0, num_points) | |
| eps_grid = np.exp(log_eps_grid) | |
| # compute log(f(eps)) = log( eps + exp(-2*N*eps^2) ) | |
| log_term1 = log_eps_grid | |
| log_term2 = -2.0 * N * eps_grid * eps_grid | |
| log_f_grid = logsumexp(log_term1, log_term2) | |
| # find minimizing index | |
| idx_min = np.argmin(log_f_grid) | |
| log_eps_star = log_eps_grid[idx_min] | |
| eps_star = eps_grid[idx_min] | |
| log_val = log_f_grid[idx_min] | |
| val = np.exp(log_val) | |
| return eps_star, val | |
| if __name__ == "__main__": | |
| Ns = [10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000] | |
| for N in Ns: | |
| eps, val = solve(N) | |
| print(f"N=10^{int(math.log10(N))} eps={eps:.3e} f(eps)={val:.6f}") |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import numpy as np | |
| # Parameters | |
| N = 30 # sheep per group | |
| sigma = 20 # std deviation in cm | |
| diff_obs = 82.1 - 77.2 # observed difference | |
| S = 1_000_000 # number of simulations | |
| # Simulate null: no real difference (mean=0) | |
| sim_diffs = np.mean(np.random.normal(0, sigma, (S, N)), axis=1) - np.mean( | |
| np.random.normal(0, sigma, (S, N)), axis=1 | |
| ) | |
| # Two-sided p-value | |
| p = np.mean(np.abs(sim_diffs) >= abs(diff_obs)) | |
| print(f"Observed difference: {diff_obs:.2f} cm") | |
| print(f"p-value ≈ {p:.3g}") | |
| if p < 0.05: | |
| print("Reject the null hypothesis: heights differ significantly.") | |
| else: | |
| print("Fail to reject the null hypothesis: no significant difference.") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment