GS_sandbox.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "GS_sandbox.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/TomMaullin/BLMM-sandbox/blob/master/GS_sandbox.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TIMtrhFKB3Ay",
        "colab_type": "text"
      },
      "source": [
        "# GS implementation in python\n",
        "\n",
        "This code implements the Gibbs Sampling algorithm for estimating the parameters of linear mixed effects models as described in [my pdf](https://www.overleaf.com/read/wytrsgqgxzsk)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GNFQ0-MpQuBm",
        "colab_type": "text"
      },
      "source": [
        "## Pip Installations\n",
        "\n",
        "Pip install everything.**bold text**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UX02P0KvBWJr",
        "colab_type": "code",
        "outputId": "4d368562-a588-458f-dac3-206d4bf36d9f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 298
        }
      },
      "source": [
        "!pip install numpy\n",
        "!pip install pandas\n",
        "!pip install scipy\n",
        "!pip install matplotlib"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (1.16.4)\n",
            "Requirement already satisfied: pandas in /usr/local/lib/python3.6/dist-packages (0.24.2)\n",
            "Requirement already satisfied: pytz>=2011k in /usr/local/lib/python3.6/dist-packages (from pandas) (2018.9)\n",
            "Requirement already satisfied: numpy>=1.12.0 in /usr/local/lib/python3.6/dist-packages (from pandas) (1.16.4)\n",
            "Requirement already satisfied: python-dateutil>=2.5.0 in /usr/local/lib/python3.6/dist-packages (from pandas) (2.5.3)\n",
            "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.6/dist-packages (from python-dateutil>=2.5.0->pandas) (1.12.0)\n",
            "Requirement already satisfied: scipy in /usr/local/lib/python3.6/dist-packages (1.3.0)\n",
            "Requirement already satisfied: numpy>=1.13.3 in /usr/local/lib/python3.6/dist-packages (from scipy) (1.16.4)\n",
            "Requirement already satisfied: matplotlib in /usr/local/lib/python3.6/dist-packages (3.0.3)\n",
            "Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (2.5.3)\n",
            "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (2.4.2)\n",
            "Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (0.10.0)\n",
            "Requirement already satisfied: numpy>=1.10.0 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (1.16.4)\n",
            "Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (1.1.0)\n",
            "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.6/dist-packages (from python-dateutil>=2.1->matplotlib) (1.12.0)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.6/dist-packages (from kiwisolver>=1.0.1->matplotlib) (41.0.1)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q4JYRICVBtjl",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "## Python Imports\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YkTBWbRKQ5ah",
        "colab_type": "text"
      },
      "source": [
        "We need:\n",
        " - `numpy` for matrix handling.\n",
        " - `scipy` for sparse matrix functions.\n",
        " - `pandas` for quick reading and writing of csv files.\n",
        " - `os` and `sys` for basic commandline functions\n",
        " - `time` for timing functions.\n",
        " - `matplotlib` for making displays."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tebSlxvBBruv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import cvxopt\n",
        "import pandas as pd\n",
        "import os\n",
        "import time\n",
        "import scipy.sparse\n",
        "import scipy.sparse.linalg\n",
        "import scipy.stats\n",
        "from matplotlib.pyplot import imshow\n",
        "import sys"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ilAB3qmDMHa8",
        "colab_type": "text"
      },
      "source": [
        "## Helper Functions\n",
        "\n",
        "This section contains miscellaneous functions used to help the `GS` function including functions to work out the duplication matrix.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "jt3iahAIYzuL"
      },
      "source": [
        "## Toy Dataset\n",
        "\n",
        "This section read ins and formats a toy dataset. The files used here were generated in `R` and with **True** values (those with postfix `True`) being those used to generate the data and **Estimated** (those with postfix `REst`) values being the estimates `R`'s `lmer` package generated from this data."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qy39zwuhkn4C",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Make a data directory\n",
        "if not os.path.isdir('/Data'):\n",
        "  os.mkdir('/Data')\n",
        "  \n",
        "os.chdir('/Data')\n",
        "\n",
        "# Clone small git repo containg some csv files.\n",
        "if not os.path.isdir('/Data/BLMM-testdata'):\n",
        "  !git clone https://github.com/TomMaullin/BLMM-testdata.git\n",
        "  "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2EWsCZjCQcc9",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "### Z matrix\n",
        "\n",
        "The below reads in Z and makes an image of Z transpose.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aoKOwHqcEKDT",
        "colab_type": "code",
        "outputId": "e4d71593-58c6-403e-9f8b-310c4ad741aa",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 286
        }
      },
      "source": [
        "# Read in random effects design matrix and convert it into it's sparse format in\n",
        "# cvxopt.\n",
        "Z_3col=pd.read_csv('/Data/BLMM-testdata/Z_3col.csv',header=None).values\n",
        "Z = scipy.sparse.csr_matrix((Z_3col[:,2].tolist(), \\\n",
        "                            ((Z_3col[:,0]-1).astype(np.int64), \\\n",
        "                             (Z_3col[:,1]-1).astype(np.int64))))\n",
        "\n",
        "# Create an image of Z'\n",
        "imshow(Z.toarray().transpose(), \\\n",
        "       interpolation='nearest', vmin=-5, vmax=5, aspect='auto')\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.image.AxesImage at 0x7fe14a6bdda0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 5
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFnNJREFUeJzt3X2QXFWZx/Hfk5kkM3khLxrHJENM\nhoQgIhswQCIuOyVi2Ihi1VIqo27czW5E1zWIL4C6peyqBborYolAMO5EpAwIrLCRZWAx6KIYSEyD\ngQA9BJAkkPAOCQKZ5Nk/+nbP7Z5+f5mZzPl+qqa677n3nnPu7dPP3LndzxxzdwEAwjBqqDsAABg8\nBH0ACAhBHwACQtAHgIAQ9AEgIAR9AAgIQR8AAkLQB4CA1BT0zexUM3vIzHrN7Lx6dQoA0BhWbUau\nmTVJeljSKZK2S7pH0pnu/kChfZomjvfmaZOrai+jjz9OAITl9Se2P+Pu0+pRV3MN+x4vqdfdt0mS\nma2VdLqkgkG/edpktX/r0zU0KR3Y1VLT/gBwsHns7C88Xq+6arlsninpidjy9qgMADBMNfxeiZmt\nMLONZrbxwMt7G90cAKCIWoL+DkmHxpbbo7Is7r7K3Re6+8JRE8fX0BwAoFa13NO/R9I8M5ujVLD/\niKSuYju8feKzGv+T0Rpzyz3q2ZnQfa+/qqueW6z7jnX17ExoyYwFkqSenQlJ0pIZCzLP08u9Fy+q\nocsAELaqg76795nZZyT1SGqS9GN3v79uPQMA1F0tV/py95sl3VynvgAAGowvvQNAQKpOzqrG2FmH\n+ozPn11THaPaXq1pf77nD+Bg89jZX9jk7gvrURdX+gAQEII+AASEoA8AASHoA0BAavrK5lBou2Gs\nxl+3QZIyCV09OxN63wmnqe+J7ZnytCUzFmj3p9+pzV/9IcldAILHlT4ABISgDwABIegDQEAI+gAQ\nkIMuI7dWtWb0SmT1AhhcZOQCAKpC0AeAgBD0ASAgBH0ACEhwQb910zh1dCU0/q7U44eP3KSOroSS\nnd3q6ErouLc8rv1Pt2TKJt3emlmXLgOAg1VwQR8AQkbQB4CAEPQBICDBJWfVAwleAAYTyVkAgKoQ\n9AEgIAR9AAgIQR8AAnLQTZc4HMz73E7dvPlWLZmxQFJqesbP7DhBP5i5ITN9Y/pRkj65fbF+88tj\nNOuC32XKDrvmrCHrP4BwcaUPAAEh6ANAQAj6ABAQgj4ABISM3CFCVi+AcpGRCwCoCkEfAAJC0AeA\ngBD0ASAgZOQOkY6uxICydLaupEy2b/vvJ2j1rDslSae+5Xj5vtclSa+edry2nzwIHQUwonClDwAB\nIegDQEBKBn0z+7GZ7TazLbGyqWZ2m5klo8cpje0mAKAeSiZnmdlJkvZI+om7HxWVfVvSc+5+oZmd\nJ2mKu59bqjGSs+qr1gQvkruAg8OgJme5+28kPZdTfLqkNdHzNZI+WI/OAAAaq9p7+m3u/mT0/ClJ\nbXXqDwCggWr+INdT94cK3iMysxVmttHMNu7fs7fW5gAANag26O8ys+mSFD3uLrShu69y94XuvrBp\nwvgqmwMA1EO1yVk3SVom6cLo8ca69Qhlm/nT0brjR1dqn+/X0b/7hLaeeFXWdI3SwISv+LreixcN\nSb8BDJ1yvrL5M0l3SZpvZtvNbLlSwf4UM0tKek+0DAAY5kpe6bv7mQVW8U8AAOAgQ0YuAASEoA8A\nAWG6xIAxZSNwcGC6RABAVQj6ABAQgj4ABISgDwABYbrEgCU7uzPPl8xYIN3erp63rutfVn9G75IZ\nC7Ti4W1adXhHpuydnztLu0jqBQ4qXOkDQEAI+gAQEII+AASE5CzUhAQvoPFIzgIAVIWgDwABIegD\nQEAI+gAQEJKzUJOOrlSiVs/OhE7Z+n6NOvmJrCka04740ae0+mOX6l87jpUkXfmnOzWreYKWzFjA\ntI3AIOJKHwACQtAHgIAQ9AEgIAR9AAgIGbkYcmT1AsWRkQsAqApBHwACQtAHgIAQ9AEgIGTkYsgl\nO7v1rWfm69dHt6pnZ0JLZizIPKb1XnWM5n58s6RU9u/Sh5bq5vk3Z7YhqxcoD1f6ABAQgj4ABISg\nDwABITkLI0KtCV4kd2E4IzkLAFAVgj4ABISgDwABIegDQEBIzsKIcGh3s16b1KTfXnJFWdunE8DS\nz0nuQii40geAgBD0ASAgBH0ACEjJoG9mh5rZejN7wMzuN7OVUflUM7vNzJLR45TGdxcAUIuSGblm\nNl3SdHf/g5lNlLRJ0gclfULSc+5+oZmdJ2mKu59brC4ycjFcMWUjhrNBzch19yfd/Q/R85clbZU0\nU9LpktZEm61R6hcBAGAYq+ievpnNlnSMpA2S2tz9yWjVU5LaCuyzwsw2mtnG/Xv21tBVAECtyg76\nZjZB0vWSznb3l+LrPHWPKO99Indf5e4L3X1h04TxNXUWAFCbsoK+mY1WKuBf7e43RMW7ovv96fv+\nuxvTRQBAvZTMyDUzk7Ra0lZ3/25s1U2Slkm6MHq8sSE9BAZBR1cqOzedpXvqrIV69OojNfvD92Vl\n7v759OPVeuPdWdM6ph12zVmD33GgQuX8G4YTJX1c0h/NLD3Cv6xUsL/WzJZLelzShxrTRQBAvZQM\n+u5+pyQrsPrk+nYHANBIZOQCQECYLhGoExK80ChMlwgAqApBHwACQtAHgIAQ9AEgIEyXCNRJsrNb\nUv9UjPFHqT/x68m+PfrErHdl9iPBC4OJK30ACAhBHwACQtAHgIAQ9AEgIGTkAsMIWb3Ih4xcAEBV\nCPoAEBCCPgAEhKAPDCPJzm4lO7vV0ZXI/CQ7u9X3WlOmfN43XpHvaM2sO7DfNG/lDiU7uzX5ttah\nPgQMcwR9AAgIQR8AAkLQB4CAEPQBICAkZwEjTK0JXiR3DT8kZwEAqkLQB4CAEPQBICAEfQAICNMl\nAiPMzxdfoXPnnKDP9j6o7889Qj07E7p2zyR9aMKLkpSZvjFt1FFH6H9uXZsp77140aD3GYOHK30A\nCAhBHwACQtAHgIAQ9AEgIGTkAsjClI3DDxm5AICqEPQBICAEfQAICMlZALIkO7uzltNJWz07E3nL\n4+tI8Br+uNIHgIAQ9AEgIAR9AAhIyaBvZi1mdreZ3Wtm95vZBVH5HDPbYGa9ZnaNmY1pfHcBALUo\nmZxlZiZpvLvvMbPRku6UtFLSOZJucPe1Zna5pHvd/bJidZGcBYSBBK/6GtTkLE/ZEy2Ojn5c0rsl\nXReVr5H0wXp0CADQOGXd0zezJjNLSNot6TZJj0h6wd37ok22S5rZmC4CAOqlrKDv7vvdfYGkdknH\nSzqi3AbMbIWZbTSzjfv37K2ymwCAeqjo2zvu/oKk9ZIWS5psZunkrnZJOwrss8rdF7r7wqYJ42vq\nLACgNiUzcs1smqR97v6CmbVKOkXSRUoF/zMkrZW0TNKNjewogIPH5ENe0Qltjyt53GuSpK9sS+ik\nllTGbuuv2/SLeT1aMmOBFt27TxdMu3/AFI49OxM67JqzhqLrI145/4ZhuqQ1Ztak1F8G17r7OjN7\nQNJaM/uGpM2SVjewnwCAOigZ9N39PknH5CnfptT9fQDAQYKMXAAICEEfAALCdIkAhiWyevsxXSIA\noCoEfQAICEEfAALCdIkAhqWOrtQUjPFpGuNTNy49+mRZc7P2zXmz7K57Zce9XbfceJUO7/6U5nz5\nLklM25gPV/oAEBCCPgAEhKAPAAEh6ANAQEjOAjBi1ZrgNVySu0jOAgBUhaAPAAEh6ANAQAj6ABAQ\ngj6AEaujK6FkZ7c6uhI6/Iu7M8+Tnd2acOe4zPP9z47N2jb9OBIR9AEgIAR9AAgIQR8AAkLQB4CA\nkJELAAUMlykbycgFAFSFoA8AASHoA0BAmC4RAAro6EoUnK7x5I8vV/Ptm9SzM5EpT645VpPvGqtp\nl6ema3zPlpd1+W2nDH7Hi+BKHwACQtAHgIAQ9AEgIAR9AAgIyVkA0ED1SPDaduZXSc4CAFSOoA8A\nASHoA0BACPoAEBAycgGggTq6Elr+8KP60IQXtbTzb3TzHddL6s/ujUtn97bddYh2LX4pkw3cVMf+\ncKUPAAEh6ANAQMoO+mbWZGabzWxdtDzHzDaYWa+ZXWNmYxrXTQBAPVRypb9S0tbY8kWSLnb3uZKe\nl7S8nh0DANRfWRm5ZtYuaY2kb0o6R9L7JT0t6c3u3mdmiyV93d2XFKtnbMdMn7HyHI1qe1UHdrVk\nHiVlPU8rtD5329x1kjL1p58XaqNYW7mZdKWmPSunX/mU079y+pzbx3xl8f0KrS+3L+W8JvHjLrS+\nktek1P6lXoN8fcqtq1RbpcZtNSo5jnhZ+hgqqa+S163Y8ZY6H4XGZDnv1XLayqfccVftfpW+TvmW\nc+vO149cQzFd4vckfUnSgWj5DZJecPe+aHm7pJn16BAAoHFKBn0zO03SbnffVE0DZrbCzDaa2cYD\nL++tpgoAQJ2U8z39EyV9wMyWSmqRdIikSyRNNrPm6Gq/XdKOfDu7+ypJq6TU7Z269BoAUJWSV/ru\nfr67t7v7bEkfkfQrd/+opPWSzog2WybpxlJ1vX3is5Kkud9K3dc6/EtPZ9YlO7sHbP/zxVfkXZ+7\nbUdXImtdS+vrmfL4tu+e+1DBvsW3a+ptzdp/8ZxtWW0U0vdK9u/Qf3tH/ylJdnaroyuhWaubNGbL\nOCU7uzVpfWtWu60bx5VsIy7ep46uxIA+5iuL75fs7Na8Tz+at+63zXyyZPvt057P25d03ZJ02Hf2\n5e1Hen2x83rIxFf6n/+6Ne/+uWNh80mXZy1P/E3/fvHjjvepo+2ZrLpGPZrd1vRrx2Ttnzve6iG3\nntzllsTAsZHs7Nbh5z9bVn2F+tzx/QN5t+voSqjlD+MG7JvveaH1fS+M0YTx2Z9l5fZj/3NjM89z\nx0LfS2OytpWkE2Y/pmLSdbRsHqf3zd8y4DzEY0pcfDxNvaVVb535VN56873Hii3nex2KjZl6jadi\navme/rmSzjGzXqXu8a+uT5cAAI1S0b9hcPc7JN0RPd8m6fj6dwkA0Chk5AJAQAj6ABCQYTddYr6E\nqtzlQskV5WwX3za+fal9Cu2fr45y+1FO3em6Kkk6K7S+VFJOvP58SWXFEpgKtVPu+nKSr8pN0MrX\nbr5jq0Y1bebuF5d7TnP3KbRdsUStUolG8f5Uc2zFjrnUe65UIlmpRLlS7RRan+/9E1dOW+WO60Jj\nN378lb6HhyI5CwAwAhD0ASAgBH0ACAhBHwACMuyCfm7GWjqjM23uRa/l3e+QO7KzW4tleqbbKJU9\nV0jfy6PV9+dmdXQl9C/HrhtQd6EMxXKk+xV/jJfnayO33/P++U9Z66fdlPowaP27flBwn6wM1Zy6\n563cMSADNt+5WnH0/w3YP7eNQvvmOuw7+wbUtXpR9n7pelo3Fc9k7uhKqO36VObn/qdr+4+Yb7ly\n1IB+FWpz9g+zl9OPJ8x+TMnO7gHnctL61qx6t/zVlfrliZeqbepLkqTf/uWlWdvHxV+jN/3X2Exb\nyc7urKzk3P7Ezf326wXrL0excxI/3nzr0pofzM4CnvWm5ypqp9D6ePsdXQm1r2keEAdO6ugt2rdC\n7bb/ZHTWtnMvfHXA9p9fcJvm/mMyszzv6y9nrd///NisOufP2FVx7CjXsAv6AIDGIegDQEAI+gAQ\nEII+AARk0DNy2799VsEsxFLZb8UyRYsplpFZTvZqoW2LlVcypV4lU6hVM61eOcdV6VR65WQWlttW\nJVPclfMaVDo2CmVTxusolcGaK19/y80uLZXZWuy9UklWbKFM4GLHWOg9W+6YyLe+UJuljjmfUhm/\nlSj0nsk3voq9D8t9TYod57Yzv0pGLgCgcgR9AAgIQR8AAjLoQT83QSeegPTGdf33vT4w/74B+65d\nvCqrnrEt+8pKYIgnYOSblm7mVaML9q9QPbnHlDud38fednfWfqX6V+4UapVOz5ZbXui40kY90pq3\nvFAduX3/xTsvyzzPl0j3H8ddm3n+9Xf8d1Z99lirDuxuKXpMcz/5SOEpILen+n7IHa0a/cA4zbx6\n9IDtcuUeU/z1P+zi/UWndYz3Pz2OOw9LasqkvVn9Su8779/2Dqjji8fcmv9YSvS3YALg9oGJWPF9\n8yX1JTu71fxw/gS3eFJZvO10+cyfjh7Qh3F3jyuaoFjofZQ275t/ztuHeHLTsqN+r46uhG498Qfy\nnS1FpyGsNtmsddxrefdPt3X4V/qnC537Dw9nbeNumeezLzMVkuzs1qTbU6/Zm6/rnzoy3la9k7S4\n0geAgBD0ASAgBH0ACAhBHwACMqjJWWb2sqSHBq3B4e2Nkp4Z6k4ME5yLfpyLfpyLfvPdfWI9Kmqu\nRyUVeKheWWUHOzPbyLlI4Vz041z041z0M7ON9aqL2zsAEBCCPgAEZLCD/qrSmwSDc9GPc9GPc9GP\nc9GvbudiUD/IBQAMLW7vAEBABiXom9mpZvaQmfWa2XmD0eZQMrNDzWy9mT1gZveb2cqofKqZ3WZm\nyehxSlRuZvb96PzcZ2bHDu0R1J+ZNZnZZjNbFy3PMbMN0TFfY2ZjovKx0XJvtH72UPa73sxsspld\nZ2YPmtlWM1sc6rgws89F748tZvYzM2sJZVyY2Y/NbLeZbYmVVTwOzGxZtH3SzJaV03bDg76ZNUm6\nVNJfSzpS0plmdmSj2x1ifZI+7+5HSlok6Z+iYz5P0u3uPk/S7dGylDo386KfFZIuG1jlQW+lpK2x\n5YskXezucyU9L2l5VL5c0vNR+cXRdiPJJZJucfcjJP2FUuckuHFhZjMlfVbSQnc/SlKTpI8onHHR\nLenUnLKKxoGZTZX0NUknSDpe0tfSvyiKcveG/khaLKkntny+pPMb3e5w+pF0o6RTlEpMmx6VTVcq\nb0GSrpB0Zmz7zHYj4UdSezSI3y1pnSRTKummOXeMSOqRtDh63hxtZ0N9DHU6D5MkPZp7PCGOC0kz\nJT0haWr0Oq+TtCSkcSFptqQt1Y4DSWdKuiJWnrVdoZ/BuL2TfnHTtkdlQYj+DD1G0gZJbe7+ZLTq\nKUlt0fORfo6+J+lLkg5Ey2+Q9IK790XL8ePNnIto/YvR9iPBHElPS/rP6FbXj8xsvAIcF+6+Q9K/\nS/qTpCeVep03KcxxkVbpOKhqfPBBbgOZ2QRJ10s6291fiq/z1K/mEf/VKTM7TdJud9801H0ZBpol\nHSvpMnc/RtJe9f8JLymocTFF0ulK/SKcIWm8Bt7uCFYjx8FgBP0dkg6NLbdHZSOamY1WKuBf7e43\nRMW7zGx6tH66pN1R+Ug+RydK+oCZPSZprVK3eC6RNNnM0v8GJH68mXMRrZ8k6dnB7HADbZe03d03\nRMvXKfVLIMRx8R5Jj7r70+6+T9INSo2VEMdFWqXjoKrxMRhB/x5J86JP5cco9WHNTYPQ7pAxM5O0\nWtJWd/9ubNVNktKfsC9T6l5/uvxvo0/pF0l6MfZn3kHN3c9393Z3n63Ua/8rd/+opPWSzog2yz0X\n6XN0RrT9iLjydfenJD1hZvOjopMlPaAAx4VSt3UWmdm46P2SPhfBjYuYSsdBj6T3mtmU6C+n90Zl\nxQ3SBxZLJT0s6RFJXxnqD1AG4XjfpdSfZvdJSkQ/S5W6B3m7pKSk/5U0NdrelPqG0yOS/qjUNxqG\n/DgacF46Ja2LnndIultSr6SfSxoblbdEy73R+o6h7nedz8ECSRujsfELSVNCHReSLpD0oKQtkq6S\nNDaUcSHpZ0p9lrFPqb8Al1czDiT9fXROeiX9XTltk5ELAAHhg1wACAhBHwACQtAHgIAQ9AEgIAR9\nAAgIQR8AAkLQB4CAEPQBICD/D6dLkr+Yz5X4AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rkd2LdwKLWtP",
        "colab_type": "text"
      },
      "source": [
        "### Estimated Random Effects matrix\n",
        "\n",
        "The below reads in the Random effects variance predicted by `R`'s `lmer`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tOqQRAROqdkv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Read in estimated variance\n",
        "RFXVar_REst = pd.read_csv('/Data/BLMM-testdata/estd_rfxvar.csv',header=None).values"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RStOLwF_LlTE",
        "colab_type": "text"
      },
      "source": [
        "### Y vector\n",
        "\n",
        "The response vector is read in here:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DG8eWNdpPQOs",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "Y=pd.read_csv('/Data/BLMM-testdata/Y.csv',header=None).values"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RHm6DjdbPTvg",
        "colab_type": "text"
      },
      "source": [
        "### X matrix\n",
        "\n",
        "The fixed effects design matrix is read in here. It consists of an intercept and two random (Gaussian) columns.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xqTKH4n_Po8e",
        "colab_type": "code",
        "outputId": "b69b2631-3f21-47cd-ad40-62de037e5264",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 286
        }
      },
      "source": [
        "X=pd.read_csv('/Data/BLMM-testdata/X.csv',header=None).values\n",
        "\n",
        "# Image of the first 20 rows of X\n",
        "imshow(X[1:20,:])"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.image.AxesImage at 0x7fe147e4fe80>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 8
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFYAAAD8CAYAAADt0VN/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADBNJREFUeJztnXuMXVUVh78fbXnVhjeUIm9Kk0Kw\nkKaCqCkWEUZClaC2MVAERVQUjGhQIxiMRkOUaCCQSsvDFGjkZYXyqECCGEBK00JLqRQCgQIFWuiD\nAnXo8o+zBy+393b23HNXe+/p+pLJnHvOuufs+WbPvuexZm2ZGUH72WZLN6CqhFgnQqwTIdaJEOtE\niHUixDoRYp0IsU4M3tINaMQOu2xnO40Ymh2/z+B3smMXvr1Hdmzvirf4YO07yn5DDR0pdqcRQzn9\nxgnZ8b/Zc1527CF/Pzc79rVf/zE7tp5SQ4GkEyUtkbRU0kUNtm8naWba/pikA8ocr5toWaykQcCV\nwEnAaGCypNF1YWcDb5nZIcDlwO9aPV63UabHjgOWmtnzZrYeuBmYWBczEbg+Ld8CTJDU0pjVbZQR\nuw/wUs3rl9O6hjFm1gusAnYrccyuoWNOtySdI2mupLnr3np/SzenNGXELgP2rXn98bSuYYykwcBO\nwIpGOzOzqWY21szG7rjLdiWa1RmUEfs4MFLSgZK2BSYBs+piZgFT0vJpwAO2lTyyaPk81sx6JZ0H\n3AsMAqab2SJJlwJzzWwWMA34i6SlwEoK+VsFpS4QzGw2MLtu3cU1y+8BXylzjG6lI6+8BkrPPkdl\nx5427/Hs2BlD17XSHKCDzgqqRoh1IsQ6EWKdCLFOhFgnQqwTIdaJEOtEiHUixDpRiXsFM176V3bs\ncX/8cXbsqpWPtNIcIHqsGyHWiRDrRIh1IsQ6EWKdCLFOlMnd2lfSg5KelrRI0vkNYsZLWiVpfvq6\nuNG+qkiZC4Re4EdmNk/SMOAJSXPM7Om6uH+a2ckljtOVtNxjzexVM5uXltcAi9k4d2urpS2XtCnv\n9UjgsQabj5G0AHgFuNDMFrXjmLX0/OLC7Nj3R+Un4lgJO6XFSvoYcCtwgZmtrts8D9jfzNZK6gHu\nAEY22c85wDkAw/besWyztjhlM7qHUEidYWa31W83s9VmtjYtzwaGSNq90b4iKS6REoinAYvN7A9N\nYob3JRpLGpeO1zDbsGqUGQqOBU4HnpI0P637GbAfgJldTZFh+B1JvcC7wKTINuwHM3sY2GTau5ld\nAVzR6jG6mbjyciLEOhFinQixToRYJ0KsE5V4/P2Z7ze6RdGYOx4alx1rJbpd9FgnQqwTIdaJEOtE\niHUixDoRYp0IsU6EWCdCrBOVuKS955ajs2M3HPTf/B1v0/pTpOixToRYJ0qLlfSCpKdS0tvcBtsl\n6U+pWtyTkvKrNnQx7RpjjzOzN5tsO4ki+2Uk8EngqvS90myOoWAicIMVPArsLGnvzXDcLUo7xBpw\nn6QnUv5VPTkV5SpX0KwdQ8GnzWyZpD2BOZKeMbOHBroTM5sKTAUYftiuXZ8tU7rHmtmy9P114HaK\nYpK15FSUqxxlsw2HpmxuJA0FTgAW1oXNAs5IZwdHA6vM7NUyx+0Gyg4FewG3p4TCwcCNZnaPpHPh\nw8S42UAPsBRYB3yj5DG7grKV4p4HPtFg/dU1ywZ8r8xxupFK3CtYd2j+WcQZRz2aHTtt6NpWmgPE\nJa0bIdaJEOtEiHUixDoRYp0IsU6EWCdCrBMh1olKXNIOWb5tduwNjxybHbti7UaP8LKJHutEiHUi\nxDoRYp0IsU6EWCdCrBNlSpeMqilUNl/SakkX1MVEQbOBYmZLgDHw4YyfyyjyCuqJgmYlmAA8Z2Yv\ntml/XU+7xE4Cbmqy7RhJCyTdLemwNh2v42lHQbNtgVOAnzbYvFkKmvWOyH/8PeiN/PsKbGh9Ct12\n9NiTgHlmtrx+QxQ0K8dkmgwDUdCsRVIi3OeBb9esq83bioJmrWBm71A35XRd3lYUNAvaS4h1IsQ6\nEWKdCLFOhFgnKvH4+67P5p/RrdywfXbst659vZXmANFj3QixToRYJ0KsEyHWiRDrRIh1IsQ6EWKd\nCLFOhFgnKnGvoOf+H2THnnREfZ2K5qzond1Kc4DosW5kiZU0XdLrkhbWrNtV0hxJz6bvuzR575QU\n86ykKe1qeKeT22OvA06sW3cRcL+ZjQTuT68/gqRdgUsoCpiNAy5p9guoGlliU7mnlXWrJwLXp+Xr\ngS81eOsXgDlmttLM3gLmsPEvqJKUGWP3qqlG9BpF4Z16soqZVZG2fHil7JZSGS5VqxRXRuzyvhqF\n6Xuj5xjZxcwiKe7/zAL6PuWnAH9rEHMvcIKkXdKH1glpXeXJPd26CXgEGCXpZUlnA78FPi/pWeD4\n9BpJYyVdA2BmK4FfAY+nr0vTusqTdeVlZpObbJrQIHYu8M2a19OB6S21roupxCXtoLfzf4yHXjo4\nO3bt+tbH+rikdSLEOhFinQixToRYJ0KsEyHWiRDrRIh1IsQ6EWKdqMS9gtFjX8iOXfjk/tmxG9YP\naqE1BdFjnQixToRYJ0KsEyHWiRDrRL9im+RtXSbpmTS55O2Sdm7y3k1OVFllcnrsdWycFjQHONzM\njgD+Q+MKRn0cZ2ZjzGxsa03sTvoV2yhvy8zuM7Pe9PJRikSMoIZ2jLFnAXc32dbfRJWVpWwVo58D\nvcCMJiHZE1WWKWj21OL9smMPHZ0/XePKHdYPqB21lKnGeSZwMvD1ZiWfMiaqrI2N3C1JJwI/AU4x\ns3VNYnImqqwsOadbjfK2rgCGUfx5z5d0dYodIanvPyL2Ah6WtAD4N3CXmd3j8lN0IP2OsU3ytqY1\niX2FYkbPphNVbi3ElZcTIdaJEOtEiHUixDoRYp0IsU5U4vH3zgvzf4zn1uzbf1Di/XUDKJReR/RY\nJ0KsEyHWiRDrRIh1IsQ6EWKdCLFOhFgnQqwTlbikHfRuftWUH37xzuzYy65d1UpzgOixbrSaFPdL\nSctqZu/safLeEyUtkbRU0kZ1uapMq0lxAJenZLcxaXa5j5Bm/rySYoa60cBkSaPLNLabaCkpLpNx\nwFIze97M1gM3UxRB2yooM8ael/Jjpzcpq7fVFjOD1sVeBRxMMeHvq8DvyzYkCpoBZrbczD4wsw3A\nn2mc7JZdzCztM5Li+irEJb5M42S3x4GRkg5Mc9dOoiiCtlXQ7wVCSoobD+wu6WWKsqXjJY2hSCx+\ngTTLp6QRwDVm1mNmvZLOo6gMNwiYbmaLXH6KDsQtKS69ng20Xo+5i4krLycqca9gzYH5sUvWDc+O\nfW/DkBZaUxA91okQ60SIdSLEOhFinQixToRYJ0KsEyHWiRDrRCUuaZecdVV27EH/OCs7dtV7Df9R\nPYvosU6EWCdCrBMh1okQ60SIdSLEOpHzlHY6RVGd183s8LRuJjAqhewMvG1mYxq89wVgDfAB0Ls1\nFTXLuUC4jqIGzA19K8zsa33Lkn4PbCqR9Dgze7PVBnYrOY+/H5J0QKNtkgR8Ffhce5vV/ZQdYz8D\nLDezZ5tsz64UV7XcrbL3CiYDN21ie3alODObCkwFGH7YrgOaMfTYJ0/ND149gEfaH2ggzfgIZSrF\nDQZOBWY2ixlIpbiqUWYoOB54xsxebrQxKsX1Q5NKcVBkD95UFxuV4hKtJsVhZmc2WBeV4hJx5eVE\niHUixDoRYp0IsU6EWCfUpLz2FkXSG8CLdat3Bzb3XbJRZjaslTd2ZF6Bme1Rv07S3M19P7fMbCMx\nFDgRYp3oJrFTu+mYHfnhVQW6qcd2FR0ntr9yJ5K2kzQzbX+s2fO4ARxvX0kPSnpa0iJJ5zeIGS9p\nVU2plov73bGZdcwXxT8zPwccBGwLLABG18V8F7g6LU8CZpY85t7AUWl5GMX8ZPXHHA/cOZD9dlqP\nzSl3MhG4Pi3fAkxIT4tbwsxeNbN5aXkNsJg2VALpNLE55U4+jLFisrZVwG7tOHgaVo4EHmuw+RhJ\nCyTdLemw/vbVkVdeWwJJHwNuBS4ws9V1m+cB+5vZ2lQK6w5g5Kb212k9NqfcyYcx6UnxTsCKMgeV\nNIRC6gwzu61+u5mtNrO1aXk2METS7pvaZ6eJzSl3MguYkpZPAx6wEifjaXyeBiw2sz80iRneN45L\nGkfhbdO/zC19JtDgU7qH4pP5OeDnad2lFJOyAWwP/BVYSvH096CSx/s0RcbOk8D89NUDnAucm2LO\nAxZRnKU8Cnyqv/3GlZcTnTYUVIYQ60SIdSLEOhFinQixToRYJ0KsE/8D9NqNWg2rPVoAAAAASUVO\nRK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "W_sXGD8qzskp",
        "colab_type": "text"
      },
      "source": [
        "### Number of Levels and Parameters\n",
        "\n",
        "The number of levels is given by a vector with one entry for each grouping factor. e.g. nlevels=[10,2] means there are 10 levels for factor 1 and 2 levels for factor 2. \n",
        "\n",
        "The number of parameters is given by a vector with one entry for each grouping factor. e.g. nparams=[3,4] means there are 3 variables for factor 1 and 4 variables for factor 2."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_e5bUA2DztCf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "nlevels = np.array([20,3])\n",
        "nparams = np.array([2,2])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0byKXygpCa1P",
        "colab_type": "text"
      },
      "source": [
        "### True b values\n",
        "\n",
        "The true recorded values of the random effects b vector in this example."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LzwuuWSdCbCJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "b_True=pd.read_csv('/Data/BLMM-testdata/true_b.csv',header=None).values"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Bvfns6c-CbPd",
        "colab_type": "text"
      },
      "source": [
        "### True beta values\n",
        "\n",
        "The true fixed effects parameters used to generate this example."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mlvolWwvCbbP",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "beta_True=pd.read_csv('/Data/BLMM-testdata/true_beta.csv',header=None).values"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bVyBWvQK2Gw5",
        "colab_type": "text"
      },
      "source": [
        "### Product Matrices\n",
        "\n",
        "All products of matrices are calculated beforehand as it is both more computationally efficient and also similar to the setting we are interested in. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "e53HYCsj2G7I",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Z transpose Z\n",
        "ZtZ = Z.transpose()*Z # This works for products involving sparse\n",
        "\n",
        "# Z transpose X\n",
        "XtZ = X.transpose()*Z\n",
        "\n",
        "# XtX\n",
        "XtX = np.matmul(X.transpose(),X)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RhOHZ6F2DU5e",
        "colab_type": "text"
      },
      "source": [
        "## Perform Gibb's Sampling\n",
        "\n",
        "This section contains the function which performs Gibb's sampling for the linear mixed model\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PR3bUshSts74",
        "colab_type": "code",
        "outputId": "2d8895d9-282b-489f-88e4-592e001196c5",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        }
      },
      "source": [
        "def GS(XtX,XtY,XtZ,ZtY,YtY, beta0, u0, sigma_ffx0, sigma_rfx0, nit, n, nparams, nlevels):\n",
        "  \n",
        "  # Estimators of the population parameters\n",
        "  uhat = 0\n",
        "  betahat = 0\n",
        "  sigmahat_ffx = 0\n",
        "  \n",
        "  # Number of grouping factors\n",
        "  r = nlevels.shape[0]\n",
        "  \n",
        "  # We have a matrix estimate for each grouping factor\n",
        "  # variance matrix\n",
        "  sigmahat_rfx = dict()\n",
        "  for i in range(r):\n",
        "    sigmahat_rfx[i] = np.zeros((nparams[i],nparams[i]))\n",
        "  \n",
        "  # Current values\n",
        "  beta_curr = beta0\n",
        "  u_curr = u0\n",
        "  sigma_ffx_curr = sigma_ffx0\n",
        "  sigma_rfx_curr = sigma_rfx0  \n",
        "  \n",
        "  # Average values\n",
        "  beta_avg = np.zeros(beta0.shape)\n",
        "  u_avg = np.zeros(u0.shape)\n",
        "  sigma_ffx_avg = 0\n",
        "  sigma_rfx_avg = dict()\n",
        "  for i in range(r):\n",
        "    sigma_rfx_avg_i = np.zeros((nparams[i],nparams[i]))\n",
        "    sigma_rfx_avg[i] = sigma_rfx_avg_i\n",
        "  \n",
        "  \n",
        "  \n",
        "  # Constants needed for priors\n",
        "  b = 2\n",
        "  v = nparams+2\n",
        "  A = dict()\n",
        "  for i in range(r):\n",
        "    A_tmp = np.eye(nparams[i])\n",
        "    A_tmp[A_tmp==0]=0.5\n",
        "    A[i] = A_tmp\n",
        "  \n",
        "  # Run the number of iterations\n",
        "  for i in range(nit):\n",
        "    \n",
        "    # Current mean/variance of beta\n",
        "    betamean = np.matmul(np.linalg.inv(XtX),(XtY-np.matmul(XtZ,u_curr)))\n",
        "    betamean = betamean.reshape(betamean.shape[0])\n",
        "    betavar = sigma_ffx_curr*np.linalg.inv(XtX)\n",
        "    \n",
        "    # New estimate of beta\n",
        "    beta_curr = np.random.multivariate_normal(betamean, betavar)\n",
        "    beta_curr = beta_curr.reshape(beta_curr.shape[0],1)\n",
        "    \n",
        "    # Update average\n",
        "    beta_avg = beta_avg + beta_curr\n",
        "    \n",
        "    # D inverse\n",
        "    Dinv = np.zeros((np.dot(nparams,nlevels),np.dot(nparams,nlevels)))\n",
        "      \n",
        "    # Top right hand index of next block to be filled\n",
        "    TRHind=0\n",
        "    \n",
        "    for j in range(r):\n",
        "      \n",
        "      # Current sigma_rfx block\n",
        "      sigma_rfx_curr_j = sigma_rfx_curr[j]\n",
        "      \n",
        "      # Work out the next block\n",
        "      DinvBlock = np.linalg.inv(sigma_rfx_curr_j)\n",
        "      \n",
        "      # This gets repeated for the number of levels for factor j\n",
        "      for k in range(nlevels[j]):\n",
        "        \n",
        "        # Indices for current block\n",
        "        blockIndices = np.ix_(TRHind+np.arange(nparams[j]),TRHind+np.arange(nparams[j]))\n",
        "        \n",
        "        Dinv[blockIndices] = DinvBlock\n",
        "        \n",
        "        # Move Top Right hand index to next block\n",
        "        TRHind = TRHind + nparams[j]\n",
        "    \n",
        "    # Current mean/variance for u\n",
        "    umean = np.matmul(np.linalg.inv(ZtZ.toarray()+sigma_ffx_curr*Dinv),(ZtY-np.matmul(XtZ.transpose(),beta_curr).reshape(ZtY.shape)))\n",
        "    umean = umean.reshape(umean.shape[0])\n",
        "    uvar = np.linalg.inv(ZtZ.toarray()+sigma_ffx_curr*Dinv)*sigma_ffx_curr\n",
        "    \n",
        "    # New estimate of u\n",
        "    u_curr = np.random.multivariate_normal(umean, uvar)\n",
        "    u_curr = u_curr.reshape(u_curr.shape[0], 1)\n",
        "    \n",
        "    # Update average\n",
        "    u_avg = u_avg + u_curr\n",
        "    \n",
        "    # Residual sum of squares (Y-Xbeta-Zu)'(Y-Xbeta-Zu)\n",
        "    rss_curr = YtY + np.matmul(np.matmul(u_curr.transpose(), ZtZ.toarray()),u_curr) + \\\n",
        "          np.matmul(np.matmul(beta_curr.transpose(), XtX),beta_curr) - \\\n",
        "          2*np.matmul(beta_curr.transpose(), XtY) - \\\n",
        "          2*np.matmul(u_curr.transpose(), ZtY) + \\\n",
        "          2*np.matmul(np.matmul(beta_curr.transpose(), XtZ),u_curr)\n",
        "    \n",
        "    # New estimate of fixed effects vairance\n",
        "    sigma_ffx_curr = scipy.stats.invgamma.rvs(a=(b+(n/2)), scale=rss_curr[0]/2, size=1)\n",
        "    \n",
        "    # Update average\n",
        "    sigma_ffx_avg = sigma_ffx_avg + sigma_ffx_curr\n",
        "    \n",
        "    # Index of the first value in the next uj\n",
        "    uj_ind = 0\n",
        "    for j in range(r):\n",
        "      \n",
        "      # Current sigma_rfx block\n",
        "      sigma_rfx_curr_j = sigma_rfx_curr[j]\n",
        "      \n",
        "      # Work out the indices of uj we need\n",
        "      uj_indices = uj_ind + np.arange(nparams[j]*nlevels[j]) \n",
        "      \n",
        "      # Get u for this grouing factor\n",
        "      uj_curr = u_curr[uj_indices]\n",
        "      \n",
        "      # Initialis sum of uut \n",
        "      suut = np.zeros(nparams[j],nparams[j])\n",
        "      \n",
        "      # Now we need sum ujl*ujl' over each level l\n",
        "      for l in range(nlevels[j]):\n",
        "      \n",
        "        # Work out the uj observations for level k\n",
        "        ujl_curr = uj_curr[np.arange(nparams[j])+nparams[j]*l]\n",
        "        suut = suut + np.matmul(ujl_curr,ujl_curr.transpose())\n",
        "      \n",
        "      # Work out parameters to sample new sigma_rfx block from inverse wishart\n",
        "      sigma_rfx_curr_j_scalematrix = A[j] + suut\n",
        "      sigma_rfx_curr_j_df = v[j] + nlevels[j]\n",
        "      \n",
        "      # Sample new variance block from inverse wishart\n",
        "      sigma_rfx_curr[j] = scipy.stats.invwishart.rvs(df=sigma_rfx_curr_j_df, scale=sigma_rfx_curr_j_scalematrix)\n",
        "      \n",
        "      # Work out average\n",
        "      sigma_rfx_avg[j] = sigma_rfx_avg[j] + sigma_rfx_curr[j]\n",
        "      \n",
        "      # Update index\n",
        "      uj_ind = uj_ind + nparams[j]*nlevels[j]\n",
        "      \n",
        "  \n",
        "  # Convert sums to averages\n",
        "  beta_avg = beta_avg/nit\n",
        "  sigma_ffx_avg = sigma_ffx_avg/nit\n",
        "  u_avg = u_avg/nit\n",
        "  for j in range(r):\n",
        "    sigma_rfx_avg[j] = sigma_rfx_avg[j]/nit\n",
        "    \n",
        "  return(beta_avg, sigma_ffx_avg, u_avg, sigma_rfx_avg)\n",
        "  \n",
        "\n",
        "# Example\n",
        "\n",
        "# Product matrices\n",
        "XtX = np.matmul(X.transpose(),X)\n",
        "XtY = np.matmul(X.transpose(),Y)\n",
        "XtZ = np.matmul(X.transpose(),Z.toarray())\n",
        "ZtY = np.matmul(Z.toarray().transpose(),Y)\n",
        "YtY = np.matmul(Y.transpose(),Y)\n",
        "\n",
        "# Number of subjects,n\n",
        "n = 1000\n",
        "\n",
        "# Initial values\n",
        "\n",
        "# RFX variance\n",
        "r = nlevels.shape[0]\n",
        "sigma0_rfx = dict()\n",
        "for i in range(r):\n",
        "  sigma0_rfx[i] = np.eye(nparams[i],nparams[i])\n",
        "\n",
        "# FFX variance\n",
        "sigma0_ffx = 1\n",
        "\n",
        "# Beta\n",
        "beta0 = np.zeros((X.shape[1],1))\n",
        "\n",
        "# u\n",
        "u0 = np.zeros((Z.shape[1],1))\n",
        "\n",
        "t1=time.time()\n",
        "estimates=GS(XtX,XtY,XtZ,ZtY,YtY,beta0, u0, sigma0_ffx, sigma0_rfx, 1000, n, nparams, nlevels)\n",
        "t2=time.time()\n",
        "print(t2-t1)\n",
        "\n",
        "print(\"Results after 1000 iterations\")\n",
        "print(\"-----------------------------\")\n",
        "print(\"Beta estimates (Gibbs)\")\n",
        "print(estimates[0])\n",
        "print(\"Beta estimates (R)\")\n",
        "print(pd.read_csv('/Data/BLMM-testdata/estd_beta.csv',header=None).values)\n",
        "print(\"Beta true\")\n",
        "print(beta_True)\n",
        "print(\"-----------------------------\")\n",
        "print(\"U estimates (Gibbs)\")\n",
        "print(estimates[2].reshape(23,2))\n",
        "print(\"U estimates (R)\")\n",
        "print(pd.read_csv('/Data/BLMM-testdata/estd_b.csv',header=None).values.reshape(23,2))\n",
        "print(\"U true\")\n",
        "print(pd.read_csv('/Data/BLMM-testdata/true_b.csv',header=None).values.reshape(23,2))\n",
        "print(\"-----------------------------\")\n",
        "print(\"Random effects variance estimates (Gibbs)\")\n",
        "print(\"Factor 1\")\n",
        "print(estimates[3][0])\n",
        "print(\"Factor 2\")\n",
        "print(estimates[3][1])\n",
        "print(\"Random effects variance estimates (R)\")\n",
        "estd_rfxvar=pd.read_csv('/Data/BLMM-testdata/estd_rfxvar.csv',header=None).values\n",
        "print(\"Factor 1\")\n",
        "print(estd_rfxvar[np.ix_(np.arange(2),np.arange(2))])\n",
        "print(\"Factor 2\")\n",
        "print(estd_rfxvar[np.ix_(np.arange(2)+2,np.arange(2)+2)])\n",
        "print(\"Random effects variance true\")\n",
        "true_rfx_3col=pd.read_csv('/Data/BLMM-testdata/true_rfxvar_3col.csv',header=None).values\n",
        "true_rfx = scipy.sparse.csr_matrix((true_rfx_3col[:,2].tolist(), \\\n",
        "                            ((true_rfx_3col[:,0]-1).astype(np.int64), \\\n",
        "                             (true_rfx_3col[:,1]-1).astype(np.int64))))\n",
        "print(\"Factor 1\")\n",
        "print(true_rfx[np.ix_(np.arange(2),np.arange(2))].toarray())\n",
        "print(\"Factor 2\")\n",
        "print(true_rfx[np.ix_(np.arange(2)+40,np.arange(2)+40)].toarray())"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "4.513150453567505\n",
            "Results after 1000 iterations\n",
            "-----------------------------\n",
            "Beta estimates (Gibbs)\n",
            "[[0.5464067 ]\n",
            " [1.98512161]\n",
            " [3.03201756]]\n",
            "Beta estimates (R)\n",
            "[[1.04778749]\n",
            " [1.98321039]\n",
            " [3.03156038]]\n",
            "Beta true\n",
            "[[1]\n",
            " [2]\n",
            " [3]]\n",
            "-----------------------------\n",
            "U estimates (Gibbs)\n",
            "[[-0.16641273 -4.2168499 ]\n",
            " [ 1.4136973   1.41212541]\n",
            " [ 1.59112192 -1.58888116]\n",
            " [ 1.5582383   0.40464328]\n",
            " [ 1.3003867   0.93790526]\n",
            " [-1.20813631  2.19750368]\n",
            " [-0.62719399 -1.18432693]\n",
            " [-0.90428168 -1.34791901]\n",
            " [ 0.20105858  2.53792982]\n",
            " [-0.58205751 -3.43663566]\n",
            " [ 1.33329767  0.10709062]\n",
            " [-1.21503839 -2.43700162]\n",
            " [ 0.63534816  1.33610731]\n",
            " [-0.99885399 -0.08067159]\n",
            " [ 1.44451179  2.19572301]\n",
            " [-1.26277561  4.30747634]\n",
            " [ 0.2658458   0.87146507]\n",
            " [-0.7912232  -3.16064885]\n",
            " [-1.48857717 -0.24107696]\n",
            " [-0.47916287  3.13871743]\n",
            " [ 2.31250477  0.64337653]\n",
            " [ 2.57262529 -0.15742232]\n",
            " [-6.8652865   0.4375333 ]]\n",
            "U estimates (R)\n",
            "[[-0.16266309 -4.21724427]\n",
            " [ 1.41643004  1.41243783]\n",
            " [ 1.59527575 -1.58921871]\n",
            " [ 1.57314783  0.40497904]\n",
            " [ 1.29727291  0.93796984]\n",
            " [-1.20537013  2.19715236]\n",
            " [-0.62177368 -1.1845146 ]\n",
            " [-0.8982196  -1.34761407]\n",
            " [ 0.2105669   2.5381855 ]\n",
            " [-0.57561587 -3.43679857]\n",
            " [ 1.33808317  0.10724134]\n",
            " [-1.20676845 -2.43712662]\n",
            " [ 0.64262693  1.33629341]\n",
            " [-0.9916963  -0.08043681]\n",
            " [ 1.44782818  2.19618068]\n",
            " [-1.25815704  4.3074814 ]\n",
            " [ 0.27424557  0.87149959]\n",
            " [-0.78704685 -3.16062475]\n",
            " [-1.48916827 -0.24116661]\n",
            " [-0.47444214  3.13847411]\n",
            " [ 1.80633411  0.64368738]\n",
            " [ 2.06756827 -0.15742136]\n",
            " [-7.37272861  0.43752657]]\n",
            "U true\n",
            "[[-0.19301143 -4.22666242]\n",
            " [ 1.43735476  1.40883489]\n",
            " [ 1.34711832 -1.58788242]\n",
            " [ 1.81825209  0.3910735 ]\n",
            " [ 1.25654173  0.94063103]\n",
            " [-1.22775617  2.19293973]\n",
            " [-0.79408731 -1.19281546]\n",
            " [-0.9898057  -1.35221872]\n",
            " [-0.03319243  2.53178969]\n",
            " [-0.65604937 -3.43139032]\n",
            " [ 1.42910277  0.11018046]\n",
            " [-0.96244844 -2.43071144]\n",
            " [ 0.70088294  1.33487455]\n",
            " [-1.3137399  -0.07784482]\n",
            " [ 1.45666816  2.19521647]\n",
            " [-1.3663195   4.30629966]\n",
            " [ 0.22757187  0.87208539]\n",
            " [-0.89571743 -3.14622337]\n",
            " [-1.51292639 -0.23242295]\n",
            " [-0.5588196   3.15054584]\n",
            " [ 1.85205881  0.64080503]\n",
            " [ 2.10133112 -0.1570577 ]\n",
            " [-7.27154738  0.44280999]]\n",
            "-----------------------------\n",
            "Random effects variance estimates (Gibbs)\n",
            "Factor 1\n",
            "[[1.23276805 0.37282448]\n",
            " [0.37282448 4.85900331]]\n",
            "Factor 2\n",
            "[[15.4790992  -0.29463425]\n",
            " [-0.29463425  0.40732196]]\n",
            "Random effects variance estimates (R)\n",
            "Factor 1\n",
            "[[1.24022707 0.36034323]\n",
            " [0.36034323 5.0719078 ]]\n",
            "Factor 2\n",
            "[[20.68631165 -0.79595202]\n",
            " [-0.79595202  0.21015464]]\n",
            "Random effects variance true\n",
            "Factor 1\n",
            "[[2.  0.5]\n",
            " [0.5 4. ]]\n",
            "Factor 2\n",
            "[[6.  0.1]\n",
            " [0.1 1. ]]\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}