example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mrrce import MrRCE\n",
    "from simulations.simulation_utils import generate_data, model_error\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.linear_model import LinearRegression, RidgeCV, MultiTaskLassoCV\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Navigate\n",
    "\n",
    "1. [Generate random dataset](#Generate-random-dataset)\n",
    "1. [Fit with CV](#Fit-with-CV)\n",
    "1. [Estimations](#Estimations)\n",
    "1. [Compare to other estimators](#Compare-to-other-estimators)\n",
    "1. [Fit with no CV](#Fit-with-no-CV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Generate random dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# settings\n",
    "params = dict(\n",
    "n = 50,                          # num obs\n",
    "p = 20,                          # num predictors\n",
    "q = 5,                           # num tasks\n",
    "sigma = 1,                       # coeff variance\n",
    "corr_x = .7,                     # grid of values for rho (correlation coefficient)\n",
    "sigma_err = 1,                   # correlation coefficient for predictors\n",
    "err_corr = .9,                   # correlation coefficient for errors\n",
    "g_sparse_level = .0,             # group sparsity level\n",
    "sparse_level = .0,               # sparsity level\n",
    "err_cov_type = 'ar'             # error covariance form. One of ['ar', 'equi', 'fgn', 'identity']\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "rho = .4 # correlation coefficient for coefficients\n",
    "X, Y, B, Sigma, Sigma_X = generate_data(rho = rho, **params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Fit with CV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mrrce = MrRCE(verbose=True) # init MrRCE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2019-09-11 08:37:44,104: iter 1, loss 2.398151\n",
      "2019-09-11 08:37:44,147: iter 2, loss 0.615972\n",
      "2019-09-11 08:37:44,187: iter 3, loss -0.729327\n",
      "2019-09-11 08:37:44,227: iter 4, loss -1.576275\n",
      "2019-09-11 08:37:44,272: iter 5, loss -2.032753\n",
      "2019-09-11 08:37:44,314: iter 6, loss -2.249618\n",
      "2019-09-11 08:37:44,356: iter 7, loss -2.344970\n",
      "2019-09-11 08:37:44,401: iter 8, loss -2.385307\n",
      "2019-09-11 08:37:44,446: iter 9, loss -2.402189\n",
      "2019-09-11 08:37:44,490: iter 10, loss -2.409225\n",
      "2019-09-11 08:37:44,532: iter 11, loss -2.412099\n",
      "2019-09-11 08:37:44,576: iter 12, loss -2.413288\n",
      "2019-09-11 08:37:44,619: iter 13, loss -2.413782\n"
     ]
    }
   ],
   "source": [
    "mrrce.fit(X, Y) # fit MrRCE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "selected regularization parameter is 0.0015\n"
     ]
    }
   ],
   "source": [
    "print(f\"selected regularization parameter is {mrrce.alpha_best:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Estimations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.11, -1.43,  1.22,  1.23, -0.35],\n",
       "       [ 1.02, -0.1 , -0.17,  0.84,  1.56],\n",
       "       [ 0.32, -0.79, -2.04,  0.69, -0.52],\n",
       "       [-1.18,  1.34,  0.73,  0.72, -0.75],\n",
       "       [-0.08,  0.18, -0.13,  0.75,  0.57],\n",
       "       [-1.47, -1.71,  0.03, -1.74, -0.26],\n",
       "       [-0.72, -0.75, -0.58, -0.36,  1.02],\n",
       "       [-0.28,  1.49, -0.34, -1.32, -0.24],\n",
       "       [ 1.65,  1.62,  2.62,  0.96,  1.01],\n",
       "       [ 2.22,  0.67,  1.74,  0.02,  0.93],\n",
       "       [ 0.34,  1.5 , -0.28,  1.79,  0.37],\n",
       "       [-0.89, -1.72, -0.84, -2.02, -0.41],\n",
       "       [ 0.03,  0.32,  0.46,  0.11, -0.3 ],\n",
       "       [-0.35,  0.39, -0.11,  0.37, -1.24],\n",
       "       [ 3.12,  0.81,  0.63,  0.28,  2.55],\n",
       "       [-0.75,  1.23,  0.82, -0.06,  1.91],\n",
       "       [-0.41, -0.41, -1.34, -0.58,  0.48],\n",
       "       [ 0.48,  0.04,  1.29,  0.26, -0.63],\n",
       "       [-0.86, -1.14, -0.74, -1.39, -0.87],\n",
       "       [ 1.06,  0.07, -0.1 , -0.04,  0.3 ]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# coefficient matrix\n",
    "B.round(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.18, -1.52,  1.24,  1.27, -0.15],\n",
       "       [ 0.79, -0.17,  0.59,  0.89,  1.59],\n",
       "       [ 0.4 , -0.73, -1.82,  0.59, -0.43],\n",
       "       [-1.23,  1.24,  0.32,  0.84, -0.57],\n",
       "       [ 0.13,  0.29, -0.05,  0.67,  0.47],\n",
       "       [-1.42, -1.61, -0.62, -1.53, -0.31],\n",
       "       [-0.81, -0.7 , -0.4 , -0.4 ,  0.8 ],\n",
       "       [-0.46,  1.35,  0.17, -1.3 , -0.  ],\n",
       "       [ 1.74,  1.59,  2.33,  0.86,  0.89],\n",
       "       [ 2.37,  0.93,  1.56,  0.15,  1.01],\n",
       "       [ 0.17,  1.29,  0.48,  1.64,  0.44],\n",
       "       [-0.85, -1.69, -1.22, -2.01, -0.59],\n",
       "       [-0.01,  0.36,  0.84,  0.21,  0.15],\n",
       "       [-0.34,  0.48, -0.47,  0.37, -1.39],\n",
       "       [ 3.16,  0.68,  0.82,  0.14,  2.41],\n",
       "       [-0.71,  1.27,  0.59, -0.06,  1.89],\n",
       "       [-0.6 , -0.46, -0.88, -0.45,  0.64],\n",
       "       [ 0.48,  0.03,  1.26,  0.34, -0.54],\n",
       "       [-0.82, -1.  , -0.89, -1.35, -0.71],\n",
       "       [ 1.01, -0.01, -0.13, -0.14,  0.19]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mrrce.Gamma.round(2) # coefficient matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True (sigma, rho) = (1, 0.4)\n",
      "Estimated (sigma, rho) = (1.056, 0.426)\n"
     ]
    }
   ],
   "source": [
    "print(f\"True (sigma, rho) = ({params['sigma']}, {rho})\\nEstimated (sigma, rho) = ({mrrce.sigma:.3f}, {mrrce.rho:.3f})\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.28,  0.02, -0.84,  0.  , -0.26],\n",
       "       [ 0.02,  0.13, -0.01,  0.09,  0.05],\n",
       "       [-0.84, -0.01,  3.82, -0.11,  0.93],\n",
       "       [ 0.  ,  0.09, -0.11,  0.22,  0.13],\n",
       "       [-0.26,  0.05,  0.93,  0.13,  0.55]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# transform Sigma\n",
    "Sigma.round(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.18, -0.01, -0.47, -0.  , -0.17],\n",
       "       [-0.01,  0.13,  0.09,  0.1 ,  0.1 ],\n",
       "       [-0.47,  0.09,  3.56, -0.38,  0.71],\n",
       "       [-0.  ,  0.1 , -0.38,  0.29,  0.14],\n",
       "       [-0.17,  0.1 ,  0.71,  0.14,  0.49]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# estimated transformed Sigma\n",
    "mrrce.Sigma.round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Compare to other estimators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OLS\n",
    "lm = LinearRegression(fit_intercept=False).fit(X, Y)\n",
    "B_ols = np.matrix(lm.coef_.transpose())\n",
    "# Ridge\n",
    "ridge = RidgeCV(fit_intercept=False).fit(X, Y)\n",
    "B_ridge = np.matrix(ridge.coef_.transpose())\n",
    "# Group Lasso\n",
    "gl = MultiTaskLassoCV(fit_intercept=False).fit(X, Y)\n",
    "B_gl = np.matrix(gl.coef_.transpose())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Model Error:\n",
      "============\n",
      "MrRCE: 1.882\n",
      "OLS: 3.299\n",
      "GL: 2.645\n",
      "Ridge: 2.792\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Model Error\n",
    "print(\n",
    "f\"\"\"\n",
    "Model Error:\n",
    "============\n",
    "MrRCE: {model_error(B, mrrce.Gamma, Sigma_X):.3f}\n",
    "OLS: {model_error(B, B_ols, Sigma_X):.3f}\n",
    "GL: {model_error(B, B_gl, Sigma_X):.3f}\n",
    "Ridge: {model_error(B, B_ridge, Sigma_X):.3f}\n",
    "\"\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Fit with no CV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "mrrce = MrRCE(\n",
    "    alpha=3e-3,\n",
    "    use_cv=False, \n",
    "    verbose=True\n",
    ") # init MrRCE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2019-09-11 08:39:05,169: iter 1, loss 2.417206\n",
      "2019-09-11 08:39:05,217: iter 2, loss 0.663573\n",
      "2019-09-11 08:39:05,261: iter 3, loss -0.640970\n",
      "2019-09-11 08:39:05,304: iter 4, loss -1.443197\n",
      "2019-09-11 08:39:05,345: iter 5, loss -1.863802\n",
      "2019-09-11 08:39:05,384: iter 6, loss -2.059526\n",
      "2019-09-11 08:39:05,423: iter 7, loss -2.144564\n",
      "2019-09-11 08:39:05,461: iter 8, loss -2.180230\n",
      "2019-09-11 08:39:05,503: iter 9, loss -2.194865\n",
      "2019-09-11 08:39:05,542: iter 10, loss -2.201070\n",
      "2019-09-11 08:39:05,584: iter 11, loss -2.203620\n",
      "2019-09-11 08:39:05,625: iter 12, loss -2.204570\n",
      "2019-09-11 08:39:05,666: iter 13, loss -2.204956\n"
     ]
    }
   ],
   "source": [
    "mrrce.fit(X, Y) # fit MrRCE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.18, -1.52,  1.24,  1.27, -0.15],\n",
       "       [ 0.79, -0.17,  0.59,  0.89,  1.6 ],\n",
       "       [ 0.4 , -0.73, -1.82,  0.59, -0.43],\n",
       "       [-1.22,  1.23,  0.32,  0.84, -0.57],\n",
       "       [ 0.13,  0.29, -0.05,  0.67,  0.47],\n",
       "       [-1.42, -1.61, -0.62, -1.53, -0.31],\n",
       "       [-0.81, -0.7 , -0.39, -0.4 ,  0.8 ],\n",
       "       [-0.45,  1.35,  0.17, -1.3 , -0.  ],\n",
       "       [ 1.74,  1.59,  2.33,  0.86,  0.89],\n",
       "       [ 2.37,  0.93,  1.56,  0.15,  1.  ],\n",
       "       [ 0.17,  1.29,  0.48,  1.64,  0.45],\n",
       "       [-0.85, -1.69, -1.22, -2.01, -0.59],\n",
       "       [-0.01,  0.36,  0.84,  0.21,  0.15],\n",
       "       [-0.34,  0.48, -0.47,  0.37, -1.39],\n",
       "       [ 3.16,  0.68,  0.82,  0.14,  2.41],\n",
       "       [-0.71,  1.27,  0.59, -0.06,  1.89],\n",
       "       [-0.6 , -0.46, -0.88, -0.45,  0.64],\n",
       "       [ 0.48,  0.03,  1.26,  0.34, -0.54],\n",
       "       [-0.81, -1.  , -0.89, -1.35, -0.71],\n",
       "       [ 1.01, -0.01, -0.12, -0.14,  0.19]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mrrce.Gamma.round(2) # coefficient matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Model Error:\n",
      "============\n",
      "MrRCE (no GLASSO cv): 1.883\n",
      "OLS: 3.299\n",
      "GL: 2.645\n",
      "Ridge: 2.792\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Model Error\n",
    "print(\n",
    "f\"\"\"\n",
    "Model Error:\n",
    "============\n",
    "MrRCE (no GLASSO cv): {model_error(B, mrrce.Gamma, Sigma_X):.3f}\n",
    "OLS: {model_error(B, B_ols, Sigma_X):.3f}\n",
    "GL: {model_error(B, B_gl, Sigma_X):.3f}\n",
    "Ridge: {model_error(B, B_ridge, Sigma_X):.3f}\n",
    "\"\"\"\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}