{"cells":[{"metadata":{},"cell_type":"markdown","source":"## 1. Kullback-Leibler divergence \n\n**Introduction.**\n\n`'KL'-divergence` (or `Kullback-Leibler` divergence) between 2 probability distributions of $P$ from $Q$ is defined as\n\n$$ d_{KL} \\left( P \\Vert Q \\right) = \\int \\log \\left( \\frac{dP}{dQ} \\right) d P $$\n\nIn the `discreate-form`, we can rewrite the above formula as\n\n$$ d_{KL} \\left( P \\Vert Q \\right) = \\sum_{k=1}^n \\log \\left( \\frac{p(x_k)}{q(x_k)} \\right) p(x_k). $$\n\nThe Kullback–Leibler divergence is `defined (exists)` only if for all ${\\displaystyle x},$ then $P$ is absolutely continuous with respect to $Q;$ that meant $Q(x) = 0$ implies $P(x) = 0$, or $P \\ll Q$. Moreover, $d_{KL} \\left( P \\Vert Q \\right) = 0$ if $P(x) = Q(x), \\; \\forall x.$ \n\nOtherwise,  \n\n$$ d_{KL} \\left( P \\Vert Q \\right) = \\infty, \\quad $$\nif $P$ is not absolutely continuous w.r.t. to $Q$.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport scipy.stats as sst\nimport seaborn as sns\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**`Example 1.1.`** `Approximate the KL-divergence of` $\\mathcal{N}(0, 1)$ from $\\mathcal{N}(\\mu, 1)$ for all $\\mu \\in \\mathbb{R}$\n\nDefine a function of the approximation for KL-divergence by using the availabel function on `scipy.stats`. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def kl_divergence(p, q):\n    return sst.entropy(p, q)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Set the partitions then display**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let the partitions = 201 points.\nn_pat = 201\na, b = -30, 30\nx = np.arange(a, b, 1/n_pat)\n\n## Fix p to the random-variables generated by the standard-normal distribution\np = sst.norm.pdf(x, 0, 1)\n\n## Initialize the KL-divergence of P from Q\ndPQ = []\n\n## Create the array of the mean of the normal distributions's sequence\nmuy = np.linspace(-8, 8, 17)\nfor mu in muy:\n    qk = sst.norm.pdf(x, mu, 1)\n    dPQ.append(kl_divergence(p, qk))\n\n## Set the size of figure    \nsns.set(rc = {\"figure.figsize\" : (15, 12)})\n\n## Show the d_KL on the top of the figure\nplt.subplot(211)    \nplt.plot(muy, dPQ)\n\n## Add the title and labels\nplt.title('KL divergences of $\\mathcal{N}(0,1)$ from $\\mathcal{N}(\\mu, 1)$')\nplt.ylabel('$d(P \\Vert Q)$')\nplt.xlabel('$\\mu$')\n\n## Show a few distributions in the normal distribution sequences N(mu_k, 1)\nplt.subplot(212)\nfor mu in muy[::5]:\n    plt.plot(x, sst.norm.pdf(x, mu, 1), label = '$\\mu = $'+str(mu))\n    plt.legend(loc = 'upper right')\n\n## Add the x,y_labels    \nplt.ylabel('$\\mathcal{N}(\\mu, 1)$')\nplt.xlabel('some normal distribution in the list')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**`Comment / intuition:`** The KL-divergences $d_{KL}(P \\Vert Q)$ looks like a parabol when the normal distributions take the same variance.\n\n**`Explain the intuition`**\n\nLet $\\phi(\\mu, \\sigma)$ be the density of $\\mathcal{N}(\\mu, \\sigma)$ and $\\Vert f - g \\Vert_{L^2}^2 \\triangleq \\displaystyle \\int (f(t) - g(t))^2 dt$, then\n\n$i) \\; \\displaystyle{ d_{KL}\\left( \\phi(0, 1) \\Vert \\phi(t, 1) \\right) = \\dfrac{t^2}{2}, \\quad \\forall t \\in \\mathbb{R} }$\n\n$ii) \\; \\displaystyle \\Vert \\phi(t, 1) - \\phi(0, 1) \\Vert_{L^2} \\leq \\dfrac{t^2}{2}, \\quad \\forall t \\in \\mathbb{R}$\n\n$iii) \\; \\displaystyle{ d_{KL}\\left( \\phi(\\mu_1, \\sigma_1) \\Vert \\phi(\\mu_2, \\sigma_2) \\right) = \\dfrac{\\sigma_1^2 - \\sigma_2^2 + (\\mu_1 - \\mu_2)^2}{2 \\sigma_2^2} + \\log \\dfrac{\\sigma_1}{\\sigma_2}, \\quad \\forall \\mu_1, \\mu_2 \\in \\mathbb{R} }$ and $\\forall \\sigma_1, \\sigma_2 > 0$","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Entropy, cross-entropy via KL-divergences\n\nWe have these formulas:\n\n#### Entropy of the information X\n$$ \\displaystyle H(p) = \\mathbb{E}(I(X)) = -\\sum_{i=1}^n p_i \\log(x_i),$$ \nwhere $p=(p_i)_i$ and $p_i = p(x_i)$ is the probability of that $x = x_i$, (remember that $p_i \\in [0, 1]$ and $\\sum_{i=1}^n p_i = 1).$\n\n### Cross-entropy\n$$ H(p, q) = H(p) + d_{KL}(p \\Vert q)$$\n\n**The available function `entropy` in `scipy.stats`**\n\nSee: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.entropy.html\n\n**Entropy example.** Imagine that we have an unfair coin which, when flipped, has a 99% chance of landing heads and only 1% chance of landing tails. \n\nThen p = [0.99, 0.01] and \n$$H(p) = -\\left[ 0.99 \\log(0.99) + 0.01 \\log(0.01) \\right] \\approx 0.056$$\n\nNow, let's see how it works on scipy.stats","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"p = [0.99, 0.01]\nHp = sst.entropy(p)\nprint('H(p) =\\t \\t' + str(Hp))\nprint('Check again: \\t'+ str(-(0.99*np.log(0.99) + 0.01*np.log(0.01))))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"One thing to note; if we are dealing with information expressed in `bits` (i.e. each `bit` is either 0 or 1) the logarithm has a base of 2, then\n$$H(p) = -\\left[ 0.99 \\log_2(0.99) + 0.01 \\log_2(0.01) \\right] \\approx 0.08$$\n\nSo, in `scipy.stats`; we must choose `base = 2` in `entropy(....)`. Let's check this again","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"Hp2 = sst.entropy(p, base = 2)\nprint('H_2(p) =\\t' + str(Hp2))\nprint('check again: \\t'+ str(-(0.99*np.log2(0.99) + 0.01*np.log2(0.01))))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Cross-entropy example.**\n\nIn the syntax `scipy.stats.entropy(p_k, q_k, base, axis)`, they conclude that if `qk != None` then compute the Kullback-Leibler divergence S = sum(pk * log(pk / qk), axis=axis); this meant the `cross-entropy` will be define as\n\n`H(p, q) = entropy(p_k) + entropy(p_k, q_k)`","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## values of true distribution\np = [0.1, 0.2, 0.2, 0.3, 0.2]\n## values of the prediction\nq = [0.12, 0.28, 0.25, 0.2 ,0.15]\n## cross-entropy\nsst.entropy(p, q) + sst.entropy(p)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Another examples of KL-divergences. :D","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## KL divergences between normal & exponential distribution\nx = np.arange(0, 6, 0.01)\np = sst.norm.pdf(x, 0, 1)\nq = sst.expon.pdf(x, 1)\nsst.entropy(p, q), sst.entropy(q, p)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Explain.** Let $p = \\mathcal{N}(0, 1)$ and $q = \\mathcal{E}(1)$, then\n$$d_{KL}(q \\Vert p) = \\sqrt{2\\pi} \\int_0^{\\infty} \\left( \\frac{x^2 - 2x}{2} \\right) e^{-x} dx \\approx 1.608$$\nis well-defined, while\n$$d_{KL}(p \\Vert q) = \\infty $$\nsince $p$ is not absolutely continuous w.r.t to $q.$ Indeed,\n$$q(x) = \\left \\lbrace \\begin{array}{ll} e^{-x} & \\text{ if } x \\geq 0 \\\\ 0 & \\text{ otherwise } \\end{array} \\right. \\text{ and } p(x) = \\frac{1}{\\sqrt{2\\pi}} e^{-x^2/2},$$\ntake $A = (-1, 0) \\in \\mathcal{B}(\\mathbb{R})$ then\n$$Q(A) = \\int_A q(x) dx = 0$$\nbut\n$$P(A) =\\int_A p(x) dx \\approx 0.3413 \\neq 0$$","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## Calculate P(A) by using the c.d.f of N(0, 1) on (-1, 0)\nsst.norm.cdf(0) - sst.norm.cdf(-1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`KL divergence between 2 Bernoulli distributions` $B(1, p)$ `and` $B(1, q).$\n\nIt is easy to verify that\n$$ d_{KL}(p \\Vert q) = p \\log\\left( \\frac{p}{q} \\right) + (1 - p) \\log \\left( \\frac{1-p}{1-q} \\right) $$\nwhile\n$$ d_{KL}(q \\Vert p) = q \\log\\left( \\frac{q}{p} \\right) + (1 - q) \\log \\left( \\frac{1-q}{1-p} \\right) $$\n\nIn the following example, we take $p \\sim B(1, 0.05)$ and $q \\sim B(1, 0.99)$","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"p, q = 0.05, 0.99\nP = sst.bernoulli.cdf(k = [0, 1], loc = 0, p = p)\nQ = sst.bernoulli.cdf(k = [0, 1], loc = 0, p = q)\nprint(sst.entropy(P, Q), sst.entropy(Q, P))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Summary\n\n- 1) The `KL` divergences is not a distance since it is not symmetric.\n- 2) But $d_{KL}(p \\Vert q) \\geq 0$, indeed\n$$ d_{KL}(p \\Vert q) = \\int_{ \\lbrace pq > 0 \\rbrace}  p \\log \\left( \\frac{p}{q} \\right) = - 2\\int_{ \\lbrace pq > 0 \\rbrace}  p \\log \\left( \\sqrt{ \\frac{q}{p} } - 1 + 1 \\right) \\geq -2 \\int_{ \\lbrace pq > 0 \\rbrace}  p \\left( \\sqrt{ \\frac{q}{p} } - 1 \\right) = \\int_{ \\lbrace pq > 0 \\rbrace} \\left( \\sqrt{p} - \\sqrt{q} \\right)^2 $$\n- 3) The `JS` (or Jensen-Shannon) divergence, is the symmetrized version of $d_{KL}$ defined as\n$$ d_{JS}(p, q) = \\frac{1}{2} \\left( d_{KL}\\left(p \\left \\Vert \\frac{p+q}{2} \\right. \\right) + d_{KL}\\left( q \\left \\Vert \\frac{p+q}{2} \\right. \\right) \\right) $$","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## 2. Jensen-Shannon divergence\n\nUp to the section `KL divergence` has been refered in the previous section!\n\nThe `JS-`divergence between 2 probability distribution $P, Q$ defined as\n$$ d_{JS}(p, q) = \\frac{1}{2} \\left( d_{KL}\\left(p \\left \\Vert \\frac{p+q}{2} \\right. \\right) + d_{KL}\\left( q \\left \\Vert \\frac{p+q}{2} \\right. \\right) \\right) $$\n\nThe `JS-`divergence between 2 probability distribution $P, Q$ defined as\n$$ d_{JS}(p, q) = \\frac{1}{2} \\left( d_{KL}\\left(p \\left \\Vert \\frac{p+q}{2} \\right. \\right) + d_{KL}\\left( q \\left \\Vert \\frac{p+q}{2} \\right. \\right) \\right) $$","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## Define a function\ndef JS_divergence(p, q):\n    return ( sst.entropy(p, (p+q)/2 ) + sst.entropy(q, (p+q)/2 ))/2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let the partitions = 201 points.\nn_pat = 201\na, b = -30, 30\nx = np.arange(a, b, 1/n_pat)\n\n## Fix p to the random-variables generated by the standard-normal distribution\np = sst.norm.pdf(x, 0, 1)\n\n## Initialize the KL-divergence of P from Q\ndPQ = []\n\n## Create the array of the mean of the normal distributions's sequence\nmuy = np.linspace(-10, 10, 21)\nfor mu in muy:\n    qk = sst.norm.pdf(x, mu, 1)\n    dPQ.append(JS_divergence(p, qk))\n\n## Set the size of figure    \nsns.set(rc = {\"figure.figsize\" : (15, 12)})\n\n## Show the d_JS on the top of the figure\nplt.subplot(211)    \nplt.plot(muy, dPQ)\n\n## Add the title and labels\nplt.title('JS divergences between $\\mathcal{N}(0,1)$ and $\\mathcal{N}(\\mu, 1)$')\nplt.ylabel('$d_{JS}(p \\Vert q)$')\nplt.xlabel('$\\mu$')\n\n## Show a few distributions in the normal distribution sequences N(mu_k, 1)\nplt.subplot(212)\nfor mu in muy[::6]:\n    plt.plot(x, sst.norm.pdf(x, mu, 1), label = '$\\mu = $'+str(mu))\n    plt.legend(loc = 'upper right')\n\n## Add the x,y_labels    \nplt.ylabel('$\\mathcal{N}(\\mu, 1)$')\nplt.xlabel('some normal distribution in the list')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Comment.** \n- 1) When the value of mean $\\mu = 0,$ that mean $q(x) = \\mathcal{N}(0, 1) = p(x)$ and hence $d_{JS}(p, q) = 0$. \n\n- 2) When $\\mu = \\pm 7.5$ then $d_{JS}$ is around $0.7$, When $\\vert \\mu \\vert > 7.5$, we can see that the `JS-divergence` in this case is stable and bounded by $0.7 ??$ \n\n**Explain.**\n\n2) Well, this is OK. In fact, the upper-bound of `JS-divergence` is $\\log(2) \\approx 0.693$ (the natural logarithm) for all distributions.\n\nIndeed, for any probability distributions, we have\n$$ \\begin{array}{cll} d_{JS}(p, q) &=& \\displaystyle \\frac{1}{2} \\int \\left[ p \\log \\left( \\frac{2p}{p+q} \\right) + q \\log\\left( \\frac{2q}{p+q} \\right) \\right] d\\nu \\\\ & \\leq & \\displaystyle \\frac{1}{2} \\max \\left\\lbrace \\log \\frac{2p}{p+q}, \\log \\frac{2q}{p+q} \\right\\rbrace \\left( \\int (p+q) d\\nu \\right) \\leq \\log 2. \\end{array} $$\n\nThis distance is always bounded even if the `KL-divergence` is infinity. As we known, for $p = \\mathcal{N}(0, 1)$ and $q = \\mathcal{E}(1)$, then\n$$d_{KL}(q \\Vert p) = \\sqrt{2\\pi} \\int_0^{\\infty} \\left( \\frac{x^2 - 2x}{2} \\right) e^{-x} dx \\approx 1.608$$\nis well-defined, while\n$$d_{KL}(p \\Vert q) = \\infty $$\nbut \n$$ d_{JS}(p, q) \\approx 0.47 $$","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"p = sst.norm.pdf(x, 0, 1)\nq = sst.expon.pdf(x, 1)\nJS_divergence(p, q)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Remark.** \n- 1) In fact, the `Kullback Leibler divergence` is used more often than the `Jensen Shannon` even if the `Kullback-Leibler` has a few nice properties. \n- 2) The square root of `JS-divergence` is defined a metric.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## 3. Wasserstein distance\n\n**Introduction.**\n\nThe `Earth-Move, (EM)` distance or `Wasserstein_1` distance, between the distributions $p, q$ is defined as\n$$ d_W(p, q) = \\inf_{\\tau \\in J(p, q)} \\mathbb{E}_{(X, Y) \\sim \\tau} \\Vert X - Y \\Vert. $$\nwhere $J(p, q)$ is the set of (probability) distributions $\\tau(x, y)$ whose marginals are respectively $p, q.$\n\nIn case the 1D-distributions, then the previous formula becomes\n$$ d_W(p, q) = \\inf_{\\tau \\in J(p, q)} \\mathbb{E}_{(X, Y) \\sim \\tau} \\vert X - Y \\vert. $$\n\nIntuitively, the terms $\\tau(x, y)$ indicates how much *`mass`* must be transported from $x$ to $y$ in order to transform the distribution $p$ to $q$. So, this distance is also knowns as the `EM`distance since it can be seen as the minimum amount of `“work”` required to transform $p$ into $q$; where the `“work”` is measured as the amount of distribution weight that must be moved, multiplied by the distance it has to be moved.\n\nRecall that, we defined $W_K$ distance as\n$$ d_{W_K}(\\mu, \\nu) = \\left( \\inf_{\\tau \\in J(p, q)} \\mathbb{E}_{(X, Y) \\sim \\tau} \\Vert X - Y \\Vert^K \\right)^{1/K}. $$","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.stats import wasserstein_distance as ssw\nimport numpy as np\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Example 1.** Use the function `wasserstein_distance` in `scipy.stats` to `approximate` the EM distance between 2 *independent* `uniform distributions`, such as $U([1, 2])$ and $U([0, 1])$.\n\n*Note that this function only using for calculating the distance of 1D distributions.*\nReference: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wasserstein_distance.html#scipy.stats.wasserstein_distance","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Fix the range of values\nx = np.linspace(-5, 5, 10)\n\n# Generate the empirical distribution values of the 2 uniform distributions\np = np.random.uniform(1, 2, 1000)\nq = np.random.uniform(0, 1, 1000)\n\n# Approximate the W1_distance\nssw(p, q)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Explain EX1.** There is only one possible coupling of these two measures, that is $\\tau(x, y) \\sim U([0, 1] \\times [1, 2])$ and hence we get\n$$d_W(p, q) = \\int_0^1 \\int_1^2 |x - y| dy dx = 1.$$","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**Example 2.** Now, we will work now with `discrete uniform distributions` in 2D space (instead of 1D space as above). Let $Z \\sim U([0, 1]), p$ be the distribution of $(0, Z)$ and $q$ be the distribution of $(\\eta, Z).$\n\nThis is easy to check that the EM distance  between them is equal to $\\vert \\eta \\vert$.\n\nFirst, look at the illustration bellow, $p, q$ are uniform on the straight vertical line passing the origin and the point $\\theta$ respectively.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"np.random.seed(42)\neta = 0.5\nn_points = 5\n\nZ = np.random.uniform(0,1, n_points)\np = np.array([[0, z] for z in Z ])\nq = np.array([[eta, z] for z in Z ])\n\nplt.figure(figsize=(6, 3))\nplt.scatter(p[:, 0], p[:, 1], label='$p$')\nplt.scatter(q[:, 0], q[:, 1], label='$q$')\nplt.xlabel('$\\eta = $' + str(eta))\nplt.legend(loc = 'upper right');\nplt.axis([-0.05, 1, 0, 1])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we will approximate this distance by the empirical distribution.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"p = np.array([[0, z] for z in Z ])\nq = np.array([[eta, z] for z in Z ])\n\n# calculate the distance for 2D distributions\nfrom scipy.spatial.distance import cdist\nfrom scipy.optimize import linear_sum_assignment\nd = cdist(p, q)\nassignment = linear_sum_assignment(d)\nprint(d[assignment].sum() / n_points)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, I take $\\theta \\in [-1, 1]$","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"Theta = [-1, -0.2, 0, 0.25, 1]\nn_points = 500\nD_Theta = []\nfor theta in Theta:\n    Z = np.random.uniform(0,1, n_points)\n    p = np.array([[0, z] for z in Z ])\n    q = np.array([[theta, z] for z in Z ])\n    d = cdist(p, q)\n    d_theta = d[linear_sum_assignment(d)].sum() / n_points\n    D_Theta.append(d_theta)\n\nplt.plot(Theta, D_Theta);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"","execution_count":null}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}