{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Asymmetric Uncertainties\n",
    "\n",
    "A model for asymmetric uncertainties [(Eq 21 in Barlow, 2003)](https://www.slac.stanford.edu/econf/C030908/papers/WEMT002.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import PchipInterpolator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logL(x, mu, sigma_lo, sigma_hi):\n",
    "    \"\"\"Log-likelihood for a 1D confidence interval with asymmetric errors.\n",
    "\n",
    "    Based on eq. 21 from R. Barlow, \"Asymmetric Errors,\" PHYSTAT 2003.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : float or ndarray\n",
    "        Value(s) to evaluate the likelihood.\n",
    "    mu : float\n",
    "        Maximum likelihood value for the interval.\n",
    "    sigma_lo : float\n",
    "        Lower edge of central 68% confidence interval.\n",
    "    sigma_hi : float\n",
    "        Upper edge of central 68% confidence interval.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    logL : float or ndarray\n",
    "        Log likelihood.\n",
    "    \"\"\"\n",
    "    beta = sigma_hi / sigma_lo\n",
    "    gamma = (sigma_hi*sigma_lo) / (sigma_hi - sigma_lo)\n",
    "    cut = (x - mu)/gamma > -1\n",
    "    return -0.5 * (np.log(1. + (x - mu)/gamma, where=cut, out=np.inf*np.ones_like(x)) / np.log(beta))**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Log Likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize = (3,5))\n",
    "m = 33.82 # central value, an example\n",
    "a = np.linspace(25, 40, 1000)\n",
    "# Negative and Positive uncertainties\n",
    "sig_pos = 0.78\n",
    "sig_neg = 0.76\n",
    "\n",
    "lnL = logL(a, mu=m, sigma_lo=sig_neg, sigma_hi=sig_pos)\n",
    "#print( lnL)\n",
    "\n",
    "ax.plot(a, lnL)\n",
    "ax.set(xlabel = \"m\", ylabel = \"lnL\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probability Distribution Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(5,5))\n",
    "L = np.exp(lnL)\n",
    "L/=np.trapz(L, a)\n",
    "# print(L)\n",
    "ax.plot(a, L)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cumulative PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(5,5))\n",
    "cdf = np.zeros_like(a)\n",
    "for i in range(a.size):\n",
    "    cdf[i] = np.trapz(L[:i], a[:i])\n",
    "ax.plot(a, cdf)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random Sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(5,5))\n",
    "cut = cdf > 0 \n",
    "Cdf = PchipInterpolator(cdf[cut], a[cut])\n",
    "u = np.random.uniform(0.,1., 50000)\n",
    "A = Cdf(u)\n",
    "ax.hist(A, bins=np.linspace(25,40,100), alpha=0.7, density=True)\n",
    "ax.plot(a, L)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ValueCI(object):\n",
    "\n",
    "    def __init__(self, val, ci_lo, ci_hi):\n",
    "        \"\"\"Initialize a value with potentially asymmetric error bars.\n",
    "        Assume the CI refers to 68% region around the central value.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        val : float\n",
    "            Centra value.\n",
    "        ci_lo : float\n",
    "            Lower range of 68% C.I. around central value.\n",
    "        ci_hi : float\n",
    "            Upper range of 68% C.I. around central value.\n",
    "        \"\"\"\n",
    "        self.value = val\n",
    "        self.ci_lo = ci_lo\n",
    "        self.ci_hi = ci_hi\n",
    "    \n",
    "    def likelihood(self, x):\n",
    "        \"\"\"Likelihood for a 1D confidence interval with asymmetric errors.\n",
    "\n",
    "        Based on eq. 21 from R. Barlow, \"Asymmetric Errors,\" PHYSTAT 2003.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        x : float or ndarray\n",
    "            Value(s) to evaluate the likelihood.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        L : float or ndarray\n",
    "            Likelihood.\n",
    "        \"\"\"\n",
    "        if self.ci_lo == self.ci_hi:\n",
    "            pass\n",
    "        else:\n",
    "            beta = self.ci_hi / self.ci_lo\n",
    "            gamma = (self.ci_hi*self.ci_lo) / (self.ci_hi - self.ci_lo)\n",
    "            cut = (x - self.value)/gamma > -1\n",
    "            return np.exp(-0.5 * (np.log(1. + (x - self.value)/gamma, where=cut, out=np.inf*np.ones_like(x)) / np.log(beta))**2)\n",
    "\n",
    "    def get_pdf(self, a):\n",
    "        if self.ci_hi == self.ci_lo:\n",
    "            pass\n",
    "        else:\n",
    "            L = self.likelihood(a)/np.trapz(self.likelihood(a), a)\n",
    "            return L\n",
    "    \n",
    "    def get_cdf(self, a):\n",
    "        if self.ci_hi == self.ci_lo:\n",
    "            pass\n",
    "        else:\n",
    "            cdf = np.zeros_like(a)\n",
    "            for i in range(a.size):\n",
    "                cdf[i] = np.trapz(self.get_pdf(a)[:i], a[:i])\n",
    "            return cdf\n",
    "    \n",
    "    \n",
    "    def get_random(self, n=1):\n",
    "        \"\"\"Randomly sample values from the distribution.\n",
    "        If the distribution is asymmetric, treat it as a 2-sided Gaussian.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        n : int\n",
    "            Number of random draws.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        Returns n random draws from a symmetric/asymmetric Gaussian about\n",
    "        the central value.\n",
    "        \"\"\"\n",
    "        if self.ci_lo == self.ci_hi:\n",
    "            return np.random.normal(loc=self.value, scale=self.ci_lo, size=n)\n",
    "        else:\n",
    "#             return np.random.normal(loc=self.value, scale=self.ci_lo, size=n)\n",
    "            a = np.linspace(self.value - 5*self.ci_lo, self.value + 5*self.ci_hi, n)\n",
    "            L = self.get_pdf(a)\n",
    "            cdf = np.zeros_like(a)\n",
    "            for i in range(n):\n",
    "                cdf[i] = np.trapz(L[:i], a[:i])\n",
    "            cut = self.get_cdf(a) > 0 \n",
    "            Cdf = PchipInterpolator(self.get_cdf(a)[cut], a[cut])\n",
    "            u = np.random.uniform(0.,1., 50000)\n",
    "            A = Cdf(u)\n",
    "#             cut = self.get_cdf(a) > 0 \n",
    "#             Cdf = PchipInterpolator(self.get_cdf(a)[cut], a[cut])\n",
    "#             u = np.random.uniform(0.,1., 50000)\n",
    "            return A\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def asym(arr, val, sig_hi, sig_lo):\n",
    "#     beta = sigma_hi / sigma_lo\n",
    "#     gamma = (sigma_hi*sigma_lo) / (sigma_hi - sigma_lo)\n",
    "#     cut = (x - mu)/gamma > -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a = ValueCI(20, 3, 2)\n",
    "# x = np.linspace(0, 50, 100)\n",
    "# y = a.get_random(2000)\n",
    "# # print(y)\n",
    "# # ln = a.get_cdf(x)\n",
    "# # # ln/=np.trapz(ln, x)\n",
    "# # # print(ln)\n",
    "# fig, ax = plt.subplots(1, 1, figsize=(5,5))\n",
    "# ax.hist(y, bins=x, alpha=0.7, density=True)\n",
    "# ax.plot(x, a.get_pdf(x))\n",
    "# # ax.plot(x,y)\n",
    "# fig.tight_layout()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}