Skip to content

Commit

Permalink
Add numerical calculations.
Browse files Browse the repository at this point in the history
  • Loading branch information
dp-rice committed Mar 12, 2024
1 parent 6745feb commit 6be286e
Show file tree
Hide file tree
Showing 188 changed files with 392 additions and 852 deletions.

Large diffs are not rendered by default.

Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 2 additions & 2 deletions docs/index.html
Original file line number Diff line number Diff line change
Expand Up @@ -145,9 +145,9 @@ <h1 class="title">Dan’s NAO Notebook</h1>

<div class="quarto-listing quarto-listing-container-default" id="listing-listing">
<div class="list quarto-listing-default">
<div class="quarto-post image-right" data-index="0" data-listing-date-sort="1709528400000" data-listing-file-modified-sort="1709652000850" data-listing-date-modified-sort="NaN" data-listing-reading-time-sort="18">
<div class="quarto-post image-right" data-index="0" data-listing-date-sort="1709528400000" data-listing-file-modified-sort="1710270528186" data-listing-date-modified-sort="NaN" data-listing-reading-time-sort="17">
<div class="thumbnail">
<p><a href="./notebooks/2024-02-22_StochasticMVP.html"> <p class="card-img-top"><img src="notebooks/2024-02-22_StochasticMVP_files/figure-html/cell-2-output-1.png" class="thumbnail-image card-img"/></p> </a></p>
<p><a href="./notebooks/2024-02-22_StochasticMVP.html"> <p class="card-img-top"><img src="notebooks/2024-02-22_StochasticMVP_files/figure-html/cell-5-output-1.png" class="thumbnail-image card-img"/></p> </a></p>
</div>
<div class="body">
<a href="./notebooks/2024-02-22_StochasticMVP.html">
Expand Down
743 changes: 219 additions & 524 deletions docs/notebooks/2024-02-22_StochasticMVP.html

