Skip to content

Commit

Permalink
[pre-commit.ci] auto fixes from pre-commit.com hooks
Browse files Browse the repository at this point in the history 
for more information, see https://pre-commit.ci
  • Loading branch information
@pre-commit-ci
pre-commit-ci[bot] committed Aug 21, 2024
1 parent 505d89a commit 22787da
Showing 1 changed file with 94 additions and 21 deletions.
115 changes: 94 additions & 21 deletions doc/notebooks/weighted_histograms.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,25 +109,30 @@
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"def log_or_zero(x):\n",
" x = np.atleast_1d(x).copy()\n",
" ma = x > 0\n",
" x[ma] = np.log(x[ma])\n",
" x[~ma] = 0\n",
" return x\n",
"\n",
"\n",
"def l1(w, w2, mu):\n",
" s = np.abs(w) / w2\n",
" return 2 * s * (w * (log_or_zero(w) - log_or_zero(mu)) + mu - w)\n",
"\n",
"\n",
"def l2(w, w2, mu):\n",
" s = np.abs(w) / w2\n",
" return 2 * s * (w * log_or_zero(np.abs(w)) - w * log_or_zero(mu) + mu - w)\n",
"\n",
"\n",
"def l3(w, w2, mu):\n",
" s = w / w2\n",
" return 2 * s * (w * (log_or_zero(np.abs(w)) - log_or_zero(mu)) + mu - w)\n",
"\n",
"\n",
"w2 = 1\n",
"fig, ax = plt.subplots(1, 3, figsize=(10, 4))\n",
"for i, (axi, w) in enumerate(zip(ax, (-10, -5, 5))):\n",
Expand Down Expand Up @@ -190,13 +195,23 @@
"x = expon.rvs(size=npoints, random_state=rng)\n",
"w = norm.rvs(0.1, x, size=len(x), random_state=rng)\n",
"\n",
"h = bh.Histogram(bh.axis.Regular(bins, np.min(x), np.max(x)), storage=bh.storage.Weight())\n",
"h = bh.Histogram(\n",
" bh.axis.Regular(bins, np.min(x), np.max(x)), storage=bh.storage.Weight()\n",
")\n",
"h.fill(x, weight=w)\n",
"\n",
"plt.stairs(h.values(), h.axes[0].edges)\n",
"ma = h.values() > 0\n",
"plt.errorbar(h.axes[0].centers[ma], h.values()[ma], h.variances()[ma] ** 0.5, fmt=\"o\", color=\"C0\")\n",
"plt.errorbar(h.axes[0].centers[~ma], h.values()[~ma], h.variances()[~ma] ** 0.5, fmt=\"s\", color=\"C3\")\n",
"plt.errorbar(\n",
" h.axes[0].centers[ma], h.values()[ma], h.variances()[ma] ** 0.5, fmt=\"o\", color=\"C0\"\n",
")\n",
"plt.errorbar(\n",
" h.axes[0].centers[~ma],\n",
" h.values()[~ma],\n",
" h.variances()[~ma] ** 0.5,\n",
" fmt=\"s\",\n",
" color=\"C3\",\n",
")\n",
"plt.axhline(0, ls=\":\", color=\"0.5\", zorder=0)\n",
"xm = np.linspace(0, h.axes[0].edges[-1])\n",
"plt.plot(xm, expon.pdf(xm) * npoints * np.mean(w) * h.axes[0].widths[0]);"
Expand Down Expand Up @@ -931,6 +946,7 @@
"def model1(x, n, lambd):\n",
" return n * expon(0, lambd).cdf(x)\n",
"\n",
"\n",
"c1 = ExtendedBinnedNLL(np.transpose((n, vn)), xe, model1)\n",
"m = Minuit(c1, sum(n), 1)\n",
"m.limits = (0, None)\n",
Expand Down Expand Up @@ -1621,6 +1637,7 @@
"def model2(x, lambd):\n",
" return expon(0, lambd).cdf(x)\n",
"\n",
"\n",
"c2 = BinnedNLL(np.transpose((n, vn)), xe, model2)\n",
"m = Minuit(c2, 1)\n",
"m.limits = (0, None)\n",
Expand Down Expand Up @@ -1648,7 +1665,9 @@
" x = expon.rvs(size=rng.poisson(npoints), random_state=rng)\n",
" w = norm.rvs(0.1, x, size=len(x), random_state=rng)\n",
"\n",
" h = bh.Histogram(bh.axis.Regular(bins, np.min(x), np.max(x)), storage=bh.storage.Weight())\n",
" h = bh.Histogram(\n",
" bh.axis.Regular(bins, np.min(x), np.max(x)), storage=bh.storage.Weight()\n",
" )\n",
" h.fill(x, weight=w)\n",
" xe = h.axes[0].edges\n",
" n = h.values()\n",
Expand All @@ -1665,7 +1684,17 @@
" m2.limits = (0, None)\n",
" m2.migrad()\n",
"\n",
" return ntot, m1.valid, m1.values[0], m1.values[1], m1.fval, m2.valid, m2.values[0], m2.fval\n",
" return (\n",
" ntot,\n",
" m1.valid,\n",
" m1.values[0],\n",
" m1.values[1],\n",
" m1.fval,\n",
" m2.valid,\n",
" m2.values[0],\n",
" m2.fval,\n",
" )\n",
"\n",
"\n",
"result = Parallel(n_jobs=8)(delayed(run)(seed) for seed in range(1000))\n",
"ntot, valid1, ntot1, lambd1, minval1, valid2, lambd2, minval2 = np.transpose(result)"
Expand Down Expand Up @@ -1700,13 +1729,25 @@
"source": [
"plt.figure()\n",
"plt.title(\"$\\\\lambda$ estimate\")\n",
"plt.hist(lambd1, bins=20, alpha=0.5, label=f\"ExtendedBinnedNLL\\nvalid={np.mean(valid1)*100:.0f}%\\nmean = {np.mean(lambd1):.2f}\\nstd.dev. = {np.std(lambd1):.2f}\")\n",
"plt.hist(lambd2, bins=20, alpha=0.5, label=f\"BinnedNLL\\nvalid={np.mean(valid2)*100:.0f}%\\nmean = {np.mean(lambd2):.2f}\\nstd.dev. = {np.std(lambd2):.2f}\")\n",
"plt.hist(\n",
" lambd1,\n",
" bins=20,\n",
" alpha=0.5,\n",
" label=f\"ExtendedBinnedNLL\\nvalid={np.mean(valid1)*100:.0f}%\\nmean = {np.mean(lambd1):.2f}\\nstd.dev. = {np.std(lambd1):.2f}\",\n",
")\n",
"plt.hist(\n",
" lambd2,\n",
" bins=20,\n",
" alpha=0.5,\n",
" label=f\"BinnedNLL\\nvalid={np.mean(valid2)*100:.0f}%\\nmean = {np.mean(lambd2):.2f}\\nstd.dev. = {np.std(lambd2):.2f}\",\n",
")\n",
"plt.legend()\n",
"\n",
"plt.figure()\n",
"plt.hist(ntot1 / ntot - 1)\n",
"plt.title(f\"ExtendedBinnedNLL: $n_\\\\mathrm{{tot}}$ estimate\\n(mean-truth)/truth = {np.mean(ntot1) / np.mean(ntot) - 1:.3f}\")\n",
"plt.title(\n",
" f\"ExtendedBinnedNLL: $n_\\\\mathrm{{tot}}$ estimate\\n(mean-truth)/truth = {np.mean(ntot1) / np.mean(ntot) - 1:.3f}\"\n",
")\n",
"plt.xlabel(\"(estimate - truth) / truth\")\n",
"plt.xlim(-0.2, 0.2);"
]
Expand Down Expand Up @@ -1748,14 +1789,18 @@
"plt.sca(ax[0])\n",
"plt.hist(minval1, bins=50, density=True)\n",
"x = np.linspace(0, 100)\n",
"plt.plot(x, chi2(bins-2).pdf(x))\n",
"plt.title(f\"ExtendedBinnedNLL minimum value\\nndf = {bins-2} mean = {np.mean(minval1):.2f} median = {np.median(minval1):.2f}\")\n",
"plt.plot(x, chi2(bins - 2).pdf(x))\n",
"plt.title(\n",
" f\"ExtendedBinnedNLL minimum value\\nndf = {bins-2} mean = {np.mean(minval1):.2f} median = {np.median(minval1):.2f}\"\n",
")\n",
"\n",
"plt.sca(ax[1])\n",
"plt.hist(minval2, bins=50, density=True)\n",
"x = np.linspace(0, 100)\n",
"plt.plot(x, chi2(bins-2).pdf(x))\n",
"plt.title(f\"BinnedNLL minimum value\\nndf = {bins-2} mean = {np.mean(minval2):.2f} median = {np.median(minval2):.2f}\");"
"plt.plot(x, chi2(bins - 2).pdf(x))\n",
"plt.title(\n",
" f\"BinnedNLL minimum value\\nndf = {bins-2} mean = {np.mean(minval2):.2f} median = {np.median(minval2):.2f}\"\n",
");"
]
},
{
Expand Down Expand Up @@ -1784,7 +1829,9 @@
" x = expon.rvs(size=rng.poisson(npoints), random_state=rng)\n",
" w = rng.normal(0.1, 0.01, size=len(x))\n",
"\n",
" h = bh.Histogram(bh.axis.Regular(bins, np.min(x), np.max(x)), storage=bh.storage.Weight())\n",
" h = bh.Histogram(\n",
" bh.axis.Regular(bins, np.min(x), np.max(x)), storage=bh.storage.Weight()\n",
" )\n",
" h.fill(x, weight=w)\n",
" xe = h.axes[0].edges\n",
" n = h.values()\n",
Expand All @@ -1801,7 +1848,17 @@
" m2.limits = (0, None)\n",
" m2.migrad()\n",
"\n",
" return ntot, m1.valid, m1.values[0], m1.values[1], m1.fval, m2.valid, m2.values[0], m2.fval\n",
" return (\n",
" ntot,\n",
" m1.valid,\n",
" m1.values[0],\n",
" m1.values[1],\n",
" m1.fval,\n",
" m2.valid,\n",
" m2.values[0],\n",
" m2.fval,\n",
" )\n",
"\n",
"\n",
"result = Parallel(n_jobs=8)(delayed(run)(seed) for seed in range(1000))\n",
"ntot, valid1, ntot1, lambd1, minval1, valid2, lambd2, minval2 = np.transpose(result)"
Expand Down Expand Up @@ -1836,13 +1893,25 @@
"source": [
"plt.figure()\n",
"plt.title(\"$\\\\lambda$ estimate\")\n",
"plt.hist(lambd1, bins=20, alpha=0.5, label=f\"ExtendedBinnedNLL\\nvalid={np.mean(valid1)*100:.0f}%\\nmean = {np.mean(lambd1):.2f}\\nstd.dev. = {np.std(lambd1):.3f}\")\n",
"plt.hist(lambd2, bins=20, alpha=0.5, label=f\"BinnedNLL\\nvalid={np.mean(valid2)*100:.0f}%\\nmean = {np.mean(lambd2):.2f}\\nstd.dev. = {np.std(lambd2):.3f}\")\n",
"plt.hist(\n",
" lambd1,\n",
" bins=20,\n",
" alpha=0.5,\n",
" label=f\"ExtendedBinnedNLL\\nvalid={np.mean(valid1)*100:.0f}%\\nmean = {np.mean(lambd1):.2f}\\nstd.dev. = {np.std(lambd1):.3f}\",\n",
")\n",
"plt.hist(\n",
" lambd2,\n",
" bins=20,\n",
" alpha=0.5,\n",
" label=f\"BinnedNLL\\nvalid={np.mean(valid2)*100:.0f}%\\nmean = {np.mean(lambd2):.2f}\\nstd.dev. = {np.std(lambd2):.3f}\",\n",
")\n",
"plt.legend()\n",
"\n",
"plt.figure()\n",
"plt.hist(ntot1 / ntot - 1, bins=20)\n",
"plt.title(f\"ExtendedBinnedNLL: $n_\\\\mathrm{{tot}}$ estimate\\n(mean-truth)/truth = {np.mean(ntot1) / np.mean(ntot) - 1:.3f}\")\n",
"plt.title(\n",
" f\"ExtendedBinnedNLL: $n_\\\\mathrm{{tot}}$ estimate\\n(mean-truth)/truth = {np.mean(ntot1) / np.mean(ntot) - 1:.3f}\"\n",
")\n",
"plt.xlabel(\"(estimate - truth) / truth\");"
]
},
Expand Down Expand Up @@ -1879,14 +1948,18 @@
"plt.sca(ax[0])\n",
"plt.hist(minval1, bins=50, density=True)\n",
"x = np.linspace(0, 100)\n",
"plt.plot(x, chi2(bins-2).pdf(x))\n",
"plt.title(f\"ExtendedBinnedNLL minimum value\\nndf = {bins-2} mean = {np.mean(minval1):.2f} median = {np.median(minval1):.2f}\")\n",
"plt.plot(x, chi2(bins - 2).pdf(x))\n",
"plt.title(\n",
" f\"ExtendedBinnedNLL minimum value\\nndf = {bins-2} mean = {np.mean(minval1):.2f} median = {np.median(minval1):.2f}\"\n",
")\n",
"\n",
"plt.sca(ax[1])\n",
"plt.hist(minval2, bins=50, density=True)\n",
"x = np.linspace(0, 100)\n",
"plt.plot(x, chi2(bins-2).pdf(x))\n",
"plt.title(f\"BinnedNLL minimum value\\nndf = {bins-2} mean = {np.mean(minval2):.2f} median = {np.median(minval2):.2f}\");"
"plt.plot(x, chi2(bins - 2).pdf(x))\n",
"plt.title(\n",
" f\"BinnedNLL minimum value\\nndf = {bins-2} mean = {np.mean(minval2):.2f} median = {np.median(minval2):.2f}\"\n",
");"
]
},
{
Expand Down

0 comments on commit 22787da

Please sign in to comment.