| 10 | |
| 11 | |
| 12 | def test_tutorial(): |
| 13 | def model(params, t): |
| 14 | _, _, amp, loc, sig2 = params |
| 15 | return amp * np.exp(-0.5 * (t - loc) ** 2 / sig2) |
| 16 | |
| 17 | def lnlike(p, t, y, yerr, solver=BasicSolver): |
| 18 | a, tau = np.exp(p[:2]) |
| 19 | gp = GP(a * kernels.Matern32Kernel(tau) + 0.001, solver=solver) |
| 20 | gp.compute(t, yerr) |
| 21 | return gp.lnlikelihood(y - model(p, t)) |
| 22 | |
| 23 | def lnprior(p): |
| 24 | lna, lntau, amp, loc, sig2 = p |
| 25 | if (-5 < lna < 5 and -5 < lntau < 5 and -10 < amp < 10 and |
| 26 | -5 < loc < 5 and 0 < sig2 < 3): |
| 27 | return 0.0 |
| 28 | return -np.inf |
| 29 | |
| 30 | def lnprob(p, x, y, yerr, **kwargs): |
| 31 | lp = lnprior(p) |
| 32 | return lp + lnlike(p, x, y, yerr, **kwargs) \ |
| 33 | if np.isfinite(lp) else -np.inf |
| 34 | |
| 35 | np.random.seed(1234) |
| 36 | x = np.sort(np.random.rand(50)) |
| 37 | yerr = 0.05 + 0.01 * np.random.rand(len(x)) |
| 38 | y = np.sin(x) + yerr * np.random.randn(len(x)) |
| 39 | p = [0, 0, -1.0, 0.1, 0.4] |
| 40 | assert np.isfinite(lnprob(p, x, y, yerr)), "Incorrect result" |
| 41 | assert np.allclose(lnprob(p, x, y, yerr), |
| 42 | lnprob(p, x, y, yerr, solver=HODLRSolver)), \ |
| 43 | "Inconsistent results" |