{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "atmospheric-upper",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import nest_asyncio\n",
    "nest_asyncio.apply()\n",
    "import stan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "warming-fields",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mixture parameters\n",
    "\n",
    "N=1000\n",
    "\n",
    "K=4\n",
    "p = np.array([0.25,0.25,0.25,0.25])\n",
    "\n",
    "mu= np.array([1,3,5,7])\n",
    "\n",
    "sigma = np.array([0.5,0.5,0.5,0.5])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "academic-funds",
   "metadata": {},
   "outputs": [],
   "source": [
    "#sample the mixture to create artificial data\n",
    "\n",
    "np.random.seed(0)\n",
    "\n",
    "Zt = np.random.randint(0,K,N)\n",
    "\n",
    "X=np.zeros(N)\n",
    "\n",
    "for k in range(K):\n",
    "    pts=np.where(Zt==k)[0]\n",
    "    npts=len(pts)\n",
    "    X[pts]=np.random.normal(mu[k],sigma[k],npts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "lined-giving",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prior parameters\n",
    "mu0 = np.array([0.,0.,0.,0.])\n",
    "sigma0 = np.array([10.,10.,10.,10.])\n",
    "\n",
    "c = np.array([1.,1.,1.,1.])\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "innovative-sunset",
   "metadata": {},
   "outputs": [],
   "source": [
    "import stan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "internal-malawi",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Building...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Building: found in cache, done.Messages from stanc:\n",
      "Warning in '/tmp/httpstan_4zz45_p4/model_6pyh3mdl.stan', line 6, column 2: Declaration\n",
      "    of arrays by placing brackets after a variable name is deprecated and\n",
      "    will be removed in Stan 2.32.0. Instead use the array keyword before the\n",
      "    type. This can be changed automatically using the auto-format flag to\n",
      "    stanc\n",
      "Warning: The parameter mu has no priors. This means either no prior is\n",
      "    provided, or the prior(s) depend on data variables. In the later case,\n",
      "    this may be a false positive.\n",
      "Sampling:   0%\n",
      "Sampling:   0% (1/8000)\n",
      "Sampling:   0% (2/8000)\n",
      "Sampling:   0% (3/8000)\n",
      "Sampling:   0% (4/8000)\n",
      "Sampling:   3% (203/8000)\n",
      "Sampling:   5% (402/8000)\n",
      "Sampling:   9% (702/8000)\n",
      "Sampling:  11% (901/8000)\n",
      "Sampling:  14% (1101/8000)\n",
      "Sampling:  15% (1200/8000)\n",
      "Sampling:  19% (1500/8000)\n",
      "Sampling:  22% (1800/8000)\n",
      "Sampling:  25% (2000/8000)\n",
      "Sampling:  28% (2200/8000)\n",
      "Sampling:  42% (3400/8000)\n",
      "Sampling:  61% (4900/8000)\n",
      "Sampling:  79% (6300/8000)\n",
      "Sampling: 100% (8000/8000)\n",
      "Sampling: 100% (8000/8000), done.\n",
      "Messages received during sampling:\n",
      "  Gradient evaluation took 0.000273 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 2.73 seconds.\n",
      "  Adjust your expectations accordingly!\n",
      "  Gradient evaluation took 0.000245 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 2.45 seconds.\n",
      "  Adjust your expectations accordingly!\n",
      "  Gradient evaluation took 0.000381 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 3.81 seconds.\n",
      "  Adjust your expectations accordingly!\n",
      "  Gradient evaluation took 0.00049 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 4.9 seconds.\n",
      "  Adjust your expectations accordingly!\n"
     ]
    }
   ],
   "source": [
    "cluster_code = \"\"\"\n",
    "\n",
    "data {\n",
    "  int<lower=0> N;         // number of points\n",
    "  int<lower=0> K;         // number of components\n",
    "  real X[N];              \n",
    "  vector [K] sigma;\n",
    "  vector [K] sigma0;\n",
    "  vector [K] mu0;\n",
    "  vector [K] c;\n",
    "}\n",
    "\n",
    "parameters {\n",
    "  vector[K] mu;                // population treatment effect\n",
    "  simplex [K] p;\n",
    "}\n",
    "\n",
    "model {\n",
    " \n",
    "  for (k in 1:K) {     \n",
    "     target += log(p[k]) + normal_lpdf(mu[k] | mu0[k], sigma0[k]);       // prior log-density\n",
    "  }\n",
    "\n",
    "  for (k in 1:K) {     \n",
    "     target += dirichlet_lpdf(p | c);    // prior log-density\n",
    "  }\n",
    "\n",
    "  \n",
    "  \n",
    "  for (n in 1:N) {\n",
    "      vector[K] lps;\n",
    "      for (k in 1:K) {\n",
    "          lps[k] = normal_lpdf(X[n] | mu[k], sigma[k]);\n",
    "      }\n",
    "      target += log_sum_exp(lps);\n",
    "  }\n",
    "}\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "cluster_data = {\"N\": N,\"K\": 4, \"X\": X, \"mu0\": mu0, \"sigma0\": sigma0, \"sigma\": sigma, \"c\": c}\n",
    "\n",
    "\n",
    "posterior = stan.build(cluster_code, data=cluster_data, random_seed=1)\n",
    "\n",
    "fit = posterior.sample(num_chains=4, num_samples=1000)\n",
    "\n",
    "mu = fit[\"mu\"] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "threatened-guarantee",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.94646088 2.93110097 5.01200022 7.01476528]\n",
      "[0.25 0.25 0.25 0.25]\n",
      "[[0.97063669 0.89666359 0.91029708 ... 0.9866595  0.98852443 0.91273853]\n",
      " [2.96792948 2.87518591 2.88268633 ... 2.93130679 2.90934216 2.91433174]\n",
      " [5.02782928 4.99904749 4.93856443 ... 4.99127792 4.96533382 4.99514043]\n",
      " [7.03952579 6.98701906 7.01649147 ... 6.94257406 7.0003452  7.04601038]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0/UlEQVR4nO3dd3gVVfrA8e8hgYQaiIReQgcNghBAVAQLKKCI2HBV1FX5rWVtK2LHLqsrYm+oYAVBkSKogKAgzYQOoSa0AEkgvd72/v6Yy00uuQkJpMyu7+d58jhz5sw575w582bunQkaEUEppZR91ajuAJRSSpVOE7VSStmcJmqllLI5TdRKKWVzmqiVUsrmgiuj0caNG0tkZGRlNK2UUv+TYmNjj4pIRKBtlZKoIyMjiYmJqYymlVLqf5IxZl9J2/SrD6WUsjlN1EopZXOaqJVSyuY0USullM1polZKKZvTRK2UUjaniVoppWzO1ok6e/kKHAcOVHcYSilVrWydqA/cdRd7Bg8pc31xucjfti3wNocDx4EDJL/5JnFdu5G/c2exOp6cHFwpKeWKURwOxOMp1z6VqaL/fXFPQUGp293p6bjS0nzr+Tt24Dx8uFg9x4EDZP7yiy9GZ2IiIoK43WWKI2/zZjx5eYFjyMrCk59fvDwjg6RXX0McjhLbLdi1C3dWVpliKC+Pw8G+W28jbfqMUuslvfoaadOn49gX+O8dXGlpZP/2G+J24zx0qEx9iwgFu3fjcTh8x+/Jz8ednVPiPs7ERHLWri1T+6X1m/z6JPI2bwm4Pf372WT8+ONp9VHRPAUFxeZW/o6dOBMTqymi4k76l4nGmC5A0ZnWHnhGRCZXVlDgn3Dc6elkL19O2JVX4kpLI+Wtt/BkZRPUoAFpX38NQJ1+/QiOiCBz/nzazfmB4PBwatStS9aSXzk0blyx9o88/zyt3niDXQMupMGwoTR94gl2XTAAgLCrr6beRYOoe9751AgNgaAgUqdNw5OTQ8PrrqNmkya+draf3cO33PTxx0h6ZSL1hwyhRv16uJJTyFm+nMjvZhHauTPujAychw6R9NLLOJOTcR0+TKcVy8nbsoV6550HNWuS/euvpH/3PU3HP0rNli0xwcG40tKoERKCOJ0k/ftVGo2+geAmTTC1auFKSuLoBx9SKzKSvPXryV27lg4//0TOqlVk/vwzLSZOZP+YW6k/ZDD1LrwQd1Y2B++5B4AuGzeQPuNbsn75hcb3/5Madety7IMPcSUnEzljOilvvc3R996j/cIFZM7/kezffydyxnRchw9Ts2VLAHae2x+AZi88T8OrriLhqpEANLxxNHX79qXB0KEAvl+49TZt5MgzE8j44QffuHXdugVxOqkRGkpBfDzxw4bTacVyMubMIfm1/xDctCmupCQAmk98hYYjR1pJvqCAGqGh7OzT19dnzWbN8WRnk/nLzzj37QegVmRbGl1/PY6DiWTOn88Z/zcWYwy5sbHsu+lmALptjyNp4r9JnTqVVu+9R/2LLyLl7XcI6dqFBoMHkzZ9Brlr19Lo5pvZ97e/AdDx1yXk79hBSPv21lxbvJiMOXNpMHw4SS++6Du+3DVraDjqakytWhQkJBBUrx4H7/sneRs3FpuXLSe/QYPLL/et5+/YyaFx4yjYuZOG111H+syZtHjtVRpccQXutDQce/dSp1cvjn3yKRnz51MQF0fYqFGIy0nm3HkABIWH03Hpr+zoeQ4AZ4wdS5OHHyJn5UpcaWmEDR8OwO5LLvX122HxYtxpaaR9+SUmJAQpyOeM//sHtSLbIi4XpmZNUia9gacgn8w5c+m0aiVpX3/jO+5jH39Ms+efw3kwkbARVxJ/xZXUatcOR0ICgDUW23dw9N13CT3zTBrdfDMhHTvgOHCAsOHDyV23joy5c2k2YQK5a9ZSu3sUntxcgsLDyY2NpVZkJEFhYaTPnEXSiy/SfsECHPF7qBMdDUFB5G+Lo07fPhhjyNu8BXd6Gu60NIKbNMWZeJDa5/QCA+kzZ5H66afeY16EOByYoCASrrrKNy8ACuLjqVG7NsER1l94p339DSEdOxAUHk7e+vU0HD0aY0yx81lRTHnuwIwxQUAi0E9ESvxzx+joaDmVPyHPXLSII88+h/vYsYDba7ZqhfPgwXK3WxlCo6IoiI9HcnOrOxRVAZo++SRJL73kW6/ZokWZ717LqkadOnjKMF9q9+yJqR2KY98+XIeKfzqpDDXq1sWTU/Ld9nG1e/Uib906Gt5wA+kzCu/f6g8eTNaiRZUZ4n+FOv360Xba1FPa1xgTKyLRAbeVM1EPASaIyPml1TvVRB3XtVu591FKKTs5fhdeXqUl6vJ+Rz0a+OaUolBKKXVKypyojTG1gBHAzBK2jzXGxBhjYlLK+UBOKaVUycpzRz0UWCciSYE2ishHIhItItEREQH/SdWTavLoowB0+mMFXTdtJPTsswFo/tJLNH3i8VL3rdOvH13jttF+3lzazf6e+kMvp+nTTwHWg4uigsLDaTJuHE0ff8xXVm/gQN9y7V69/Oq3+PdEmj7z9EnjbzdnDm2mTaPhDTf4lYd07kyHn3+izefT6BwTQ7u5c2j71ZdEzphOsxeeByDiwQdpMm4cYd6HGM2enUCrD973tdF6yhTq9Onj127n1asCxtFw9A2EnnUW4beOoWvcNkI6dQpYr95FFxF29dU0GHGlr+yMu+6k/YIFNHv2Wf9j6NSx1GNv+uSTgcufeJzIGdMxISG+srCrRtBlXSyd/lhBm8+nEfHgA7T57FPf9to9evi1UbNFC7/1+pddhgkN9SurfY71oKzlm2/SZf062kydCoCpXZs2U6fS9MknaT9vLkHh4b59Go25haBGjfzaaXz/P6nT/1y6bY+jds+eAY+p3Zw5NH/pRWqEhdHgyitp9f57tP74I786rd55u9h+zV96kRb/nkjYNaMCtgv+H5ubPPIv2n7zdcB6XeO20ezZCcXOE0D4rbfSZNwjJfYR0q3wK8Zmzz9H/cGDab/gR1p/+EFh348+St3z+tNy0ut0WrWSWh060G7OD77tLd96k4iHHvKt1yry78+3njKFlpNeJ6RrV2p17ED7BT/S9OmnqDvwQoKbN6dJgIf7x9UICyNy5reFx3L77b68UKNePWr37k2NsDDf9rCRI4l46CEa/e1GTGgobb/4nJAzuxHUuDENhg2j5aTXaTBsKF02baTb9jga3XQTLSdPJuzqq6l73nl+fTe+/5+E3/F3v7Kghg3ptj2OyJnW/WnNtm0Iu2pEifFXGhEp0w8wHbi9LHV79+4tFSF38xaJi+oujiNJIiKSt2OH5O/cKQX79omIiDs3V3YOuFAylyw5aVvHvvhStnXpKolPPFFsmzMlRUREMhctkvzdu62y1FTZPeQyyVq+ImB9j8MhBQcOSM66dXLgn/eLx+n0bc/6fbls69JVjn76mXhcrvIfeAnc+fniTE0VR1KSuPPzRUTE4/FI/q5dcuippyV75UrJi4srtp/H7ZbMxYvF43BI5pJfJXvVKsnduLFMfebHxxf25XbL0SlTxJWZKXlbt0ruhg2S8+efknDTTeIpKBBPQYE4Dh2y6jockrl0qXg8Hl9browMceflldjX9l69ZVuXriIikj5njqR++62kz5kjztRU8bhc4khM9Kuft327ZP/xhziSksp0LL5j2r1bHIcP+44p+c23ZFvXbuLKyvar50hMlIPjxkn+nj3iTE0VZ2pqqe06k5PFlZHhV7Zz4CDZ1qWrb1yOy161WlLee0/yduyQbV26yr4775Ls1WtERCRz8WLJ27atMI5DhyQ/Pl4chw/L7uHDJWfduoD9p3z0kd+2zMWLZVuXrpI+b74cePBBOfz8C5Iw+sZS52T6/PlSsH9/idu3dekq8aOu8TvmgoQEX3/5u3aVuG9R7pwcyd+zRzwulyQ++aRs69JV4nqe49fP8blQWTwej6TPnSe5m7dI4qOP+q5hj9stye++K1nLV/hyg4iIu6BAPG63eDwe8TidxeaDOztb3Nn+c6g8gBgpIaeW6WGiMaYusB9oLyIZJ6t/qg8TK1vmL79Q78ILqXHC3VhlyI2NpfY552Bq2PpVdVtxHEzEkRBPvQEDqjuUKlWwZw+1IiMxQUHVHcpJicsFxlR4rK60NIKLfLrJWb2aGnXqUNv7qfqvoMLe+igruyZqpZSyq4p860MppVQV00StlFI2p4laKaVsThO1UkrZnCZqpZSyOU3USillc5qolVLK5jRRK6WUzWmiVkopm9NErZRSNqeJWimlbE4TtVJK2ZwmaqWUsjlN1EopZXOaqJVSyuY0USullM1polZKKZvTRK2UUjaniVoppWxOE7VSStmcJmqllLI5TdRKKWVzmqiVUsrmypSojTENjTGzjDHbjTFxxpj+lR2YUkopS3AZ670J/CQi1xpjagF1KjEmpZRSRZw0URtjwoALgdsARMQBOCo3LKWUUseV5auPdkAK8JkxZr0xZooxpu6JlYwxY40xMcaYmJSUlAoPVCml/qrKkqiDgV7A+yJyDpADPHZiJRH5SESiRSQ6IiKigsNUSqm/rrIk6oPAQRFZ412fhZW4lVJKVYGTJmoROQIcMMZ08RZdAmyr1KiUUkr5lPWtj38CX3nf+IgHbq+8kJRSShVVpkQtIhuA6MoNRSmlVCD6l4lKKWVzmqiVUsrmNFErpZTNaaJWSimb00StlFI2p4laKaVsThO1UkrZnCZqpZSyOU3USillc5qolVLK5jRRK6WUzWmiVkopm9NErZRSNqeJWimlbE4TtVJK2ZwmaqWUsjlN1EopZXOaqJVSyuY0USullM1polZKKZv7Syfq1PzUau3/sy2f8fef/16tMdhFQkYC18+7nkxHJnmuPJwep2+by+MiPiO+GqMrzu1x82P8j3jEU2KdDckbyHZkV2FU/9tWHVpFWn5adYdRLWyfqOOOxbE+eX2xcpfHRXp+OkdyjgCQlp9GbFIsSTlJdJ/WndHzRwdsb8vRLWQUZLA5ZTMDZwykx+c92Hpsq2/7syufZfr26SXGk+vMxeF2FCu//afbeWvdW7y85mWcbic703ayL3MfR3KOcDj7cMC2JsVO4s8jf5KUk8SG5A28sOoFYpNiuXvx3axIXMHGlI0ALD+4nAHTBzBj+wx+O/AbIlLygJVi+PfDuXfJvexO2827G95FRPCIh0xHJun56XjEw5TNU9iXuY8ViSsCtnEo+xDf7vgWgCeWP0H3ad3JKMjwqzNj+wxWHlrpV+b0OMl0ZPrWlx1YxqHsQ7g9bmbvms2YhWOIS41j+cHl9P2qL0NmDQEg25HNm+ve5KofrmJ32m5fYvSIh/iMeL++XR6X79zkufLId+Xj9rj588ifZR6j9ze+z897f2ZS7CSyHFl+23478Bvvb3wfgHG/j+Ox5Y/x3a7vAOsXzexds9mRuoMDmQfIceZwy8Jb6P9NfzIdmcXmjNPtxOF28NqfrzF712zfPBYRXzLakLyB3Wm7fcc2ZNYQluxfQq4zl3fWv0P3ad2Zt2ceR/OO0n1ad7pP68721O2+dibHTmZX2i7ffHG4HXyw8QO2Ht3qF0uBu4A8V17A8fCIh8+3fk7MkZgSx0xEAs7Jn/b+xOxds0vdb0/6HnKcOcXm0ImWH1zO2EVjuXDGhXSf1t3vmnJ6nPy6/1e/GLanbifuWJxv3e1x43Q7WX5wOf9Y9A9EBKfbicvjKrVfsK75GdtnlHjdZToymbJ5Sqm/tE+XKctFb4zZC2QBbsAlItGl1Y+OjpaYmJJPbEk84uGb7d9wKPsQl0dezqxds/h+1/cAvHD+C4zsOBIRYdC3g/zuhmNvjuXaedeSkJHg196DvR6kce3GPPXHUyfte+GohUzdOpUZO2YAMLLjSOrWrMslbS7h400fM7z9cAC/tj4f+jnN6zbn7fVvM3fP3JP2MarTKN/xlNXrA1/nX7/9q1j5+S3OZ0SHETg8DjambGRP+h7WJ69nWLthtA9rT/N6zXlr3Vsk5SaV2PaHgz9k2tZpxZLqcQ/3fphJsZMAmD58OqN/LPzld+uZtzJt2zTf+tkRZ7MpZRP9mvVjzZE1vvJW9VpxMPugb33JdUtoFNqIXl/0Cthno5BGpBWU/64pskEkaQVpZBRk8PbFb/PPX/8JwLWdr2XWzllc0+kaekT04Ptd3/NY38f4ZMsnNAppxNB2Q5m1axb39byP5YnLeXnNy37t/nLNL9StVZd31r/DN9u/AeD5857nmZXP+OrMvHIm1827zm+/4+NRdN3hdtAurB0vX/Ay53xxTrFj2DRmE2+vf5uPN3/M/Kvnc8XsK8o9DiX5YugX/Lz3Z76M+xKAK9pfwaN9HsXlcXHxzIt99cJDw5k+fDr7s/ZTK6gWYxaO8W2bPGgy3SO6YzAczjlMeGg4k2InsWjfImoH12btTWsREY7lH+No3lHfmLx0wUu0rNeSmCMxvLPhHV97Y88ey0ebPioW65yRc9h6dCubUjZxcZuLuWfJPcUS6siOI7m7x93Mj5/PT3t/YlfaLr9jvWXhLYA13wyGG3+8scRrYd7IeUSGRfLKmlf4evvXxbYH1wjG5XFxX8/7OPOMM7lnyT10btSZUZ1G0bNJT55Y/gTxGfEMjRzKywNeJrhG8EnPRyDGmNiScmt5EnW0iBwtS4enmqi7T+te7n2UUsouBrUexNsXv31K+5aWqG3/1YdSSv23WHZgWaW0W9ZELcAvxphYY8zYQBWMMWONMTHGmJiUlJSKi1Appf7iypqoLxCRXsBQ4F5jzIUnVhCRj0QkWkSiIyIiKjRIpZQ60d+j/jpvTJUpUYtIove/ycBsoG9lBtUtvBst67Vk2fXLALgj6g7eu+Q9ABZcvYA3Br1B0zpNWThqIVe2v5I1f1tD32ZWSDd3u7lYew/0eoAJ/SfQsWFHwHoQUdTXw4o/QDhuaLuh3NX9LhaMWlCsze9GfEfMzYXfxb950Ztc3/l6BrcdzEeDiz8kARjfZzxBJoiJAyaWPggleKLfEzx97tN8MfQL7ux+JzVr1OT5855nYKuBfvWubH8lsTfH8s7F1sObDy/9sFhbm8ZsYsXoFQyNHOor++PGP4rV+2rYV77lmVfODBjX9Cums/nWzWy+dTNLrltSbHu38G58etmnxcrvP+d+Hu/7OPNGziPm5himDJnCA70e8KtzfovzGRc9jn7N+gXsu6jjMbx7ybvWvi3PZ/aI2Xwy5BPu6XGPr154aDgNQxoCUDu4tq/8pQteYuzZ1ofGDmEd6NiwIx9e+iEfD/mY32/4nRWjVxQ7juPjd1Hri/zK69as67c++aLJrBi9gjkj5/jKHuz1IAB1gusA1vzdfOvmYufr5QteZvOtm/l+xPe0rNeSLo26+LatGL2CkKAQ33r/5v358NIPCQ0K9WujSZ0mvuVRnUZxc7ebWThqIY/3fZyXL3iZ32/4nR+u+oFAzm1+Lo/1fQywHhBPvmgys0fMZtaVs5g0aJJf3Xkj5wHQpVEXll2/jJ+v+Znw0PBibf5n4H94ql/hw/nIBpHUDq7N0uuXFqsbaP+Hej/E5ls3+9b7NevHFe39H8D+dsNvrPnbGh7r+5hvrAe0HEBoUCjRTaP5athXxea0wdCibgsmDZrE2pvW+sUzoOUAJvSfwMYxG1ly3RI2jdnEle2vBODF819k3c3rig9eRTj+ak1JP0BdoH6R5ZXA5aXt07t3bzkVM7bPkKipUZJVkHVK+x+Xnp8uKw6ukL5f9pWoqVGS68wVEZEdqTuk/1f95Uj2EUlIT5CoqVEybtk4ERH5fuf3kpCeILtSd0nU1Ci54vsrirUbNTVK7vj5jmLlW1K2yL6MfQFjyXHkyLc7vhWn2xmwvaipUbL16FbJKsiSA5kHJGpqlMzdPVcyCjLE4/HI0dyj0vuL3vLVtq/E4XKczrDInvQ9smz/Mlm0d5H8uu9XX3m+K1/m7p4rx/KOiYiIx+MRp9sp8enxAeN2e9zy8aaPJWpqlLyy5pWAff2450f57cBvIiKSkpsiHo9HRKxx3n5su2QWZMrENRMlz5kXcP/Hfn9MoqZGycKEheJyu3zlV3x/hfzfov8TEZH1SeslamqUXDX7Knlk2SOy9vBavzZWH1otDrf/mH22+TOZumVqsf6ipkbJxDUTfetHso9IgasgYGwiIg63w297Wl6ab/m3A7/JwviFJe4rIvLIskckamqUiIhkO7LF4XLI/sz9fuOdkpsiqw+tLrGNpJwkyXHk+NbXHl4rh7MP+x1zfHq8xB2LExERl9sld/58p6w6tKrU2I47lHVIoqZGydVzrvadv0A8Ho8s27/MN59L8uHGD6XPl31k3LJxfsflcrsk25HtVzchPUHG/jJW0vPTRURkzaE1vvbn75kvKbkppcY+OXayrEtaV5bDlD1peyRqapTcvejuEuvkOfNk3p55pY7D6QJipIScetK3Powx7bHuogGCga9F5KXS9jnVtz4qWkpuConZifRs0rNc+21I3kD7hu1pUKuBX3m2I5uQoBBqBtWskPhG/jCSIZFDuKfnPSevrP6niPcd9qAaQdUdSoU5lH2I5Nzkcl9vZZXpyGR/5n6iGkdVeNs/7P6Bi1pfRFhIWIW3XVan/XpeedklUSul1H8LfT1PKaX+i2miVkopm9NErZRSNqeJWimlbE4TtVJK2ZwmaqWUsjlN1EopZXOaqJVSyuY0USullM1polZKKZvTRK2UUjaniVoppWxOE7VSStmcJmqllLI5TdRKKWVzmqiVUsrmNFErpZTNaaJWSimb00StlFI2p4laKaVsThO1UkrZnCZqpZSyuTInamNMkDFmvTFmfmUGpJRSyl957qgfAOIqKxCllFKBlSlRG2NaAcOBKZUbjlJKqROV9Y56MvAo4Km8UJRSSgVy0kRtjLkCSBaR2JPUG2uMiTHGxKSkpFRYgEop9VdXljvq84ERxpi9wHTgYmPMlydWEpGPRCRaRKIjIiIqOEyllPrrOmmiFpHHRaSViEQCo4FfReTmSo9MKaUUoO9RK6WU7QWXp7KILAOWVUokSimlAtI7aqWUsjlN1EopZXOaqJVSyuY0USullM1polZKKZvTRK2UUjaniVoppWxOE7VSStmcJmqllLI5TdRKKWVzmqiVUsrmNFErpZTNaaJWSimb00StlFI2p4laKaVsThO1UkrZnCZqpZSyOU3USillc5qolVLK5jRRK6WUzWmiVkopm9NErZRSNqeJWimlbO6kidoYE2qMWWuM2WiM2WqMea4qAqt2jlxw5hWui1RfLOrU7FoM8b+Vvb4IeDyn329uKrhdp9/OX4EILHgUkrZWdyS2VpY76gLgYhHpAfQELjfGnFupUXncsGgCZCcXlh34E1yOwvW9f4DbaS3HfApp+07ebtx82L+65O07f4HtP8LSV+Dl5jCxrVUuAs81hJ+fLHlfZx6s/RgSY4v34XZZx1QWeWkw70H/XxIlObYHspL8y3KOWvvnpVlJZ/MsSNlx8rYcOf7rm2ZaY7r0ZVjzkf/YA+Qc819f/T5smwuuAvjxEchO8d/udkJeOuRnwh9vWm2eeIy5qZDw+8ljLUrEv53MQ7B9gbX81TXw+Qir3aO7rboF2fDpUEjaVrytz4bBC429x3e0fHEAHFgLydvh1XYw977S51ppx5OfaS1v/cE6L6veKzwmsOZSYmzp7eRnQuZha7+fnrCO21UA+RmF/ZT35iM1wfo50ZoPrXkWiMdjxeIqsNbz0guX45fB3uWw9kP4YlSAfd2ln4e4+dY1eTDGWj+8qfB6KMux7fy5MBaASWdZc+A4R67Vfm6qNa+qkZFynCxjTB1gBXC3iKwpqV50dLTExMScelTxy+Dzq6DzUBj5Hmz4Cn55yto28n3rgoj9DFqcAyPegQ/O99+/QUsY/ZW13VUAc++HGsGw4Utre5fhsONHa3nQExDeDuo3g2lXli2++9dbyTDhN7hkArTpD59dHrjupc/B4gnQMhruWmKVPRtWuH3k+7DybagRBEc2F5ZfMgGWeD+83PsnbPkOLngI9q+0Er8zF2be6m0vwxozjxu+9E74xl2g4yWw+r3CNtuebyXMO36xLoDFz0LfuyDmE1j3OdSNgNsWwPZ5sOR5/+M462oYNQWCgmH3YvjyGuh5k3Vu/rUDXu9i1bv6I5g91ntsH1jxDp9UmABP9HAcfHsrNIuyjjE/A4b9BxY8Ym1/JtUam4JsyEyEMzrBzDEQN8/a3rAtpO+zxiiiM0zuDun7oceNsPEb/77qRkCO9xdIrXrw+EEwxroJWPAIHN5gbWvdDw4Umd49boSjuyCxyJyOuhZc+TBwPKx4A7Z+H/j4Og6GPndYF3uPGyFpMwSFQK06Vhlinf+Bj8Gqd2DdNGu/y/8NP42HM6+CbXMK27vKez3s+wO6XWldI92vg+BaVlJM3gaNO8OLEf5xBNcGl/cX2t9/gU+HeI/jGhj+OpggmP+gdQ7a9IdLn4WQ+rDkBdi50L+tTkOs4zpzhBXHrL9b5bfMhtrh1nrthtCwDWydbW0LCgF3kaQ4bg+81qFwvXY4jE+wroGtP8A5N8OfU6wx6XEj9PsHtOhp/SKfdiXUOQNyT7hZOO742N37J+xeBNlJ1nVoDBRkQVAtawznP2TVv3cthLW2bs4Aet0Kg5+D2f+AnT8VtnvOzdYYb/wa2g2Elr3hyCZrTl46AWrUtK6PU2SMiRWR6IDbypKojTFBQCzQEXhXRMYHqDMWGAvQpk2b3vv2leEO90RvnQOp8YAB/ge/amjeAw5vrO4olFKVZcxcaD/wlHYtLVGX6WGiiLhFpCfQCuhrjIkKUOcjEYkWkeiIiIhibZRJavzx1k5tf7vTJK3U/7bPR1RKs+V660NE0oGlQAmf85VS6i+uLM+Xyqksb31EGGMaepdrA4OB7RUeiVL/c0x1B6CqWkgDqFm7wpstyx11c2CpMWYT8CewSETmV3gkJ7p8Ilw3DW7/yXpY9fRROKvIk+FLJlgPWgK5bYH1oOj8B4tvG/GOte24qGvhislQr2ngtppGWQ9iiu7z2AFKvQjvWGw9ADnRowkwfi+Miy++rSS9vA8M214Anb0fZK75BGrVh0GPQ7cR1kOl6Dush3anIzgUxszxLzvnZv/1yAHQ42+F6zXrwtmjYdTH/vUe2299X/dkEozfB/eV8HA5qBZc9W7hetS1/ttv/g4ue6V4Wf/7/Mu6jYBbyzEt74uBmnWs5Ra9/LeN31u8fngH6DsWHtkNYW2ssnPvtfrtcAlc/wVcdMJbQc+mw+0LrQdnY+ZaZV2vsB5cnejJJBj6auH6sP9Y537Q41C/hVXWuh/cvcpaPutqqy2APndZY9L7dqtOWYyZY11Xx+vfsahs+x1nSkgd9Zpax9otwFcAbc4LvM8zaYHLr5tWuDxmLjy0DZ48cvLY7lgET6XAmSMLy4JDoW4TCG0ID24pLH8qxXqYWNQ1nxQu3/KD/7bIAdbD+1FTrPWic3NcPDxajmu7PESkwn969+4tpyT2c5EJDayfkuSmiWQlF657PCL5mdY+nw4Vceb713cWiGQeLt5OwnKRlJ3F257QQGTqlSIup7VeVF66SEaitex2ibgcVvvpB0Xm3i+SsMKK5bgJDUQWjBeJ/10k/rfiMeRniSTvsPbPPmq1lbLTamvHT1b75fVSC6vf7BRvnG6RmKkiR3eLrP9a5KcnrPKV74psmyey5XurTvzvIllJ1rZ9q0R2LhJJ2WXFsPwNq83Pr7a2pyZY61MG+4/3b6+KzLkv8HiLiGz+rvD8rv5AZM+ywm1Hd1vtHR/z/Wus8Sg6ViIieRmFZXtXFh+jlF0i6QesGI735XYVbne7RBx5/v06cgvrbpppleemWuNyKo7uFknaVnodj0fkyFaRX562zkVF2jC98HgmNLDG4pu/iTwfIbLmI2tuHefMF0lcZy3vWSZyIEZk/sMib/cJfB3uXCQy607r3OxdKbJ0osiMMSVftwdjrWti3ZfWMbucVnnaPqv+nqXW+v41IjnHROb/yypf+opV/tlwb70icyVxvciKN635OKGBSFKcVZ5zzL/e8bL5DxfPCx5P8bIN061rXETku7Gl56GiEteLLH6+bHVLAcRICTm1XK/nldVpv553Kg5vhDM6Qq26Vduv3RzbY72/e85NFdvunl+t15FCva8WHloPTc6E4JDytePIsV6ZrBNesfEFcmyPdfcaXOvkdb+7y3r97v71lR9XZROxHsyHtbLGOrRB+dvIS4esw9CkW9nqO3Ig6wic0eHkdUtTkG29lnrJBAipBzsWwjejrU8ldUt4xbMyeDzgcZZ/fp+G0349r7yqJVErpdR/sdN+PU8ppVT10UStlFI2p4laKaVsThO1UkrZnCZqpZSyOU3USillc5qolVLK5jRRK6WUzWmiVkopm9NErZRSNqeJWimlbE4TtVJK2ZwmaqWUsjlN1EopZXOaqJVSyuY0USullM1polZKKZvTRK2UUjaniVoppWxOE7VSStncSRO1Maa1MWapMWabMWarMeaBqghMKaWUJbgMdVzAv0RknTGmPhBrjFkkItsqOTallFKU4Y5aRA6LyDrvchYQB7Ss7MCUUkpZyvUdtTEmEjgHWBNg21hjTIwxJiYlJaWCwlNKKVXmRG2MqQd8BzwoIpknbheRj0QkWkSiIyIiKjJGpZT6SytTojbG1MRK0l+JyPeVG5JSSqmiyvLWhwE+AeJEZFLlh6SUUqqostxRnw/cAlxsjNng/RlWyXEppZTyOunreSKyAjBVEItSSqkA9C8TlVLK5jRRK6WUzWmiVkopm9NErZRSNqeJWimlbE4TtVJK2ZwmaqWUsjlN1EopZXOaqJVSyuY0USullM1polZKKZvTRK2UUjaniVoppWxOE7VSStmcJmqllLI5TdRKKWVzmqiVUsrmNFErpZTNaaJWSimb00StlFI2p4laKaVsThO1UkrZ3EkTtTHmU2NMsjFmS1UEpJRSyl9Z7qinApdXchwAZOY72ZKY4Vfm8Qgi4ltftiOZPSnZVRHOSb0wfxtP/2D9/soucOF0e06rPbdHuOrdP/h1e1KZ6he43L6xyS5wkedwn1b/1SXf6cbtkWLl+47l4HCd3pj+tziQmkvM3tRi5W6PBBybkmQXuPyul7LweISthzJOXrEKFLjc7DuWc9J6Ho+wJv4Yiel5Abcv3Z5MclY+SZn5FR0ieQ43WxIz8JTjvJyu4JNVEJHfjTGRVRALd06NYe3eVPa8PIygGoYvVu/zJcLYpy7lSGY+t332JwB7Jw4HrAneqlFt1u1PZ1bsQV4Z1Z1Ve47RoHYwZzZvgDEGt0dYsfsot366loUPDKBVo9qE1gyiZpD1e+qRmRvZcCCdxQ8PxOMRClwegmoYPvp9D7uSsxl7YXucbmHku39YcV7Qjm7NG/DJigQAolo2YPx3mwFY//RgkrLyyXW46dmqISt2H2XroUzuHtQBAKfbQ6cnF/odw3HD31rO9iNZPDB9A3PuPZ9vYw7Sr104F3Vt4qsjIhhjSMrMp9/LSwD47u7zuOb9lQCsfeISEo7mcDTbwZCzmpKUmU/jeiEcycin7Rl1MMbw05bDiECnpvVYk5BKu8Z1qVsrmNCaQVw2+XceGdKZr9fsZ/G/BlKnljVFjmUXcCAtj3ohwTQLCyUkuAYJR3NIzXHQNzKcrHwXL/64jQs6NaZjk3q+sQeIT8km4WgOl3RrCkByVj4Gw4Mz1rMrKZvkrAJ6tWlIntPDxV0jWBOfysRrunPppN+5Ibo1vSMbsWhbEtdHt6Z/hzOoF1I4bfOdbkJrBgGwbn8aEfVCqB8aTHquk8jGdUudb79uTyLP4eHirk0IqmH4zy87SMtxcM9FHWkbXoesAhc/bjpMg9rBXNSlCb/tTKF/+zNoVLcWu5KyeH7+Nl4a2Z1ZsQcY0bMFS7en8NKCODZOGEJY7Zp4PMIXq/dxQ5/Wvhj3pGTTIqw2Lo+H+qE12Z2cRbOw2gx4dalvThxMy+W3nSks3Z7M+v3pHMtx+M2VzQcz+Hrtfi6Paka7M+rS5ow67DuWw8DXlvnqDDmzKXdc0I5ebRvx1ep9PDtvG1/e0Y8LOjUuNpfaP7EAgDn3nk9OgYu/TVnD6D6taXNGHV79aQcAvdo05Pt7zmfuxkM8Omsj+U4P0W0bMea8SIac2dR3fJsPZvDQtxvIznex+olLAMjKd5KZ7yLHezPRvGEoTeqHsuFAOvVCgkjPddK4XgjpeU7fNdavXThn1KvFgs1H+Oz2Pgzo2Bi3CJsOZrA2IZUaxvDvn7YDkPDKMNJznTSqWwuAnUlZ3D71T99xHh+7uMOZtAmvg8Pl4XBGPi0b1SbP4aaGgS9X76NJg1A6RNSjb7twDBCzL40+kY3428drEITpY/vj9gjdnvnJ1/b6pwfTqG4tnG4PMXvT6NcunBo1TKnz7lSYsvz29Sbq+SISVZZGo6OjJSYmptzBRD72Y7nqj+zZgh82HCp3P1WldXhtDqQG/o0PUDPI4HSffPyv692KmbEHqRVUA8dp3rWfio/HRHPX5+U/n4HUMFDRNyLntg9ndXzxu9GiRvZswaVnNuW+r9cD8Oq1Z/PorE2n1N+gLhEs25FSap3h3Zvz4+bDvvV/XtyRt3/dfUr9VaQ3R/fkq9X7WRvg7v1kLuwcwe87Sz/uovpGhp9SP3Z1y7lt+WL1vmLlPVqFsfFg4SeSE2/AysoYEysi0QG3VVSiNsaMBcYCtGnTpve+fcUP6GTKm6iVUspuKiNRV9hbHyLykYhEi0h0RERERTWrlFJ/efp6niqTiPoh1R2CUn9ZZXk97xtgFdDFGHPQGHNHZQXzwc29fctv3NDDt9y4Xi2m3t6nWP1re7fyLdcPtR4w3di3DQBdm9Xnt3GDuD66ld8+Z3gfOPRrFw5YDwJ3vjiUz27vQ9zzhS+3bHxmCJd2a8I/Bnbw279pgxAWPjCAK85uHvAYRvVqSXTbRrx/Uy82PzuEr+/sR//2Z3B2qzBCgmvQtVn9gPvddl5ksVhLMqx7MwAGn9mUey/qwMRR3YvV6dWmIbtfGkrNIP8HG5/d1ofNzw7xO1aAm/q1CdhX77aNmHp7H/588lLObhXmKz+3fbhvOSS4Bl/f2Y9/De4MwAOXdCLhlWF8fVc/7hnUoVibUHieKsJnt/XhzdE9fev/N7A9jerU5MNbetM8LNRX3q15A3q0CuOG6Nbsemmor/y28yJ9y6P7tCbhlWH8eP8FrHt6MA9d2pkuTQvP2bjLuvDj/Rcw/vKuAJgiw/v9Peex8IEBLH/0IhY+MMBX/vjQrr7lkuYNWOP4/FVnATCiR4sS6zUIDSYk2Lp0x/Rvy6e3RdM+wIPTr+7sx4WdrU+3NQxcdlZTRvdp7du+88WhrHr8Yq7r3Yr7L+nE40O7Mvue83huxFkl9n2iDhF1Gd498DGN7tOasRe2960vfvhCEl4Zxn+u68FDl3Yucx8A9wzqQHTbRr71Hq0b+pb7tgv3G2MoPn7tI6zx6dWmIfPuu6BY+5ef1axInANZ//RgHhnSmR/uPd831x8eXDzmr+/qx96Jw7n/4o6E1qzBoocuLNdxlZmIVPhP79695VTFp2TLhv1pIiLi8XhkSdwRcbs9IiKSmeeQd5fukjyHS15esE2y851yKD1Xcgqc5e7H5fbIc3O3yr6jOX7lRzLyJD3X4VcWuy9VHC53wDZe/3m7pOUU+OItK4/HI3GHM4qV/74zWVzuwnaOZReI2+2RmL2psv9YjuQ5XAHbS80ukH8vjJMvV+8tFn9JElKyZe/RbFm4+ZA4XW5xuz3i8XjkaFa+pGYXSFpOgeQ7i/d3/Dhj9h6T+JRsyc4v2/g7XW75bEW873wlZeTJ4m1HJO5whrQdP1/ajp8vIiK5BS5JzsyXA6k5cjAtN2Bbu5Iype34+XLRa0tP2u+nK+Kl7fj5kp5TtnEJZPa6g9J2/Hw5lF4Yz96j2aWe84WbD8sfu1KKlQca1/iUbN88L+rPhGNy17Q/JbfAJUkZeRKz91iJ/SVn5svnKxNk/f60Es/J3qPZMui1pXIkI6/EdkREEtNy/ebEDR+ulLbj54vT5ZYpy+Mlz+HyG4uibv10jbz203bfem6BS7YmFp/rItZcik/J9i3f+1Ws35gt3Z4k6735wOX2yKMzN8r+YzmSnuuQtuPnyx1T/yzW3vFzsmrPUZm2MkEOpOaIx+ORpMzSj7k8YvYeE2eAnHA6gBgpIaeW6WFieZ3qWx/qr2vh5sOE1KzBxV2blnmf469mGlP661DH5/jJ6qnSifd1PrvYkphBh4h61K4VVN2hVIjSHiae9D1qparC0BI+PpemdXidMtWzU3L5b2a3cYxqGXbySv8j9GGiUkrZnCZqpZSyOU3USillc5qolVLK5jRRK6WUzWmiVkopm9NErZRSNqeJWimlbK5S/jLRGJMClP/fObU0Bo5WYDgVReMqH42rfDSu8vlfjKutiAT8p0crJVGfDmNMTEl/RlmdNK7y0bjKR+Mqn79aXPrVh1JK2ZwmaqWUsjk7JuqPqjuAEmhc5aNxlY/GVT5/qbhs9x21Ukopf3a8o1ZKKVWEJmqllLI52yRqY8zlxpgdxpjdxpjHqqH/vcaYzcaYDcaYGG9ZuDFmkTFml/e/jbzlxhjzljfWTcaYXhUYx6fGmGRjzJYiZeWOwxhzq7f+LmPMrZUU17PGmETvmG0wxgwrsu1xb1w7jDGXFSmv0PNsjGltjFlqjNlmjNlqjHnAW16tY1ZKXNU6ZsaYUGPMWmPMRm9cz3nL2xlj1nj7mGGMqeUtD/Gu7/ZujzxZvBUc11RjTEKR8erpLa+yue9tM8gYs94YM9+7XrXjVdL/o6sqf4AgYA/QHqgFbATOrOIY9gKNTyh7FXjMu/wY8G/v8jBgIWCAc4E1FRjHhUAvYMupxgGEA/He/zbyLjeqhLieBR4JUPdM7zkMAdp5z21QZZxnoDnQy7tcH9jp7b9ax6yUuKp1zLzHXc+7XBNY4x2Hb4HR3vIPgLu9y/cAH3iXRwMzSou3EuKaClwboH6VzX1vuw8DXwPzvetVOl52uaPuC+wWkXgRcQDTgauqOSawYpjmXZ4GjCxS/rlYVgMNjTHl/39JBSAivwOppxnHZcAiEUkVkTRgEXA5p6GEuEpyFTBdRApEJAHYjXWOK/w8i8hhEVnnXc4C4oCWVPOYlRJXSapkzLzHne1dren9EeBiYJa3/MTxOj6Os4BLjDGmlHgrOq6SVNncN8a0AoYDU7zrhioeL7sk6pbAgSLrByl9UlcGAX4xxsQaY8Z6y5qKyGHv8hHg+P95tarjLW8cVRnffd6Pnp8e/3qhuuLyfsw8B+tuzDZjdkJcUM1j5v0YvwFIxkpke4B0EXEF6MPXv3d7BnBGVcQlIsfH6yXveL1hjAk5Ma4T+q+M8zgZeBTweNfPoIrHyy6J2g4uEJFewFDgXmPMhUU3ivX5pdrfZbRLHF7vAx2AnsBh4PXqCsQYUw/4DnhQRDKLbqvOMQsQV7WPmYi4RaQn0Arrrq5rVccQyIlxGWOigMex4uuD9XXG+KqMyRhzBZAsIrFV2e+J7JKoE4HWRdZbecuqjIgkev+bDMzGmsBJx7/S8P432Vu9quMtbxxVEp+IJHkvLg/wMYUf5ao0LmNMTaxk+JWIfO8trvYxCxSXXcbMG0s6sBToj/XVQXCAPnz9e7eHAceqKK7LvV8hiYgUAJ9R9eN1PjDCGLMX62uni4E3qerxOp0v2CvqBwjG+tK/HYUPTM6qwv7rAvWLLK/E+l7rNfwfSL3qXR6O/4OMtRUcTyT+D+3KFQfWnUcC1sOURt7l8EqIq3mR5YewvoMDOAv/ByfxWA/FKvw8e4/9c2DyCeXVOmalxFWtYwZEAA29y7WB5cAVwEz8H47d412+F/+HY9+WFm8lxNW8yHhOBiZWx9z3tj2IwoeJVTpeFZZcKmAQhmE9Gd8DPFnFfbf3DuJGYOvx/rG+W1oC7AIWHz/h3snxrjfWzUB0BcbyDdZHYifW91h3nEocwN+xHljsBm6vpLi+8Pa7CZiLfxJ60hvXDmBoZZ1n4AKsrzU2ARu8P8Oqe8xKiataxww4G1jv7X8L8EyRa2Ct99hnAiHe8lDv+m7v9vYni7eC4/rVO15bgC8pfDOkyuZ+kXYHUZioq3S89E/IlVLK5uzyHbVSSqkSaKJWSimb00StlFI2p4laKaVsThO1UkrZnCZqpZSyOU3USillc/8P+YEjBG0+9KUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(np.mean(np.sort(mu,axis=0),axis=1))\n",
    "print(p)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "plt.plot(np.sort(mu,axis=0).T)\n",
    "\n",
    "print(np.sort(mu,axis=0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "racial-external",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}