Large diffs are not rendered by default.

Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
Diff not rendered.
9 changes: 8 additions & 1 deletion docs/search.json
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,7 @@
"href": "notebooks/2024-02-22_StochasticMVP.html#numerical-illustrations-to-do",
"title": "NAO Cost Estimate MVP – Adding noise",
"section": "Numerical illustrations (to-do)",
"text": "Numerical illustrations (to-do)"
"text": "Numerical illustrations (to-do)\n\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom numpy.polynomial.polynomial import Polynomial\nfrom numpy.polynomial.hermite_e import HermiteE\nfrom scipy.special import factorial\nfrom scipy.stats import poisson, norm, gamma\n\n\ndef cgf_x_series(mu: float, nu: float, order: int = 4) -&gt; Polynomial:\n coeffs = np.zeros(order + 1)\n # No zeroth coefficient\n k = np.arange(1, order + 1)\n coeffs[1:] = nu * (mu / nu) ** k / k**2\n return Polynomial(coeffs)\n\n\ndef expm1_series(order: int = 4) -&gt; Polynomial:\n k = np.arange(order + 1)\n return Polynomial(1.0 / factorial(k)) - 1\n\n\ndef cgf_y_series(mu: float, nu: float, order: int = 4):\n return cgf_x_series(mu, nu, order)(expm1_series(order)).cutdeg(order)\n\n\ndef cumulant_from_cgf(cgf: Polynomial, order: int) -&gt; float:\n return cgf.deriv(order)(0)\n\nSpot-check cumulants:\n\nmu = 2.0\nnu = 10.0\ncgf = cgf_x_series(mu, nu)(expm1_series()).cutdeg(4)\ncoeffs = [0, 1, (1, 1 / 2), (1, 3 / 2, 2 / 3), (1, 7 / 2, 4, 3 / 2)]\npredicted_cumulants = [mu * Polynomial(c)(mu / nu) for c in coeffs]\nfor k in range(5):\n print(cumulant_from_cgf(cgf, k), predicted_cumulants[k])\n\n0.0 0.0\n2.0 2.0\n2.2 2.2\n2.6533333333333333 2.6533333333333333\n3.7439999999999998 3.744\n\n\nCheck variance:\n\nnu = 10.0\nmu = np.arange(1, 100)\nk2 = [cumulant_from_cgf(cgf_y_series(m, nu), 2) for m in mu]\nplt.plot(mu, k2)\nplt.plot(mu, mu + (1 / 2) * mu**2 / nu, \":\")\nplt.plot(mu, (1 / 2) * mu**2 / nu, \"--\")\n\n\n\n\nNote that it is growing quadratically, but the \\(\\mu\\) term is not negligible.\n\ndef cornish_fisher(*cumulants):\n # cumulants = (k_1, k_2, ...)\n order = len(cumulants)\n if order &lt; 2:\n raise ValueError(\"Order of approximation must be &gt;= 2\")\n if order &gt; 4:\n raise ValueError(\"Order of approximation must be &lt;= 4\")\n sigma = np.sqrt(cumulants[1])\n poly = HermiteE((0, 1))\n if order &gt;= 3:\n gamma_1 = cumulants[2] / sigma**3\n h_1 = HermiteE((0, 0, 1)) / 6\n poly += gamma_1 * h_1\n if order &gt;= 4:\n gamma_2 = cumulants[3] / sigma**4\n h_2 = HermiteE((0, 0, 0, 1)) / 24\n h_11 = -HermiteE((0, 1, 0, 2)) / 36\n poly += gamma_2 * h_2 + gamma_1**2 * h_11\n return cumulants[0] + sigma * poly\n\n\nChecking the Cornish-Fisher expansion against common distributions\nCheck against Poisson distribution:\n\norder = 4\np = np.arange(0.01, 1.0, 0.01)\nx = norm.ppf(p)\n\nfor lamb in [1, 2, 4, 8]:\n poisson_cumulants = [lamb] * order\n for o in range(2, order + 1):\n cf_poisson = cornish_fisher(*poisson_cumulants[:o])\n plt.plot(p, cf_poisson(x), label=o)\n plt.plot(p, poisson(lamb).ppf(p), color=\"k\", linestyle=\"--\", label=\"exact\")\n plt.legend()\n plt.title(f\"$\\lambda$ = {lamb}\")\n plt.show()\n\n\n\n\n\n\n\n\n\n\n\n\n\nCheck against Gamma distribution:\n\nscale = 1\nk = np.arange(1, order + 1)\n\nfor shape in [1, 2, 4, 8]:\n gamma_cumulants = factorial(k - 1) * shape * scale**k\n for o in range(2, order + 1):\n cf_gamma = cornish_fisher(*gamma_cumulants[:o])\n plt.plot(p, cf_gamma(x), label=o)\n plt.plot(\n p, gamma(shape, scale=scale).ppf(p), color=\"k\", linestyle=\"--\", label=\"exact\"\n )\n plt.legend()\n plt.title(f\"shape = {shape}\")\n plt.show()\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nCornish-Fisher expansions for the cumulative count distribution\n\norder = 4\nnu = 2.0\nfor mu in [1, 2, 4, 8]:\n cgf = cgf_y_series(mu, nu)\n cumulants = [cumulant_from_cgf(cgf, k) for k in [1, 2, 3, 4]]\n for o in range(2, order + 1):\n cf = cornish_fisher(*cumulants[:o])\n plt.plot(p, cf(x), label=o)\n plt.legend()\n plt.title(f\"$\\mu$ = {mu}\")\n plt.show()\n\n\n\n\n\n\n\n\n\n\n\n\n\n\norder = 4\nnu = 4.0\nfor mu in [1, 2, 4, 8]:\n cgf = cgf_y_series(mu, nu)\n cumulants = [cumulant_from_cgf(cgf, k) for k in [1, 2, 3, 4]]\n for o in range(2, order + 1):\n cf = cornish_fisher(*cumulants[:o])\n plt.plot(p, cf(x), label=o)\n plt.legend()\n plt.title(f\"$\\mu$ = {mu}\")\n plt.show()\n\n\n\n\n\n\n\n\n\n\n\n\n\n\norder = 4\nnu = 10.0\nfor mu in [1, 10, 100, 1000]:\n cgf = cgf_y_series(mu, nu)\n cumulants = [cumulant_from_cgf(cgf, k) for k in [1, 2, 3, 4]]\n for o in range(2, order + 1):\n cf = cornish_fisher(*cumulants[:o])\n plt.plot(p, cf(x), label=o)\n plt.legend()\n plt.title(f\"$\\mu$ = {mu}\")\n plt.show()"
},
{
"objectID": "notebooks/2024-02-22_StochasticMVP.html#implications-for-cost",
Expand Down Expand Up @@ -726,5 +726,12 @@
"title": "2023-10-10 Analyze prefilter experiment qPCR",
"section": "Amplification curves",
"text": "Amplification curves\n\namp_data |&gt;\n filter(Task == \"UNKNOWN\") |&gt;\n ggplot(mapping = aes(\n x = `Cycle Number`,\n y = dRn,\n color = as.factor(group),\n group = Well\n )) +\n geom_line() +\n geom_line(mapping = aes(\n x = `Cycle Number`,\n y = Threshold\n ), color = \"Grey\") +\n scale_y_log10() +\n facet_wrap(~Target, scales = \"free\")\n\nWarning in self$trans$transform(x): NaNs produced\n\n\nWarning: Transformation introduced infinite values in continuous y-axis\n\n\nWarning: Removed 49 rows containing missing values (`geom_line()`).\n\n\n\n\n\n\namp_data |&gt;\n filter(Task == \"NTC\") |&gt;\n ggplot(mapping = aes(\n x = `Cycle Number`,\n y = dRn,\n group = Well\n )) +\n geom_line() +\n geom_line(mapping = aes(\n x = `Cycle Number`,\n y = Threshold\n ), color = \"Grey\") +\n scale_y_log10() +\n facet_wrap(~Target, scales = \"free\")\n\nWarning in self$trans$transform(x): NaNs produced\n\n\nWarning: Transformation introduced infinite values in continuous y-axis\n\n\nWarning: Removed 149 rows containing missing values (`geom_line()`).\n\n\n\n\n\n\nplot_amp &lt;- function(data, color) {\n ggplot(data, aes(x = `Cycle Number`, y = dRn)) +\n geom_line(mapping = aes(\n color = as.factor({{ color }}),\n group = Well,\n )) +\n scale_y_log10(limits = c(1e-3, 1e1))\n}\n\nruler &lt;- function(y0_from, num_rules) {\n y0 &lt;- 10^seq(from = y0_from, by = -1, length.out = num_rules)\n rules &lt;- crossing(`Cycle Number` = amp_data$`Cycle Number`, y0 = y0) |&gt;\n mutate(dRn = y0 * 2^`Cycle Number`)\n geom_line(\n data = rules,\n mapping = aes(group = y0),\n color = \"black\"\n )\n}\n\nplot_amp_with_ruler &lt;- function(target, y0_from, num_rules) {\n amp_data |&gt;\n filter(!is.na(quantity), Target == target) |&gt;\n plot_amp(quantity) +\n ruler(y0_from, num_rules) +\n geom_line(mapping = aes(\n x = `Cycle Number`,\n y = Threshold\n ), color = \"Grey\") +\n labs(title = target)\n}\n\n\nplot_amp_with_ruler(\"16S\", -5.5, 5)\n\nWarning in self$trans$transform(x): NaNs produced\n\n\nWarning: Transformation introduced infinite values in continuous y-axis\n\n\nWarning: Removed 51 rows containing missing values (`geom_line()`).\n\n\nWarning: Removed 134 rows containing missing values (`geom_line()`).\n\n\n\n\nplot_amp_with_ruler(\"Cov2\", -6.5, 5)\n\nWarning in self$trans$transform(x): NaNs produced\n\n\nWarning: Transformation introduced infinite values in continuous y-axis\n\n\nWarning: Removed 133 rows containing missing values (`geom_line()`).\n\n\n\n\nplot_amp_with_ruler(\"CrA\", -4.5, 5)\n\nWarning in self$trans$transform(x): NaNs produced\n\n\nWarning: Transformation introduced infinite values in continuous y-axis\n\n\nWarning: Removed 34 rows containing missing values (`geom_line()`).\n\n\nWarning: Removed 133 rows containing missing values (`geom_line()`).\n\n\n\n\nplot_amp_with_ruler(\"PMMV\", -8, 5)\n\nWarning in self$trans$transform(x): NaNs produced\n\n\nWarning: Transformation introduced infinite values in continuous y-axis\n\n\nWarning: Removed 51 rows containing missing values (`geom_line()`).\n\n\nWarning: Removed 136 rows containing missing values (`geom_line()`)."
},
{
"objectID": "notebooks/2024-02-22_StochasticMVP.html#numerical-calculations",
"href": "notebooks/2024-02-22_StochasticMVP.html#numerical-calculations",
"title": "NAO Cost Estimate MVP – Adding noise",
"section": "Numerical calculations",
"text": "Numerical calculations\n\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom numpy.polynomial.polynomial import Polynomial\nfrom numpy.polynomial.hermite_e import HermiteE\nfrom scipy.special import factorial\nfrom scipy.stats import poisson, norm, gamma\n\n\ndef cgf_x_series(mu: float, nu: float, order: int = 4) -&gt; Polynomial:\n coeffs = np.zeros(order + 1)\n # No zeroth coefficient\n k = np.arange(1, order + 1)\n coeffs[1:] = nu * (mu / nu) ** k / k**2\n return Polynomial(coeffs)\n\n\ndef expm1_series(order: int = 4) -&gt; Polynomial:\n k = np.arange(order + 1)\n return Polynomial(1.0 / factorial(k)) - 1\n\n\ndef cgf_y_series(mu: float, nu: float, order: int = 4):\n return cgf_x_series(mu, nu, order)(expm1_series(order)).cutdeg(order)\n\n\ndef cumulant_from_cgf(cgf: Polynomial, order: int) -&gt; float:\n return cgf.deriv(order)(0)\n\nSpot-check cumulants:\n\nmu = 2.0\nnu = 10.0\ncgf = cgf_x_series(mu, nu)(expm1_series()).cutdeg(4)\ncoeffs = [0, 1, (1, 1 / 2), (1, 3 / 2, 2 / 3), (1, 7 / 2, 4, 3 / 2)]\npredicted_cumulants = [mu * Polynomial(c)(mu / nu) for c in coeffs]\nfor k in range(5):\n print(cumulant_from_cgf(cgf, k), predicted_cumulants[k])\n\n0.0 0.0\n2.0 2.0\n2.2 2.2\n2.6533333333333333 2.6533333333333333\n3.7439999999999998 3.744\n\n\nCheck variance:\n\nnu = 10.0\nmu = np.arange(1, 100)\nk2 = [cumulant_from_cgf(cgf_y_series(m, nu), 2) for m in mu]\nplt.plot(mu, k2)\nplt.plot(mu, mu + (1 / 2) * mu**2 / nu, \":\")\nplt.plot(mu, (1 / 2) * mu**2 / nu, \"--\")\n\n\n\n\nNote that it is growing quadratically, but the \\(\\mu\\) term is not negligible.\n\ndef cornish_fisher(*cumulants):\n # cumulants = (k_1, k_2, ...)\n order = len(cumulants)\n if order &lt; 2:\n raise ValueError(\"Order of approximation must be &gt;= 2\")\n if order &gt; 4:\n raise ValueError(\"Order of approximation must be &lt;= 4\")\n sigma = np.sqrt(cumulants[1])\n poly = HermiteE((0, 1))\n if order &gt;= 3:\n gamma_1 = cumulants[2] / sigma**3\n h_1 = HermiteE((0, 0, 1)) / 6\n poly += gamma_1 * h_1\n if order &gt;= 4:\n gamma_2 = cumulants[3] / sigma**4\n h_2 = HermiteE((0, 0, 0, 1)) / 24\n h_11 = -HermiteE((0, 1, 0, 2)) / 36\n poly += gamma_2 * h_2 + gamma_1**2 * h_11\n return cumulants[0] + sigma * poly\n\n\nChecking the Cornish-Fisher expansion against common distributions\nCheck against Poisson distribution:\n\norder = 4\np = np.arange(0.01, 1.0, 0.01)\nx = norm.ppf(p)\n\nfor lamb in [1, 2, 4, 8]:\n poisson_cumulants = [lamb] * order\n for o in range(2, order + 1):\n cf_poisson = cornish_fisher(*poisson_cumulants[:o])\n plt.plot(p, cf_poisson(x), label=o)\n plt.plot(p, poisson(lamb).ppf(p), color=\"k\", linestyle=\"--\", label=\"exact\")\n plt.legend()\n plt.title(f\"$\\lambda$ = {lamb}\")\n plt.show()\n\n\n\n\n\n\n\n\n\n\n\n\n\nWhen \\(\\lambda &gt; 1\\), you quickly converge to a Gaussian distribution, and the higher-order corrections don’t matter.\nCheck against Gamma distribution:\n\nscale = 1\nk = np.arange(1, order + 1)\n\nfor shape in [1, 2, 4, 8]:\n gamma_cumulants = factorial(k - 1) * shape * scale**k\n for o in range(2, order + 1):\n cf_gamma = cornish_fisher(*gamma_cumulants[:o])\n plt.plot(p, cf_gamma(x), label=o)\n plt.plot(\n p, gamma(shape, scale=scale).ppf(p), color=\"k\", linestyle=\"--\", label=\"exact\"\n )\n plt.legend()\n plt.title(f\"shape = {shape}\")\n plt.show()\n\n\n\n\n\n\n\n\n\n\n\n\n\nThe shape parameter controls deviations from Gaussian. For shape &gt; 1, three terms gets you close and 4 gets you very close.\n\n\nCornish-Fisher expansions for the cumulative count distribution\nFor the cumulative count distribution, there are two parameters, \\(\\mu\\) (the mean) and \\(\\nu\\) (which determines the shape).\nFor small \\(\\nu\\), the noise is quickly dominated by the latent variable so the distribution goes to a constant shape, scaled by \\(\\mu\\).\nEven with \\(\\nu = 2\\), the CF expansion converges quickly. 3 terms is quite good, 4 is indistinguisable from higher:\n\norder = 6\nnu = 2.0\nfor mu in [1, 2, 4, 8, 32]:\n cgf = cgf_y_series(mu, nu)\n cumulants = [cumulant_from_cgf(cgf, k) for k in [1, 2, 3, 4]]\n for o in range(2, order + 1):\n cf = cornish_fisher(*cumulants[:o])\n plt.plot(p, cf(x), label=o)\n plt.legend()\n plt.title(r\"$\\nu$\" f\" = {nu}, \" r\"$\\mu$ = \" f\"{mu}\")\n plt.show()\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nWith larger \\(\\nu\\), the shape of the distribution is closer to Gaussian:\n\norder = 6\nnu = 10.0\nfor mu in [1, 10, 100, 1000]:\n cgf = cgf_y_series(mu, nu)\n cumulants = [cumulant_from_cgf(cgf, k) for k in [1, 2, 3, 4]]\n for o in range(2, order + 1):\n cf = cornish_fisher(*cumulants[:o])\n plt.plot(p, cf(x), label=o)\n plt.legend()\n plt.title(r\"$\\nu$\" f\" = {nu}, \" r\"$\\mu$ = \" f\"{mu}\")\n plt.show()\n\n\n\n\n\n\n\n\n\n\n\n\n\nWith \\(\\nu = 100\\), there is a wide Poisson-dominated regime the convergence to Gaussian is quite fast.\n\norder = 4\nnu = 100.0\nfor mu in [1, 10, 100, 1000]:\n cgf = cgf_y_series(mu, nu)\n cumulants = [cumulant_from_cgf(cgf, k) for k in [1, 2, 3, 4]]\n for o in range(2, order + 1):\n cf = cornish_fisher(*cumulants[:o])\n plt.plot(p, cf(x), label=o)\n plt.legend()\n plt.title(r\"$\\nu$\" f\" = {nu}, \" r\"$\\mu$ = \" f\"{mu}\")\n plt.show()\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nPercentiles\nThe following show the mean (dashed grey lines) and the 10th percentile (black dots) of cumulative counts as a function of the mean \\(\\mu\\)” Each plot has a different shape parameter \\(\\nu\\).\n\norder = 4\n# 10th percentile\np = 0.1\nx = norm.ppf(p)\n\nmu = np.arange(1, 100)\nfor nu in [4, 10, 100]:\n cgfs = [cgf_y_series(m, nu) for m in mu]\n cumulants = [[cumulant_from_cgf(cgf, k) for k in [1, 2, 3, 4]] for cgf in cgfs]\n for m, cs in zip(mu, cumulants):\n cf = cornish_fisher(*cs)\n plt.plot(m, cf(x), \".k\")\n plt.plot(mu, mu, \"--\", color=\"grey\")\n plt.xlabel(r\"$\\mu$\")\n plt.ylabel(\"Cumulative counts\")\n plt.title(r\"$\\nu = $\" f\"{nu}\")\n plt.show()\n\n\n\n\n\n\n\n\n\n\nObservations:\n\nSmaller \\(\\nu\\) means that we have a greater reduction of the 10th percentile from the mean.\nSmaller \\(\\nu\\) means that the 10th percentile increase linearly (as expected in the latent regime).\nAs \\(\\nu\\) gets larger, there is a wider Poisson-dominated regime where the tenth percentile increases like \\(\\mu + \\mu^{1/2} const.\\)."
}
]
Loading

0 comments on commit 6be286e

Please sign in to comment.