01 - Mallows Kendall, top-k and complete rankings.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mallows Kendall\n",
    "### A python package to for Mallows Model with top-$k$ and complete rankings using Kendall's-$\\tau$ distance\n",
    "&emsp;&emsp; **By Ahmed Boujaada, Fabien Collas and Ekhine Irurozki**\n",
    "\n",
    "We present methods for inference in the Mallows Model (MM), the best-known distribution for permutations. This is short tutorial for top-$k$ rankings and complete rankings under the Kendall's-$\\tau$ distance. Theoretical details are given in\n",
    "\n",
    "> Fabien Collas and Ekhine Irurozki (2020). Concentric mixtures of Mallows Models for top-k rankings: sampling and identifiability. In: International Conference on Machine Learning (ICML 21). 2021."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import mallows_kendall as mk\n",
    "import permutil as pu\n",
    "import scipy as sp\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kendall's-$\\tau$ Distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\sigma$ and $\\pi$ be two permutations of $n$ items. The Kendall's-$\\tau$ distance $d(\\sigma, \\pi)$ counts the number of pairwise disagreements between $\\sigma$ and $\\pi$. When $\\pi = e$, the Kendall's-$\\tau$ distance between two rankings $\\sigma$ and $\\pi$ is equal to the number of inversions in the ranking $\\sigma$. The Kendall's-$\\tau$ distance between two permutations can be computed as follows: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "perm1 = np.array([3, 1, 2, 0, 4])\n",
    "perm2 = np.array([3, 1, 4, 2, 0])\n",
    "mk.distance(perm1, perm2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function builds upon Merge sort algorithm and runs in $O(n\\log n)$. If only one permutation is given as input, it will be assumed that the second permutation is the identity permutation $e = (1, 2, \\dots, n)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mk.distance(perm1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maximum value of the kendall's-$\\tau$ distance between two permutations of length $n$ is $\\frac{n(n-1)}{2}$. It is possible to get this value using the following function: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 5\n",
    "mk.max_dist(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each $\\sigma \\in S_n$ there exists a unique vector, called inversion vector of $\\sigma$ with regards to the Kendall's-$\\tau$ distance, $V(\\sigma) = (V_1(\\sigma), \\dotsc, V_{n-1}(\\sigma))$ and such that $d(\\sigma) = \\sum_{j=1}^{n-1}V_j(\\sigma)$.   \n",
    "The $j^{th}$ element $V_j(\\sigma)$ is given by: $V_j(\\sigma) = \\sum\\limits_{i=j+1}^{n}\\mathbb{1}_{\\sigma(i) < \\sigma(j)}$, $\\forall j \\in \\{1, \\dotsc, n-1\\}$. It follows that $0 \\le V_j(\\sigma) \\le n-j$. \n",
    "\n",
    "There is a bijection between each $\\sigma \\in S_n$ and each possible $V(\\sigma)$ inversion vector. In the package the conversion from one form to another is done using the two following functions:  \n",
    "\n",
    "- From $\\sigma$ to $V(\\sigma)$ (here $V$ is of the same length $n$ as the permutation, thus the final $n^{th}$ element of $V$ will always equal $0$) :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3, 1, 1, 0, 0])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V_perm1 = mk.ranking_to_v(perm1)\n",
    "V_perm1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- From $V(\\sigma)$ to $\\sigma$ :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3, 1, 2, 0, 4])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mk.v_to_ranking(V_perm1, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],\n",
       "       [ 1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],\n",
       "       [ 1,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0],\n",
       "       [ 1,  2,  2,  1,  0,  0,  0,  0,  0,  0,  0],\n",
       "       [ 1,  3,  5,  6,  5,  3,  1,  0,  0,  0,  0],\n",
       "       [ 1,  4,  9, 15, 20, 22, 20, 15,  9,  4,  1]], dtype=uint64)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mk.num_perms_at_dist(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This package also allows to sample a permutation with a given number of inversions, i.e., at a given distance, where all the possible permutations have the same probability of being generated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0, 4, 3, 2])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mk.sample_at_dist(n, dist= 4, sigma0=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*note*: We point out that this notation is taken from Meilǎ, M., Phadnis, K., Patterson, A., & Bilmes, J. (2007). Consensus ranking under the exponential model. In *Proceedings of the 23rd Conference on Uncertainty in Artificial Intelligence, UAI 2007* (pp. 285–294). Note that the definition of the inversion vector difers from the one in the original paper of MM Fligner, M. A., & Verducci, J. S. (1986). Distance based ranking models. Journal of the Royal Statistical Society, 48(3), 359–369.\n",
    "This different versions imply different expressions for the MM which will affect the capability of the MM to be used for partial permutations.\n",
    "\n",
    "# Mallows Model (MM) for complete permutations\n",
    "The probability mass function of a Mallows Model with central permutation $\\sigma_0$ and dispersion parameter $\\theta$ is $$p(\\sigma) = \\frac{exp(-\\theta d(\\sigma, \\sigma_0))}{\\psi},$$\n",
    "\n",
    "and can be computed using the following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.008161793281714687"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 5    \n",
    "sigma = np.array([3,1,2,0,4])\n",
    "sigma_0 = np.array(range(5))\n",
    "theta = 0.1\n",
    "\n",
    "mk.prob(sigma, sigma_0, theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sampling** This package includes a sampler based on the factorization of the Kendall's-$\\tau$  distance. In the later sections we will also present how this can be adapted to top-$k$  rankings. This differs to the classical sampling, usually done using the Repeated Insertion Model (RIM) and which can not be extended to top-$k$  rankings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2, 0, 3, 1, 4],\n",
       "       [2, 0, 3, 1, 4],\n",
       "       [0, 1, 2, 3, 4],\n",
       "       [1, 3, 0, 2, 4]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mk.sample(m=4,n=5,theta=1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in the package the sampling functions generates the samples considering $\\sigma_0 = e$, identity permutation by default. But any other central permutation can be given as a parameter.  \n",
    "In practice, we can draw a sample from a MM as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3, 4, 2, 1, 0],\n",
       "       [3, 4, 2, 1, 0],\n",
       "       [4, 3, 2, 1, 0],\n",
       "       [4, 3, 1, 2, 0]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mk.sample(m=4, n=5, theta=1.5, s0=np.array([4,3,2,1,0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this package, we can specify also the parameter `phi` instead of `theta` when we consider the following equivalent expression for the Mallows Model, $$p(\\sigma) \\propto \\phi^{d(\\sigma,\\sigma_0)}.$$ The relation between both expressions is given by $\\phi = \\exp(-\\theta)$. This functionality holds for most functions. The sampling, for example, is done then as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4, 0, 1, 2, 3],\n",
       "       [1, 2, 0, 3, 4],\n",
       "       [1, 0, 2, 3, 4],\n",
       "       [0, 3, 2, 4, 1]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mk.sample(m=4,n=5,phi=.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Expected distance** The expected value of Kendall's-$\\tau$ distance under the MM is given by: \n",
    "$$\\mathbb{E}[D] = \\frac{n \\cdot \\exp(-\\theta)}{1 - \\exp(-\\theta)} - \\sum_{j=1}^n\\frac{j \\cdot \\exp(-j \\theta)}{1 - \\exp(-j \\theta)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.4578024409695693"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta_mm = 0.7\n",
    "expected_dist = mk.expected_dist_mm(n, theta_mm)\n",
    "expected_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Variance** The variance of Kendall's-$\\tau$ distance under the MM can be expressed as follows: \n",
    "$$\\mathbb{V}[D] = \\dfrac{ n \\cdot \\exp(-\\theta) }{ (1 - \\exp(-\\theta))^2 } - \\sum_{j=1}^n \\dfrac{ j^2\\exp(-j \\theta) }{ (1 - \\exp(-j \\theta))^2 }$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.763301117320907"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variance_dist = mk.variance_dist_mm(n, theta_mm)\n",
    "variance_dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_mm = mk.sample(m=4,n=5,phi=.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Learning** Fitting a Mallows model Model is done in two stages. First, estimate the central permutation $\\sigma_0$ and, second, compute the dispersion parameter. \n",
    "\n",
    "Function `median` uses Borda algorithm to approximate the central permutation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 3, 2, 4, 1])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mk.median(sample_mm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `fit_mm` returns an approximation to the MLE of the parameters $\\sigma_0$ and $\\theta$. The consensus $\\sigma_0$ is approximated with the well known Borda count algorithm. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0, 3, 2, 4, 1]), 0.5943115810463766)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma_borda, phi_mm_mle = mk.fit_mm(sample_mm)\n",
    "sigma_borda, phi_mm_mle "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mallows Model for Top-$k$ rankings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sampling** In this part the sampling of the top-$k$ rankings is discussed.  \n",
    "Perhaps the most natural idea to generate top-$k$ rankings would be to generate the full rankings, using the Repeated Insertion Model for example, and to cut these obtained permutations after position $k$. This is possible, yet it seems computationally non-optimal with a complexity of $O(n^2)$.    \n",
    "\n",
    "Here we adapt the methods of previous sections to top-$k$ rankings. This method is based on sampling partially the inversion vectors. In practice in the package, a top-$k$ ranking can be sampled as follows: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  3., nan, nan,  4., nan,  1., nan,  2., nan],\n",
       "       [nan,  2.,  1.,  0., nan, nan,  4.,  3., nan, nan],\n",
       "       [ 4.,  3.,  1.,  0., nan, nan, nan, nan, nan,  2.],\n",
       "       [ 3.,  2., nan, nan,  4.,  0., nan, nan, nan,  1.],\n",
       "       [ 0.,  3., nan, nan, nan, nan, nan,  4.,  1.,  2.],\n",
       "       [ 1.,  4.,  2., nan, nan, nan, nan,  3., nan,  0.],\n",
       "       [ 4.,  0., nan, nan, nan, nan,  3.,  1., nan,  2.],\n",
       "       [nan,  3.,  1., nan, nan, nan,  0.,  4.,  2., nan],\n",
       "       [ 3., nan, nan, nan,  2., nan,  0.,  1.,  4., nan],\n",
       "       [nan,  0.,  3., nan, nan, nan,  1.,  4.,  2., nan],\n",
       "       [ 2., nan, nan,  0.,  3.,  4.,  1., nan, nan, nan],\n",
       "       [nan, nan, nan,  4.,  0., nan,  1.,  2.,  3., nan],\n",
       "       [ 3., nan,  1.,  4., nan, nan, nan, nan,  0.,  2.],\n",
       "       [ 3., nan, nan,  2., nan, nan,  0.,  4., nan,  1.],\n",
       "       [ 0., nan, nan, nan, nan,  1., nan,  3.,  2.,  4.]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = 15\n",
    "n = 10 \n",
    "k = 5\n",
    "phi = 0.9 \n",
    "sigma_0 = np.array(range(10))\n",
    "sigma = np.random.permutation(n)\n",
    "sample_top_k = mk.sample(m=m, n=n, k=k, phi=phi, s0=sigma_0)\n",
    "sample_top_k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Distance** \n",
    "\n",
    "To compare top-$k$ rankings, we include an extension of the classical Kendall's-$\\tau$ denoted as $p$-parametrized Kendall's-$\\tau$ distance (Fagin et al., 2003). In practice we can use it, choosing the $p$ parameter in $[0,1]$, as follows: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0.,  3., nan, nan,  4., nan,  1., nan,  2., nan]),\n",
       " array([nan,  2.,  1.,  0., nan, nan,  4.,  3., nan, nan]))"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_top_k[0],sample_top_k[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mk.p_distance(sample_top_k[0], sample_top_k[1], k=k, p=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mk.p_distance(sample_top_k[0], sample_top_k[1], k=k, p=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**Expected distance** The expected value of Kendall's-$\\tau$ distance for top-$k$ rankings under the MM is given by: \n",
    "$$\\mathbb{E}[D] = \\frac{k \\cdot \\exp(-\\theta)}{1 - \\exp(-\\theta)} - \\sum_{j=n-k+1}^n\\frac{j \\cdot \\exp(-j \\theta)}{1 - \\exp(-j \\theta)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.732931959531709"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expected_dist = mk.expected_dist_top_k(n,k,theta_mm)\n",
    "expected_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Variance** The variance of Kendall's-$\\tau$ distance for top-$k$ rankings under the MM can be expressed as follows: \n",
    "$$\\mathbb{V}[D] = \\dfrac{ k \\cdot \\exp(-\\theta) }{ (1 - \\exp(-\\theta))^2 } - \\sum_{j=n-k+1}^n \\dfrac{ j^2\\exp(-j \\theta) }{ (1 - \\exp(-j \\theta))^2 }$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.391580684957216"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variance_dist = mk.variance_dist_top_k(n,k,theta_mm)\n",
    "variance_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized Mallows Model (GMM) for complete permutations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Generalized Mallows Model (GMM) is an extension of the Mallows Model, for which there are $n-1$ dispersion parameters $\\theta_j$ ($1 \\le j < n$), each affecting a particular position of the permutation.   \n",
    "Formally, the GMM under Kendall's-$\\tau$ distance is expressed as follows: \n",
    "\n",
    "$$p(\\sigma)=\\dfrac{\\exp(\\sum_{j=1}^{n-1}-\\theta_j V_j(\\sigma\\sigma_0^{-1}))}{\\psi(\\theta)}$$\n",
    "where $\\psi(\\theta) = \\prod_{j=1}^{n-1} \\psi_j(\\theta_j) = \\prod_{j=1}^{n-1} \\frac{1-\\exp(-\\theta_j(n-j+1))}{1 - \\exp(-\\theta_j)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00439275005307042\n"
     ]
    }
   ],
   "source": [
    "def prob_GMM(sigma, sigma_0, theta):\n",
    "    n = len(sigma)\n",
    "    sigma_0_inv = pu.inverse(sigma_0)\n",
    "    V = mk.ranking_to_v(pu.compose(sigma, sigma_0_inv))\n",
    "    psi = np.prod(np.array([(1 - np.exp(( - n + j ) * theta[j]))/(1 - np.exp(-theta[j])) for j in range(n-1)]))\n",
    "    return np.exp( np.sum ( [ -theta[j] * V[j] for j in range(n-1) ] ) ) / psi \n",
    "    \n",
    "sigma = np.array([3,1,2,0,4])\n",
    "sigma_0 = np.array(range(5))\n",
    "theta = [0.5,0.2,0.6,0.3]\n",
    "\n",
    "print(prob_GMM(sigma, sigma_0, theta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sampling** It is also possible to sample the Generalized Mallows Models (GMM), using the same process as for the MM only with a change in the probability of $V_j$ as follows:\n",
    "\n",
    "$$p(V_j(\\sigma\\sigma_0^{-1}) = r) = \\frac{\\exp(-\\theta r)}{\\psi_j(\\theta)}, \\; \\forall r \\in {0, \\dotsc, n-j}$$\n",
    "\n",
    "In the package it is quite similar as for the classical Mallows Models, only that the dispersion parameter will be given as a list. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 4000\n",
    "n = 5\n",
    "theta_gmm = [0.5,0.2,0.6,0.3]\n",
    "identity = np.array(range(n))\n",
    "sample_gmm = np.array(mk.sample(m, n, theta = theta_gmm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expected value of each term $V_j(\\sigma)$, where $\\sigma$ is a random Mallows permutation, is expressed as follows:\n",
    "$$ \\mathbb{E}[V_j] = \\frac{\\exp(-\\theta_j)}{1 - \\exp(-\\theta_j)} - \\frac{(n-j+1)\\exp(-\\theta_j(n-j+1))}{1 - \\exp(-\\theta_j(n-j+1))}, \\;\\forall j \\in \\{1, \\dotsc, n-1\\} $$\n",
    "And the values are computed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.09436663 1.25279068 0.62226834 0.42555748]\n"
     ]
    }
   ],
   "source": [
    "expected_v_gmm = mk.expected_v(n, theta_gmm)\n",
    "print(expected_v_gmm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expected value of Kendall's-$\\tau$ distance under the GMM is then given by: \n",
    "$$\\mathbb{E}[D] = \\sum_{j=1}^{n-1} \\mathbb{E}[V_j]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.394983134634349\n"
     ]
    }
   ],
   "source": [
    "expected_dist_gmm = np.sum(expected_v_gmm)\n",
    "print(expected_dist_gmm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variance of each term $V_j(\\sigma)$ under the GMM can be expressed as follows: \n",
    "$$\\mathbb{V}[V_j] = \\dfrac{ \\exp(-\\theta_j) }{ (1 - \\exp(-\\theta_j))^2 } - \\sum_{j=1}^n \\frac{ (n-j+1)^2 \\exp(-(n-j+1) \\theta_j) }{ (1 - \\exp(-(n-j+1) \\theta_j))^2 }, \\;\\forall j \\in \\{1, \\dotsc, n-1\\} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.48213789 1.20855956 0.56066479 0.24445831]\n"
     ]
    }
   ],
   "source": [
    "variance_v_gmm = mk.variance_v(n, theta_gmm)\n",
    "print(variance_v_gmm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expected value of Kendall's-$\\tau$ distance under the GMM is then given by: \n",
    "$$\\mathbb{V}[D] = \\sum_{j=1}^{n-1} \\mathbb{V}[V_j]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.495820543795291\n"
     ]
    }
   ],
   "source": [
    "variance_dist_gmm = np.sum(variance_v_gmm)\n",
    "print(variance_dist_gmm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Learning** An approximation to the MLE for $\\theta$ and $\\sigma_0$ given a sample of i.i.d Mallows permutations is as follows: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 2 1 3 4]\n",
      "[0.5109401970612621, 0.19755236580605012, 0.616607960418165, 0.2757338026671021]\n"
     ]
    }
   ],
   "source": [
    "sigma_borda, theta_gmm_mle = mk.fit_gmm(sample_gmm)\n",
    "print(sigma_borda)\n",
    "print(theta_gmm_mle)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}