How to sample N positive integers whose sum equals M

sampling_n_positive_integers_summing_to_m.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample $N$ positive integers that sum to $M$ such that the distribution of the $N$ integers is uniform over all possible integer sets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A solution is presented by Smith and Tromble in [\"Sampling Unfiromly from the Unit Simplex\"](https://www.cs.cmu.edu/%7Enasmith/papers/smith+tromble.tr04.pdf)\n",
    "\n",
    "The algorithm works as follows:\n",
    "\n",
    "1. Sample $N-1$ random positive integers  from $1$ to $M-1$ (inclusive) and sort them to produce the sequence $x_1 \\ldots x_{N-1}$.\n",
    "2. Compute the difference between each sequence element: $y_i = x_i - x_{i-1}$ adding the additional sequence elements $x_0 = 0$ and $x_M = M$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from collections import Counter\n",
    "from tqdm import trange\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 5, 3])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def sample_ints_without_replacement(low, high, size=1):\n",
    "    rng = np.random.default_rng()\n",
    "    return rng.choice(high - low, size=size, replace=False) + low\n",
    "\n",
    "def uniform_n_integers_that_sum_to_m(n, m):\n",
    "    points = np.concatenate([[0], np.sort(sample_ints_without_replacement(1, m, size=n-1)), [m]])\n",
    "    gaps = points[1:] - points[:-1]\n",
    "    return gaps\n",
    "    \n",
    "uniform_n_integers_that_sum_to_m(4, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Verifying uniformity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000000/1000000 [01:38<00:00, 10108.54it/s]\n"
     ]
    }
   ],
   "source": [
    "n = 4\n",
    "m = 10\n",
    "n_runs = 1e6\n",
    "\n",
    "combination_counts = Counter()\n",
    "for _ in trange(int(n_runs)):\n",
    "    combination_counts[tuple(uniform_n_integers_that_sum_to_m(n, m))] += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that the sum of all the integers is $M$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert (np.array(list(map(sum, combination_counts.keys()))) == m).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check distribution of the sequences is uniform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:xlabel='seq'>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3oAAAFmCAYAAAA72OixAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnPklEQVR4nO3df7RcZX3v8fcXEk2qiYTk5GKJEMGApVhY3iClWH4o1kLA3srytnihVJEsC6UiDQopQmsBg4BitepFUeilP0SpwgJ7i1aQmnKpoYK2iT8qiIQWEoGQWAMG/d4/Zh8d55yTM2dmzpk9z7xfa+11zjz7M3s/z/6VfM/M7InMRJIkSZJUjl363QFJkiRJUm9Z6EmSJElSYSz0JEmSJKkwFnqSJEmSVBgLPUmSJEkqjIWeJEmSJBVmVr870KlFixbl0qVL+90NSZIkSeqLe+6553uZOTLevIEt9JYuXcq6dev63Q1JkiRJ6ouIeHCieb51U5IkSZIKY6EnSZIkSYWx0JMkSZKkwljoSZIkSVJhLPQkSZIkqTAWepIkSZJUGAs9SZIkSSqMhZ4kSZIkFcZCT5IkSZIKY6EnSZIkSYWZNVkgIvYALgYOysxDqrb3Aj8Avg8cBJydmY9U884F5gMLgNsy8+aq/WDgTOABYDGwKjOfiYg5wBXAw8AyYE1mfrOXg5QkSZI6tfS8W8e0fWfNij70RGrfpIUe8HLgJuDgprb/yswLACLi7cAfAWdFxKHA0Zl5XETMBtZHxJ3Ak8D1wDGZ+UhEXAmcClwDnA18NzPfHREvqdp+tSejU8/0+gJX9wtm3fsnSYPAa6lUL56Tw2XSQi8zPxURR7W0XdD0cBcar+wBHA/cVWV2RMQG4Ajg34C5o6/6AWuBk2kUdSuA1dVzvhYRB0XE/Mzc2umgJEkNw/aPer/G63Yue7waq5RjwGtGvbhdequdV/QmFBG7Ab8GnFg1LQY2NEW2Vm2bgW3jtI8+Z7x5Ywq9iFgJrATYa6+9uul61zwQpfop4bxsdwwljHUQuJ1nxrBt52Ebr1SKQTt3Oy70IuJ5wAeBN2bm41XzJmBeU2x+1TZR+86eM0ZmXg1cDbB8+fLstO8zyf+0dafX28/tPL5he2tur5UwXv+qXTavkTPD7Tf96r6N694/DZeOCr2IWARcBZybmQ9HxImZeSNwC3BRlZkFHACMfkZve0TsUb1983Bg9Ey4FTgM+MfqM3r3TdfbNofp5PM/bQL3x0TcLp1z22kYlfBvat3H4LVlMLl/662du24eCZwCPD8iLgCuBG6rnvuXEQGNt17emJl3R8TtEXEpjbtunpOZW6rlnAxcEhEPArsC11WreB9wRbXsFwGn9XB8xfGVF0kaXF5zVXelHKOljKNf6r796v6Or7osr52bsXwR+GJL80t3kr98gvZ7GaeIy8ztNL52QRo4df8LqeuVxlf3Y8q/knfH7aJSeCyXbbr3b1c3Y6mLEt7aoLJ5XI1vmLbLMI0V6vPXTGlnPK5UCo9ljaeIQk8qhRdqSZqY10hJap+FniRNwP9USqo7r1OSJrJLvzsgSZIkSeotCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMLMmC0TEHsDFwEGZeUjVtjuwBrgfWAaszsxHq3nnAvOBBcBtmXlz1X4wcCbwALAYWJWZz0TEHOAK4OFqWWsy85u9HKQkSZIkDZNJCz3g5cBNwMFNbZcCn8/MGyLiBBqF2ikRcShwdGYeFxGzgfURcSfwJHA9cExmPhIRVwKnAtcAZwPfzcx3R8RLqrZf7c3wJEmSJGn4TPrWzcz8FLCtpXkFcFf1+9rqMcDxo+2ZuQPYABwB7APMzcxHxnnOiqbnfA04KCLmdzIYSZIkSVLnn9FbzE+Lv63AgoiY1dI+Om/xTtpbl9U672dExMqIWBcR6zZv3txh1yVJkiSpbJ0WepuAedXv84EnMvOZlvbReZt20t66rNZ5PyMzr87M5Zm5fGRkpMOuS5IkSVLZOi30bgUOq34/vHoMcMtoe/UK3wHAnTRu2rK9urFL63NubXrOS4D7MnNrh/2SJEmSpKHXzl03jwROAZ4fERcAVwKrgcsiYj9gX2AVQGbeHRG3R8SlNO66eU5mbqmWczJwSUQ8COwKXFet4n3AFdWyXwSc1sPxSZIkSdLQmbTQy8wvAl9sad4OnD5B/vIJ2u9lnCIuM7fT+NoFSZIkSVIP+IXpkiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVZlY3T46Ic4GlwPeAZcBpwFxgDXB/1bY6Mx9tys8HFgC3ZebNVfvBwJnAA8BiYFVmPtNN3yRJkiRpWHVc6EXEHsD5wKLM/HFE3AS8FvhV4POZeUNEnABcAZwSEYcCR2fmcRExG1gfEXcCTwLXA8dk5iMRcSVwKnBNd0OTJEmSpOHUzVs3fwD8kMYrdADPBf4NWAHcVbWtrR4DHD/anpk7gA3AEcA+wNzMfGSc5/yMiFgZEesiYt3mzZu76LokSZIklavjQi8ztwLnAp+IiGuBjcC/03jr5bYqthVYEBGzWtpH5y3eSft467w6M5dn5vKRkZFOuy5JkiRJReu40Ks+V3cusCIzf5fG5/QuBDYB86rYfOCJ6vN2ze2j8zbtpF2SJEmS1IFu3rq5J/B4001T/hOYA9wKHFa1HV49BrhltL16he8A4E4aN23ZXn3mr/U5kiRJkqQp6uaum/8XOK66ecoW4EDgbOBp4LKI2A/YF1gFkJl3R8TtEXEpjbtunpOZWwAi4mTgkoh4ENgVuK6LfkmSJEnSUOu40MvMH9H4SoTxnD7Bcy6foP1eGl/NIEmSJEnqkl+YLkmSJEmFsdCTJEmSpMJY6EmSJElSYSz0JEmSJKkwFnqSJEmSVBgLPUmSJEkqjIWeJEmSJBXGQk+SJEmSCmOhJ0mSJEmFsdCTJEmSpMJY6EmSJElSYSz0JEmSJKkwFnqSJEmSVBgLPUmSJEkqjIWeJEmSJBXGQk+SJEmSCmOhJ0mSJEmFsdCTJEmSpMJY6EmSJElSYSz0JEmSJKkwFnqSJEmSVBgLPUmSJEkqjIWeJEmSJBXGQk+SJEmSCmOhJ0mSJEmFsdCTJEmSpMJY6EmSJElSYSz0JEmSJKkwFnqSJEmSVBgLPUmSJEkqjIWeJEmSJBXGQk+SJEmSCmOhJ0mSJEmFsdCTJEmSpMLM6ubJEbE/cBKwHTgS+GPg34E1wP3AMmB1Zj5a5c8F5gMLgNsy8+aq/WDgTOABYDGwKjOf6aZvkiRJkjSsOi70ImJX4D3ACZn544j4C+AZ4FLg85l5Q0ScAFwBnBIRhwJHZ+ZxETEbWB8RdwJPAtcDx2TmIxFxJXAqcE13Q5MkSZKk4dTNWzcPAQI4KyLOB04AvgesAO6qMmurxwDHj7Zn5g5gA3AEsA8wNzMfGec5kiRJkqQp6uatm3sDhwEnZeaTEXE98EMab73cVmW2AgsiYlbVvqHp+Vurts1N+eb2MSJiJbASYK+99uqi65IkSZJUrm5e0dsKfD0zn6wefwk4CtgEzKva5gNPVJ+3a24fnbdpJ+1jZObVmbk8M5ePjIx00XVJkiRJKlc3hd7dwMLqs3rQeIXvm8CtNF7pAzi8egxwy2h79QrfAcCdNG7asj0i9hjnOZIkSZKkKer4rZuZ+XhEvB24KiI2AyPAO4G5wGURsR+wL7Cqyt8dEbdHxKU07rp5TmZuAYiIk4FLIuJBYFfgui7GJEmSJElDrauvV8jMTwOfbmneDpw+Qf7yCdrvBU7rpi+SJEmSpAa/MF2SJEmSCmOhJ0mSJEmFsdCTJEmSpMJY6EmSJElSYSz0JEmSJKkwFnqSJEmSVBgLPUmSJEkqjIWeJEmSJBXGQk+SJEmSCmOhJ0mSJEmFsdCTJEmSpMJY6EmSJElSYSz0JEmSJKkwFnqSJEmSVBgLPUmSJEkqjIWeJEmSJBXGQk+SJEmSCmOhJ0mSJEmFsdCTJEmSpMJY6EmSJElSYSz0JEmSJKkwFnqSJEmSVBgLPUmSJEkqjIWeJEmSJBXGQk+SJEmSCmOhJ0mSJEmFsdCTJEmSpMJY6EmSJElSYSz0JEmSJKkwFnqSJEmSVBgLPUmSJEkqjIWeJEmSJBXGQk+SJEmSCmOhJ0mSJEmFsdCTJEmSpMLM6nYBETEXuBu4LTNXRcTuwBrgfmAZsDozH62y5wLzgQVV/uaq/WDgTOABYDGwKjOf6bZvkiRJkjSMui70gIuBrzQ9vhT4fGbeEBEnAFcAp0TEocDRmXlcRMwG1kfEncCTwPXAMZn5SERcCZwKXNODvkmSJEnS0OnqrZsRcQqwlsYrcaNWAHdVv6+tHgMcP9qemTuADcARwD7A3Mx8ZJznSJIkSZKmqONCLyIOAH4hM/+2ZdZiYFv1+1ZgQUTMamkfnbd4J+3jrXNlRKyLiHWbN2/utOuSJEmSVLRuXtH7TeCpiDgPeDnwsog4G9gEzKsy84Enqs/bNbePztu0k/YxMvPqzFyemctHRka66LokSZIklavjz+hl5iWjv0fEHOC5mXlVRLwYOAx4CDgcuLWK3QJcVOVnAQcAo5/R2x4Re1Rv32x+jiRJkiRpinpx180TaXzW7lkRcRKwGrgsIvYD9gVWAWTm3RFxe0RcSuOum+dk5pZqGScDl0TEg8CuwHXd9kuSJEmShlXXhV5m3gjc2NJ8+gTZyydovxc4rdu+SJIkSZL8wnRJkiRJKo6FniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhbHQkyRJkqTCWOhJkiRJUmEs9CRJkiSpMBZ6kiRJklQYCz1JkiRJKoyFniRJkiQVxkJPkiRJkgozq9MnRsS+wMXAvwBLgMcy850RsTuwBrgfWAaszsxHq+ecC8wHFgC3ZebNVfvBwJnAA8BiYFVmPtNp3yRJkiRpmHVc6AG7A3+TmTcBRMT6iLgVOB34fGbeEBEnAFcAp0TEocDRmXlcRMwG1kfEncCTwPXAMZn5SERcCZwKXNNF3yRJkiRpaHX81s3M/PJokde0rP8CVgB3VW1rq8cAx4+2Z+YOYANwBLAPMDczHxnnOZIkSZKkKerJZ/Qi4jeBv8/Mr9N46+W2atZWYEFEzGppH523eCft461nZUSsi4h1mzdv7kXXJUmSJKk4XRd6EXE0cDTw1qppEzCv+n0+8ET1ebvm9tF5m3bSPkZmXp2ZyzNz+cjISLddlyRJkqQidVXoRcQK4NXAW4A9IuIw4FbgsCpyePUY4JbR9uoVvgOAO2nctGV7ROwxznMkSZIkSVPUzV03/zvwCWAdcDvwHODPgdXAZRGxH7AvsAogM++OiNsj4lIad908JzO3VMs6GbgkIh4EdgWu63hEkiRJkjTkOi70MvMe4LkTzD59gudcPkH7vcBpnfZFkiRJkvRTfmG6JEmSJBXGQk+SJEmSCmOhJ0mSJEmFsdCTJEmSpMJY6EmSJElSYSz0JEmSJKkwFnqSJEmSVBgLPUmSJEkqjIWeJEmSJBXGQk+SJEmSCmOhJ0mSJEmFsdCTJEmSpMJY6EmSJElSYSz0JEmSJKkwFnqSJEmSVJhZ/e6AJEnTbceOHWzcuJGnnnoKgI+85vljMhs2bBjTNtO5OXPmsGTJkjEZSZKmykJPklS8jRs3Mm/ePJYuXUpEsGPjljGZX1iy25i2mcxlJo899hgbN24ck5Ekaap866YkqXhPPfUUCxcuJCL63ZUJRQQLFy78yauOkiR1w0JPkjQU6lzkjRqEPkqSBoNv3ZQkDZ3XfGBtT5d38+8f3tPlSZLULV/RkyRJkqTCWOhJkiRJUmEs9CRJmgEf//jHecVL9+djH7yKd553Nm848Ti+v21rv7slSSqUhZ4kSTPgDW94Ay/cdxn7H/ASLlxzFctefAD/7x/v6He3JEmFstCTJGkG7b3PiwBYsHAR//X97/e5N5KkUlnoSZI0g/wKBUnSTPDrFSRJQ2e8r0P4pSW7jWn76sYtHedafe5zn+M/Hn6Iz3zien7jda/nnrv/iW99fT1vPuV1jIyMtNFrSZLaZ6EnSdIMeNWrXsXf/dN9P3n80U/cDMDIyG596pEkqWS+dVOSJEmSCmOhJ0mSJEmFsdCTJA2FzOx3FyY1CH2UJA0GCz1JUvHmzJnDY489VutCKjN57LHHmDNnTr+7IkkqgDdjkSQVb8mSJWzcuJHNmzcD8OgT28dkNmybO6ZtpnNz5sxhyZIlwPoxOUmSpsJCT5JUvNmzZ/PCF77wJ4+PPe/WMZnvrFkxpq1fOUmSulWbQi8ijgFeC2wCMjP/pM9dkiRJkqSBVItCLyJ+Dvgw8IuZ+XRE3BgRr8zMf+h33yRJkiRp0NTlZiyHAQ9m5tPV47WA72WRJEmSpA5EHe5AFhEnAb+Vmf+jevwm4KjMPLkltxJYWT3cH/hGy6IWAd9rY5XmzPU6V+e+mTNnbnBzde6bOXPmBjdX576Zm1pu78wcGTedmX2fgFcC/9D0+BzgPR0sZ505c/3I1blv5syZG9xcnftmzpy5wc3VuW/mus+NTnV56+ZdwN4R8ezq8eHA2FuTSZIkSZImVYubsWTmDyLi94A/i4jNwFfTG7FIkiRJUkdqUegBZObngM91uZirzZnrU67OfTNnztzg5urcN3PmzA1urs59M9d9DqjJzVgkSZIkSb1Tl8/oSZIkSZJ6xEJPkiRJkgpjoSdJkiRJhanNzVg6ERELgZOAXwTmAt8FbszM+8yZ6yRX576ZM2ducHN17ps5c+YGN1fnvpmbudxEBvZmLBFxFHARsBZ4GNgB7A78MrA+My8wZ24quTr3zZw5c4Obq3PfzJkzN7i5OvfN3Mzldiqn8O3qdZlovOX093cy/xjgQHPmppB7SY37Zs6cucHNeW0xZ86c1xZz05KbaP7oNLCv6NVFRMzOzB397sdMKWW8dR5Hnfum7pWyf+s+jrr3T2PVfZ/VvX/Dxv0hTW6XfnegUxHxhYh4fUTsOknuhIhYEhG7RsRF1fM+HBGLWnILm37/XxFxTUScHxFzJunKW6vnrGxdXvX834iGyyLihojYvyUXEfGbEfGqprYVEXHhJOsdzf63IR9vz8bRrzG007d+9q+P653u7TzRej2HBnAcde7fEJ5DHlNem9teXhe5od8f/doXw3bs1eiYn1IOGNxX9CJiLfAR4HU03rv64cx8fJzcJ4E3AiuB3YAvAXsDv5KZv9uU+0JmviIiTgaOB24HlgG7ZeabmnI3tKxif+AbwIsz85eacp8GHgT2pfG+2keBB4AVmfm6ptz7q8wuwFbg5GrW5zLzyDa2w+rMvHSIx9uzcfRxDHXfxv1ab0+XN4X1eg4N5jhq278hPIc8prw2t728KazX/cGYa0u/9sWwHXt9WW+3ORjsu24+nZnXAtdG48OKfx4RW4D3Z+b6ptxdmbktIsjMd4w2Rstf55vslZm/3ZR7Z8v85wH/DvxzU9utjH119F8y80+rZdyQmW+uft+zJbc1M4+r5h0C/G/gdODHzaGI+E+g9dXLAOYAzTt72Mbby3H0awx138b9Wm9Pl+c5NOF4SxlHnfs3bOfQKI+p6e3fsB1X7o+x6+3Xvhi2Y68v6+3gmjvGIBd6P5GZdwB3RMQLgbOAc5pmz4uI5wA/iIi9MvO7ETEP2GuCxf14kse/DpwLzMvMD0TEHpl5XUQ81JJrfovK5U2/t+6wJ5vG8eWIuBL4AGMvXNcBNwKbmtoCOLMlN2zj7eU4+jWGum/jfq2318vzHJr+c6if46hz/4btHBrlMTW9/Ru248r9MXa9/doXw3bs9Wu9U73mjpWT3K2lrhOwEXhdG7k9qo20DvghjZeYv0njbVjNuU8DFwIfBZZXbe8APjnBco8CPgxcMsH864B3tLSd2ro84K+AlS1thwBbWtrmA8eOs56Dhny8PRtHv8YwANu4X+vt9fI8h6b5HOrnOOrcvyE8hzymvDb3/Lhyf4x7benXOT5Ux14f1zulc2O8aWA/ozdV0fjg4t403lv7rcz8URvPeS6wIzOfnmD+YuD4zPxYm31YBDyTmVua2nYDnpOZD7dk92xtm4phG2/Lsno2jpkewwBt42lf70yMY5L1ew4N6Djq3r+ZWme/z6GWdXpMeW3uGfdH75fVrZKPvTpt56kamkJPkiRJkobFwH69giRJkiRpfBZ6kiRJklQYCz1JkiRJKkxxhV5EbI6IPzRnrle5OvfNnDlzg5urc9/MmTM3uLk6983czOWgwEIP2B+IiFhqzlyPcnXumzlz5gY3V+e+mTNnbnBzde6buZnLedfNbkXEvpn57YiYm5nbe7C8yHF2SkQ8KzN/2O3yZ8pMbZc66PVY+6WUcZRimM4hKOf4K2UcJShlX5QyjlK4PzRIintFLyJe1vJ4cUS8LyIuicb3YBARr42Im1pyvxwRt0XEeyNi94i4OyI2RcQJk6zyxOrnW1uW98rq59yI+GhEfDsi/j4iXtSS+9OIODAiRiLidmBHRHw9Ipa3rOc/IuJDEXFAWxuiRcHbpWfj7ddYezmG6RhHjfrX1XrbXV6v1ztM51A/x1v3cdSkb12tt93l9Xq9HlNem6djvSXvjyHaF12tt93l1X29O83W+A+8OxUTFzwrM/PsptxfA1+hUdS+CnhDZn43Im7PzKObcp8BrgVeBBwP/AnwHeCSzHx9U24z0LzR5gJPAXMyc15T7guZ+YqIeBuwDfgCsB9wUsvyPpSZvxcR7wS+CHwJWAqcn5m/25S7HXgj8AfAC4CPZeZn3S69G28fx1r3fdav/vV6ve0ur+77o7bnUJ/HW9txeA7N2HiH5pjqc/88rmqyP4ZwX5Qy3p6ud6cycyAn4FHgfuCBlumxltx5Tb/vBlwPPB/4Qkvu3Kbfrxrv+dXjC4ArgL1o/KdpDbA3cHlL7gvVzz9saV/d8vj3q59vmajfzcurfp8HnA18FjhjyLdLz8bbx7HWfZ/1q3+9Xm+7y6v7/qjtOeR5NP44+tg3z6FCjymPK/fHkO6LUsbb0/XubJrF4HoHcHdm3tfcGBF/1JJ7zugvmbklIs4A/hzYvSW3LCKenZlPA+9qal/YHMrMiyPiSOA84Hzg8cx8MCI+PEE/nx0/+z7ukZb5z42IlwOPR8TJwN3APtU0rszcBlwVEe8DXtMye9i2Sy/H26+x1n2f9aV/07DedpdX9/0xqo7nUD/HW+dxeA51t16PKa/N07HeEvbHUO2LaVhvv8bb6/VOrN2KsG4TEMBB47Tv1vL4PcA7W9oWAfe0tP0BsKqlbRXwZxOsfzHwUeCDE8y/CLiw+vkLVdu7abwM3ZybTeOvReuBHwD/SeMvAD/fkrsFeJ3bZcx26dl4+zjWuu+zvvRvGtbb7vLqvj9qew71eby1HYfnUH32RSnHlMeV+2MY90VB4+3penc2Dexn9OogIgJ4aWbe0+++1MkwbZdSxlrKOEoxbPujlPGWMo4SlLIvShlHKdwfGjQWepIkSZJUmOK+XkGSJEmShp2FniRJkiQVprhCbyffOWHOXEe5OvfNnDlzg5urc9/MmTM3uLk6983czOWAwb3r5kQTcB9wE7CLOXO9yNW5b+bMmRvcXJ37Zs6cucHN1blv5mYul5nlFXrVBjgU2NOcuV7l6tw3c+bMDW6uzn0zZ87c4Obq3DdzM5fzrpsqSkQsAhZn5vpuc71c1hRzR2fm7RHx4sz8eg1zLwAezcwfTjKOnuUi4nnAQcD9mbmxk1xEHJGZd+6sLypHH8/fdq4tA3/utntOauYU9G9WW8epxurXvi1Fu9tloExWCQ7aRBtfKl5yDjgM+CqwFnhZ1fZa4MFCc6cB24GvAC8Dvgt8G7hgqrleLmuKuQNapiurn1fVJPeH1c8lwAYaX7z9PeDYac7dDBwIHAw8ANwF3A/8Toe5h4HPAmcAP9fO+VXHc9xrxvRfC3qdo5xzt61zre7HfEHnRin/ZrV1/Hm8jM31a98WtI3b2i79Gke3681MZjGgIuKMCWa9GvjksOaAs4A3AbsCb42Ij2Xm30bEWS3PKyV3JPB8YBnwfhr/AXkC+FAHuV4uayq5+4CHgKgeP4/GCT8fOLsGuRU0LvanAydm5vqI2B14L/B305j7amb+a0RcTOMLap+IiF1obMu/6CD3zWrdvwH8VUR8A/hAZj7UlKn9Oe41Y0auBb3OlXLutnWu1f2YL+jcKOXfrLaOP4+XcXN92bcFbeO2tku/xjGF5U1oYAs9GheEr/DTg3HUnkOe+5fM/Ofq97si4vKIeArIQnP3ZeYW4MsRcWtmPg4QEfd3kOvlsqaSOwb4n8DqzHwyIt6emZdFxB/VJDdqa1ZvZ8jMxyNiwzTnvlf93JSZT1S5H0fEwx3mMht/CvsM8JmIOAi4KCKelZm/05Sr+znuNWP6rwW9zpVy7rZ7rtX9mC/l3Cjl36xRkx1/Hi9jc/3at6Vs43a3S7/G0e7yJpZtvvRXtwl4JY2//LS2nzXkuQ8Cz2ppex/wnUJzNwLLW9qeA1wz1VwvlzWVXNU+AlwNHAK8rWqbVYccjYvMEcBHgb2qtgXAddOc+x3gAuDCKns68C7gyg5zXwLmjrMNFk/zOVn3XN3P8Rm/FkxTroRzt91zre7HfCnnRhH/Zk3h+PN4acn1a98WtI3b3X79Gkdby9vZ1FaorhMwz9yYzKnA+S1tuzD2H+JSckcCr25pezPwxqnmermsqeSa5gVwPnD9JPt4RnPVOEanRVXb24Hfms5c1X44jbdQfBb4Gxr/sZzdSQ64CDhusnOoytb2HO91bgDO8Rm/FkxHrmof6HO3am/3nKztMd9uzmN+Zv7NmuLx5/HSlOvnvi1kG0/l+j3j45jKeieavOumVDMR8ezMfLquOUnd8dxVSTxOy+U+G3wWepIkSZJUmF363QFJkiRJUm9Z6EmSJElSYYor9CLijIhYaM5cr3J17ps5c+YGN1fnvpkzZ25wc3Xum7mZy0GBhR6wP/CpiFhgzlyPcnXumzlz5gY3V+e+mTNnbnBzde6buZnLlXkzlqrKfTozv2/OXC9yde6bOXPmBjdX576ZM2ducHN17pu5mct1/L0MTk5OTqVOwO40fXcNsAw4wpy5TnN17ts05Q4Dfr76/deBM4C55sx1k2t5zqnVz1eb6+s69wFOAX5lkmUNbA7Yc/R4pPHdexcCvz0DuQuBF+ys35NNHT+xrhOwypy5Xubq3Ddzvc8BbwIeBx4F3lK1HQjcZ85cJ7k6922ach8A7gW+ReML1W8HrgWuNWeui9wZLdOnq583DWuuj337P9XPw4ANwF8DXwbeXmjuNmApcCJwB/Au4K+AK6Y592/AR4C/YJJCdaJpFgMqIm6YYNaLgSvMmZtqrs59MzdzOeBlwGIggDMjYnVmXhoRj7c8z5y5dnN17tt05H6YmQdHxHOBWzLzKICIuNScuS5yq2n8J3xL9fgFwCE0Xh0Z1ly/+jb6+FjgpZm5HSAi3l9o7o7M/E5EnAIcnVUVFhEXT3NuU2aeHhGLgTMi4nwaxegnMvNHtGFgCz1gNnArsKOlvfUGM+bMtZurc9/MzVzu65n5TPX7VRFxUkS8BUhz5jrM1blv05F7FCAzvx8RH2pq32rOXBe5lwIXAx/JzC9HxNsz87KIeMMQ5/rVt1FbR4ujysOF5qL6ua36ffSa11pH9TqXAJm5CfjjiHgW8HrgJuD4CcbSsoQOXgaswwTsC7x+nPbXmDPXSa7OfTM3o7nrgeUtbW8CnjBnrpNcnfs2TblbgNNa2g6h8Vdoc+Y6ylXtAZwHnAW8rXX+MOb6tM4HaHx+7B+orgnAq4FPFpo7FPgUjbcUrwf+ElgLvHmac/8KHDTRvmpn6viJTk5OTiVOwHLghHHaWwtCc+baytW5b9OUGwGe39J2BPAyc+Y6zbXMPwp470TzhzHXr7415V8BHFhqDpgH/BbwNuDNwP4TPL9nORo3anlRu/tg3PV08+R+TTTearVwJ/MDWGjO3BRyi2rcN3PmzA1uzmuLOXPmvLaYm5bcRPNHp9bPpgyEzPwxcEFEvLx1XkT8PPBhYJY5c1PI7VrjvpkzZ25wc15bzJkz57XF3LTkWueNMVklWNcJmEvjlqP/AXwN+ArwIPBPwEvMmZtqrs59M2fO3ODm6tw3c+bMDW6uzn0zN3O5nU1RLWhgRcRzaNxk4dnAQ5n5iDlz3eTq3Ddz5swNbq7OfTNnztzg5urcN3Mzlxv3uYNe6EmSJEmSftZAfkZPkiRJkjQxCz1JkiRJKoyFniRJkiQVxkJPkiRJkgpjoSdJkiRJhZn8i/YkSRpCEbE/sBpYDxwI/ClwMo1/O38EbMvMd1fZPwNmAw8DxwLfAv4gM7f2oeuSJFnoSZI0gWOBp4D3AnsCRwC/nJm/BhARd0TEbdW8ZZl5bNX+68C1FnmSpH6y0JMkaXwfAc4D/hH4BvA14Oci4rxq/kPACPCLNF7BG3X/THZSkqTx+Bk9SZLGdyiwJjMPBR4FAtiUmWsycw3wcRoF4Hpgv6bn7TPjPZUkqYWv6EmSNL7dgfdExP00Xrn7IPCsiHgX8Awwh8Yrfg8Bx0bER4Bv03i7pyRJfRWZ2e8+SJJUjIi4lsZn9O7oc1ckSUPMt25KktQjEfFy4JeAUyJiXr/7I0kaXr6iJ0mSJEmF8RU9SZIkSSqMhZ4kSZIkFcZCT5IkSZIKY6EnSZIkSYWx0JMkSZKkwvx/LhNk0vDUnFkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(15, 5)) \n",
    "pd.DataFrame({\n",
    "    'seq': sorted(combination_counts),\n",
    "    'n': [combination_counts[k] for k in sorted(combination_counts)]\n",
    "}).plot.bar(x='seq', y='n', ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that we have all possible integer sequences in our sample\n",
    "\n",
    "To do this we ideally want an analytical solution to the number of positive integer sequences of $N$ elements that sum to $M$. This can be seen as the problem of taking a set of $M$ elements and considering how many ways of partitioning those elements into $N$ sets. To split the $M$ elements into $N$ sets we need to determine $N-1$ division locations consequently which can be placed in any of the $M-1$ gaps between elements, consequently there are \n",
    "\n",
    "$${M - 1} \\choose {N-1} $$\n",
    "\n",
    "ways of forming these partitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import comb\n",
    "\n",
    "assert comb(m-1, n-1, exact=True) == len(combination_counts)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:pseudo-threads]",
   "language": "python",
   "name": "conda-env-pseudo-threads-py"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}