statistics-resampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "My Octave Notebook \n",
    "============\n",
    "\n",
    "Interact with Octave and the statistics-resampling package online in Jupyter Notebook with Binder. All commands are interpreted by Octave.  Since this is a [MetaKernel](https://github.com/Calysto/metakernel), a standard set of magics are available.  Help on commands is available using the `%help` magic or using `?` with a command. For the best experience, run this notebook in [Jupyter Lab](https://mybinder.org/v2/gh/acpennlab/statistics-resampling-online/master?labpath=statistics-resampling.ipynb). \n",
    "\n",
    "The code cells and their output below illustrate some examples of loading data and then performing analysis using functions from the **statistics-resampling** package. Follow the links in the **Quick start** section of the [README page](https://gnu-octave.github.io/statistics-resampling/readme.html), or on the [Function reference page](https://gnu-octave.github.io/statistics-resampling/function_reference), to find further examples. The example data files provided, **data1.tsv** and **data2.csv**, can be edited directly in this instance of JupyterLab. \n",
    "\n",
    "The statistics-resampling package already comes installed and loaded at **statistics-resampling-online**. Visit [GNU Octave Packages](https://gnu-octave.github.io/packages/) to see what other packages are available and how to install them. To load a newly installed package called `<package-name>`, use the command `pkg load <package-name>`. To preview demonstration number `X` of a function called `<function-name>`, use the command `demo ('<function-name>', X)`. To get more information about a particular function called `<function-name>`, use the command `help <function-name>`. Note that the **parallel** package cannot be installed at **statistics-resampling-online**, and that additions to the binder environment may be required to build or run some other Octave packages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some Octave packages come already installed and in this statistics-resampling-online Binder environment. If you were running Octave on your own machine, you would need to install and load the statistics-resampling package first by running the code in the following two cells:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this command once to install the package\n",
    "pkg install -forge statistics-resampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this command in each new Octave session\n",
    "pkg load statistics-resampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the code in the next code cell if you wish to know:\n",
    "1) The version of GNU Octave that you are running\n",
    "2) What packages (and their version) are already installed (and loaded*) here\n",
    "3) More information about the statistics-resampling package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ver('Octave');\n",
    "pkg list;\n",
    "pkg describe -verbose statistics-resampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example 1** (below): Bootstrap confidence intervals for the correlation coefficient using `bootci` on (bivariate) data (in wide format) loaded from **data1.tsv**. Application of `randtest2` for comparing the mean of paired samples is also illustrated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ci =\n",
      "\n",
      "   0.5189\n",
      "   0.8739\n",
      "\n",
      "pval = 2.8182e-04\n"
     ]
    }
   ],
   "source": [
    "# Load data from tab-separated-value (TSV) file (i.e. tab-delimited text file) with header\n",
    "# Formatting: each %f corresponds to a number (double precision)\n",
    "[before, after] = textread ('data1.tsv', '%f %f', 'Delimiter', '', 'HeaderLines', 1);\n",
    "\n",
    "# 95% BCa bootstrap confidence intervals for the correlation coefficient\n",
    "ci = bootci (1999, @cor, before, after)\n",
    "\n",
    "# Randomization test for comparing means of paired samples\n",
    "pval = randtest2 (before, after, true, 10000, @(A, B) mean (A) - mean (B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example 2** (below): Bootstrap two-way ANOVA (Type I) using `bootlm` with data (in long format) loaded from **data2.csv**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ANOVA SUMMARY with effect of gender averaged over levels of degree\n",
      "F(1,18) = 87.20, p = 0.0001 for the model: salary ~ 1 + degree\n",
      "F(1,18) = 10.97, p = 0.00495 for the model: salary ~ 1 + degree + gender\n",
      "F(1,18) = 0.42, p = 0.567 for the model: salary ~ 1 + degree + gender + degree:gender\n",
      "\n",
      "MODEL FORMULA (based on Wilkinson's notation):\n",
      "\n",
      "salary ~ 1 + degree + gender + degree:gender\n",
      "\n",
      "MODEL COEFFICIENTS\n",
      "\n",
      "name                                   coeff       CI_lower    CI_upper    p-val\n",
      "--------------------------------------------------------------------------------\n",
      "(Intercept)                            +25.00      +24.11      +25.92           \n",
      "degree_1                               -8.000      -9.803      -6.046           \n",
      "gender_1                               +2.000      +0.2268     +3.812           \n",
      "degree:gender_1                        +1.000      -1.649      +3.549           \n",
      "\n",
      "\n",
      "MODEL FORMULA (based on Wilkinson's notation):\n",
      "\n",
      "salary ~ 1 + degree + gender + degree:gender\n",
      "\n",
      "MODEL ESTIMATED MARGINAL MEANS\n",
      "\n",
      "name                                   mean        CI_lower    CI_upper        N\n",
      "--------------------------------------------------------------------------------\n",
      "1, f                                   +25.00      +24.03      +26.07          8\n",
      "0, f                                   +17.00      +15.20      +19.25          4\n",
      "1, m                                   +27.00      +25.12      +29.00          3\n",
      "0, m                                   +20.00      +19.02      +21.03          7\n",
      "\n",
      "\n",
      "MODEL FORMULA (based on Wilkinson's notation):\n",
      "\n",
      "salary ~ 1 + degree + gender + degree:gender\n",
      "\n",
      "MODEL POSTHOC COMPARISONS\n",
      "\n",
      "name                                   mean        CI_lower    CI_upper    p-adj\n",
      "--------------------------------------------------------------------------------\n",
      "f - m                                  -2.500      -4.421      -0.5792      .012\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAAAXNSR0IB2cksfwAAAAlwSFlzAAALEwAACxMBAJqcGAAAACB0RVh0U29mdHdhcmUAR1BMIEdob3N0c2NyaXB0IDEwLjAyLjHDX9SJAAAgAElEQVR4nO3df2wb530/8EdtjMTKTErakAQVT3WU2pJBBugsw1qotfMQBfBlwzJ4k8hhMcp/xGwrUNaBGBoDZHQiikk+LSsdIINIFKNhDxPlwasFLOdBWmIgEF11kbHAx8ZyUjn1ycM8ZJKOmRUPSaHvH8/X1yspUbR0d89zd+/XH8b5dHzuc+QdP3zunh8N6+vrBAAAgLUvsQ4AAACAECQkAADgBBISAABwAQkJAAC4gIQEAABcQEICAAAuICEBAAAXkJAAAIALSEgAAMAFJCQAAOACEhIAAHABCQkAALiAhAQAAFxAQgIAAC4gIQEAABeQkAAAgAtISAAAwAUkJAAA4AISEgAAcAEJCQAAuICEBAAAXEBCAgAALiAhAQAAF5CQAMCp5ufn5+fnjQvgaA3r6+usYwAAz5mZmXnttdfu3r0bCASOHTuWSqW2UUhvby8tSl+o/1VdXV2EkNHR0dpbXrt2jRDi8/kikQjNeca9PNR+YUuoIXGqpaWlu7tbX95eCfryrVu3AoFAS0vLgQMHzp8/X2MlfSFFA6je+/biATCanp7et2/f2bNnjx07Rgjp7e3du3dvZ2fnN7/5zUgk0t3d/corr1y6dKm3tzcQCLz00ks/+MEPurq6Ojs7U6nU/Px8b29vZ2fnRx99VF0y/St91TvvvNPb29ve3v6Nb3yDvvxb3/rWT37yk2Kx+P7771dsnEgkjhw5cuTIEWN2PHjw4PLy8pe+9KULFy5UlP87v/M7c3NzxWIxm83a8p65HxIS1/Q8cf78+UAgEAgEbt261d3d3d3dfeXKle7ublEUux8ghBw/fjwQCOiZTPfiiy9+5zvfWV5e/ru/+7vvfOc7NVZSy8vLy8vL9LdhRQCvv/766urq8ePHLTxs8ICTJ08+99xzb7755tzcXHNzMyHkq1/96j/+4z8eOHBgbm5ufX09Ho+fP3/+3r17hw4d0jTtX//1X3//93//sccem5mZGRkZefzxx//6r/96w99G9K/nzp1ramo6ceJEa2vr3/zN33zlK1/57d/+7ccee0xRlMOHD4fD4S+++GJ+ft648aVLl+gdI+Pdv2vXrvl8vrW1tb6+vory9+7d++u//uvhcDgej9vynrnfI6wDgE1NTEx84xvf6OnpIYScPHny3XffnZ2djUajH3744ZkzZ44cOfLhhx/Oz893dXXRfwkhn3zyybvvvhsMBiuK+vTTT1977TVCyJEjR2qvpFpaWu7fv3/w4EF9jR7AG2+80dTUdO7cOSsOGbxjZGTk/v37P/zhD1Op1Ntvv/3ZZ5898sgvv458Pl9XV1dzc/Pa2tqzzz5LCLl58+Z//Md/lMvl5uZmmsCOHDmSyWQ++uijvXv3lstl/XSlf/36179OCPm///u/Rx99tK2t7fPPP//pT39KX24Mw7jx559/vmvXroo4Dx48qN+RozfoKl5iwXvjXagh8evpp58eGRmJRqOEkD179jz99NM9PT3lcpkQ8vLLL+vb6P8SQpaXl994443qovbs2UMrW8ePH29sbOzu7m5paalYadx+eXl5bW1tbm7OWIIxAIAdOnTo0DvvvPPMM89cunTpF7/4xe7du3/+85//2Z/92f3795988km6zTPPPKOqajabffvttz/99NOFhYW1tbWVlZV79+698847oVBoZWUlEol8/PHHBw8evHbtGr3//Mknn/z4xz9ub2//yU9+8sILL/z4xz9+8cUX33333Y8//lh/+b179/77v/+b7kLfOBAI1BO58SUHDhy4d++efjcPdmoduNTc3EwXDh8+3NzcfPTo0dbW1t27d587d07/E10w/tva2nr48OHdu3efOHFC32x9ff2dd95pbW1tbW3t7Ow8fPhwjZXGXRv/awygtbX1xIkT1h07eNDzzz///PPPs46i0nvvvffee++xjsJD0MrOc27duqXXqDZcef78eb0GBmAP+tiG3nkGz0JCAgAALuAZEgAAcAEJCQAAuICEBAAAXEBCAgAALiAhAQAAF5CQAACAC0hIAADABVclJEVRWIcAYB9FUSRJkmWZdSAA5nBPQpJleWxsjHUUAPaRZTkajRaLRdaBAJjDJaN9V9eNOjo66MLU1JTt4YAD6GeIcyWTSU3TfD6fcWVDQ0ONl1SPzILtPbU950PzuCQhFQoFn8+nqqpx5cLCghX7WlhYsO67zLrCUXJFyVYUa7NisaiqajKZNK7cv3//Qx3dw35DYXt3b8+WSxJSOp0mhJRKJdaBANinUChomibLcj6fZx0LgAlckpAoXJbgKZlMhnUIAGZyT6MGAABwNCQkAADgAhISAABwAQkJAAC4gIQEAABcQEICAAAuICEBAAAXkJAAAIALSEgAAMAFJCQAAOCCq4YOAqgHHf/YWYNOAngBEhJ4iJ6KeBvtW9O0ycnJ27dv0/+KohgOh9mGBGA/JCTwBM5rRZOTk8FgcGBggP5XlmVFUUKhENuoAGyGhATu19DQwG0qovRURImiuO2ibt68udkcbvRN0P9a8V9s74XtOb8QkJDAzTivGFlhywn6Kt6NLd8cbO/u7bmChATu5MFUBOB0SEjgNk5MRZIkVazp6elBuwbwGiQkcBX+HxdtKJlMsg4BgD0kJHAJJ1aMAMAIIzWA4zU0NNCKkQuykaZprEMAYAYJCRzMTalIkiRFUU6dOpVIJFjHAsAGEhI4lWtSEVUulwkhkUgElSTwLCQkcB69YsQ6EDOFQqFCoSAIQk9PD+tYANhAowZwEhe3XIhEIpFIhFSN2gDgHUhI4AwuTkWEkFgsZvxvPp9nEwcAU0hIwDt3pyLKmIEKhQK7QABYckxC0jTN7/ezjgLsUzFGJCFkcXGxvb29etkdJEm6evUqIcTn89F7dwBe44BGDfU0h1UURZIkWZbtDAysUz2k8czMTF9fX/Wya5TL5eHh4cHBQUEQNttGVdWK/xaLRUVRrI8OwA4OqCHpzWGz2exm28iyHI1Gs9mscdz+jo4OQsjU1JTpIVk6vZt1hTul5M7OTrpw48aNzs5OWrggCI8++mj1smv4fL7av6iKxeLY2NjFixf1NRMTExs2yaNnvsveH/ACByQk2hw2Ho/XaA6bTCY1TfP5fMaVFl2QCwsL9IJ3VuGOKNn4uKihoaGzs9P46KixsVHfkXF5Gzj8st5yOLtwOFxxhhNCZFlua2urmMqPw6MDqIcDElI9zWGLxaKqqhih0tEquhZxONG4pfSGdn6/P5PJ1POSaDQqCEIsFkNLcXAH3hNSnc1hC4WCpmmyLKO9rBPV2Y5uZmZmw2V30E/doaGhLTeWJCmZTE5MTJTL5RrPnGrzWrMR4B/vCanO5rB1/qIE3nihSXed9N9ewWBws230y4HeDNjJLQG9PqqnpZmZmVQqNT8/X7EMYBveExJBc1iXQirSKYoyNjZm807X19cr7pH29vY2NzdXLwPYxgEJiTaHLZfLaNXtGu4biW4nQqFQPp+XJEkURb/fX6M1aT1u3rxZUfvZbFlf0Lf/t3/7tw2Xa5eDZQctc37dOaAfEm0OOzs7yzoQMIErx0U1BX0aJAgC7eewbfv379ffXuNo6NXLRnT9888/v+Fy7XKw7KBlwjcH1JDQds4dqvu6kqon5/Pz811dXTYGxRFRFOmNu6NHj7KKwd3NRoB/DkhI22gOC1ypaM214VN0QsiFCxfGx8c9+z0YDofD4TDrKABYckBCeqjmsMAb43OIGk/RPd6gS5IkQkipVKL/Re8F8CYHJKR6msOCbSpax1X3Zam4NWe8f13jAdL4+Hhzc/Pi4qJFYXNOvy+NQYTBy7hOSEyaw0INxurO+oN2wwsLC/rocxUqck+NZ6q0aZmX60m0lV0ulyPoVwdexXUrO9ocNhgMDg4OptPpHbY+AtPRnGRcY/xvPU16Kp4YefYBEjEMIqxpGutYANjgOiFRZjWHhZ3b8LYbHQJVX69jGKcT0UGEBUGoMYgwgLtxfcuO4qE5LOiq78ItLCyga9HOHT16tFwuq6qK8xw8ywEJCc1heYYRgMxCn6JpmpbNZtHKDryJ94SE5rDc0lORpyaJsE65XA6Hw36/H63swLN4T0hoDssh1IqsIIoiHc8et+ygftUX45UrV44cOUIImZ+fn5yc/OpXv/oXf/EXrMJ7WA5o1CBJkqIop06dSiQSrGMBordcYB2I2wiCkMlkMplMOBxWVZV1OOAA+sWot249ffr0iy++SJdHRkYEQfirv/ordgE+NAckJDSH5QTGRbVOIpEYGhoaGhpKJBK2/fCiH6j+3/n5ecyB5HSvvfbaY489RpdPnjw5PT29e/dutiE9FAckJDSHtY1xVgLjV5WxSTej0FxO07R8Pq+qaiaTyefz254Etn7VP67Hx8f7+/tHRkas3jWYpfboJ9/97nc7Ozvv3r1rf2DbxvszJILmsHYxDr5AfrXTa3UHWHCllZUV1JCcpcZvxHA4vLKy4vP57IxnhxxQQ8pms6urq5qmYXBV+6FKZBs6ZmMsFovFYoqiWL276h/X/f39hJC+vj6rd12h4QGb92snOw9weXmZLoyOjmazWdSQTIbmsLapuE234TJYgUl/hopfG9PT04uLi4cOHWISiYvPsYoRIFmHwzUHJCQ0h7XHZpNE4BLyiNHRUdYhgNc5ICHR5rCEEE3TVFW14XmvN6F3ESsVo9oHg0FPzZJsfHLpSnr9z8XHaBbeE1IikdA0jY6sqmlaOp1mHZEL4Wphi45qT6ef8Pv9dAwhT3H9uef6AzQL7wmJNoeNxWIYNMgi7v5x6iB0VHu/349R7cGzeE9IYB1UjLiCUe0BHJCQ9OawhJDBwcFQKMQ2HhdAKuJQMBgslUrBYBBPScGzeE9IuFNnLqQibmH6CQCuE1L9rY8URZFlWRRF1J8qVHfv2LDPB1IUc9vrb4d2p+AmXCek+lsf5XK5wcHBsbEx2kCc6ujoIIRMTU2ZHpilMwCZXviNGzc6OzvpAiGkoaGBLmy4cnuse0O8M9nS9evXS6WSIAiLi4t1vqRYLI6NjV28eLFiPT3zvfPWgWtwnZCoelofBQIBQRAqhgO36IJcWFigF7xTCqeJhzz4nqJr9CpRR0cHrTBte7/WvSGWlmxFsTvR1NQkiiLt3lDnS8Lh8IYjlXF4dAD1cEBCqqf1Ec1VuHdhZLwvR+d11R8gGYdjwFMlTgSDwUgkUiwWPdgPCYByQEKqp/VRW1tbIpEIBAJ2BsazDXsX6WswMhCH6PPRcDj87//+76xjAWDDAQmpntZHAwMDtsbEMdR4nEiSJEJIqVSi/62/lR3a44GbOGD6iXK5HAqFMNr3ljCNnnMlk8m2trZ8Pp/P50VRZB0OABsOqCFhtO8toVbkdLFYrFQqybJMCBEEIRKJsI4IgAEHJKRsNiuKYiQSicVi+PFYDYPRuUA+n1cUBb3owOMckJCoXC7HOgTuoGLkJrIs693A8WQIvMkBCQnNYSlj+qmxjMbcDtXU1JTJZPCgFLzMAQnJs81ht8w0Gy5v+Crg3+rqaiKRoMv21JD0nmo4SYATvCekbTeHdTpjBjKmHHCreDw+OTlp22jfOK+AQ7w3+0Zz2Aq0VbexhfeGwy4YV4IjZLPZ1dVVTdOGhoZs2N2GY+wy0fAA60Cs4voDNBHvNSTPNoc1fl9Up5YNh1rA+AuOtr3RvneI+a8Wj1TUcAu9TrwnJDSHxRnsEfb3t+Ph1OKnomY1jxzmDvGekIgnm8OijZynVMz7pSiKp+5Ou/48r2hqBDU4ICF5rTks6vVeU/+8X+BQuKLr5ICEZH9zWFbwM8rL6pn3C8DdHJCQbG4Oy4SeijC1mmfVM+8XgLs5ICHVM/2Ec6FWBFQ4HA6Hw6yjAGDJAQmJSXNYG1SnIvSc97JYLEYX/H5/JpNhGgsAGw5ISI6efqIi68zPz09OTvb39x86dIg8aH5j7KOwsLDQ2dnJMGBgRa/929MxFoBDXCckdzSHNWad8fHxXC53+vRpgmoQ/Cq9hhQMBpkGAsAM1wmJz+awtAZjTCeLi4vt7e3kQQXohRde6O3tJVUNuGltaXp6Oh6PLy4uVlSe0FnB4+ikX4QQej8AwIN4H8uOPGgOKwgCD81hGxoabty4YexePjMz09fXR5dHRkZeeeWV8fFx+l99TDnyIDlNT08TQuj21XONr6+v37hxw7ZjAX7EYrGxsbFYLBaLxRRFYR0OABtc15AozpvD9vb2Njc30+ULFy6kUilaW6KM9+sIIdPT04uLi/QBEoBOHyLL4wNlgcc5ICFx1Ry29o21bDY7OjpK79eRjdrRjY6O2hImOIwkSU1NTVevXr19+3ahUEin06wjAmDAAbfsYg/o4zWwRW/Zbfinn/3sZ/F4vLm52Tg9hM3hgRMtLS0NDAz4fL50Oq2qKutwANhwQA2J/+awMzMzdGF0dBQNE2AbNE0rFovlclnTNJf1twOonwMSkoOaw2JcVNieSCQyOzsbj8ez2Syfz0oBbOCAhOSI5rCoGMFOiKJI+9glk0nWsQAww3tCqnPGWEVRZFkOhUL295xFKgJ70JNcFEW9GZ6qqqqq+nw+NMwDd+C9UUM+n//7v//7fD4/ODhYo+mRLMvRaLRYLNoZG1ougJ1yuVw0Gs3lcvqaiYkJhvEAmI73GlKdzWGTyaSmaT6fz7iyo6ODEDI1NWV6VHTEOdqJ1fQJI6ybgQIlc0uSpIo1PT09Fb0dAoGAIAiaphlXyrLc1tZWUUOiZ75H3jpwE94T0tLSUjKZpKlIb91gFIvFRFEUBEFV1Yr77xZdkJbeo1tYWKDfJijZ6pKtKHbb6nl0REcqMc4KFo1GBUGIxWIDAwPGLXk7OoA68X7LbsvmsPl8PhKJFAoFWZY3zFgmovfoavRDArBOW1tbIpGgtwFojWpiYmJoaMjFE1eC1/BeQ6qzOazV88cYa0X4+QmmKxQKgiDQe3SqqhaLxXA4XJFpjNUgWqNCkzxwGd4TEvPmsGhEBzagtfyhoSFVVYPBIL0XxzooALvxnpDYQkdXsE0kEtmsVwOARyAhbQwVIwAAm/HeqMF+6F0ETOgDkfA8IgmApVBD+iXUioCVWCx2/fr1LUckAXA3JKT/D4+LgKF8Pp9Kpf7yL/8SQ32DlyEhoWIEXPiN3/gNfcYvfcoVAE/x9DMkPC4CfpRKJX05FospisIwGAAmPFpDQq0IeCMIQjwe1zTt6tWrFUMBAXiEF2tIqBUBb+gsEoIgCIJw+/Zt1uEAsOGtGhIqRsAnQRCampokSSqXy/bP6QXACa8kJKQi4By9TbfZIMIAXuD+W3ZouQCOIEmSoiinTp3S29oBeI3LExJSETgFne4oEolUTMEH4B2uTUg3b95EX1dwkFAoRCeh6OnpYR0LABuufYa0f/9+TFwEDqKP9o023+BZrq0hAQCAsyAhAQAAF5CQAACAC0hIAADABSQkAADgAhISAABwAQkJAAC4gIQEAABcQEICAAAuICEBAAAXkJAAAIALrkpIiqLYsJc/+IM/cGLhKNmeku2kKAqdtGJ7L6czs9S/fUdHx/Z2VE8kFpVsXcwPVfJDvdWclMyEexKSLMtjY2OsowCwTy6Xi0ajuVxu2yWsr69blwzq1NDQsH//fuZhWEefBMf0Y7SuZFZcMtr3hj8SOflxxE/hKNllAoGAIAjV8yfV86bdvHlz//79dMv632Qnnvk8nJz6lnW+xLqSOef4hBSLxURRVBTF5/Opqqqvx9wT4Hp0Tj9BEIwr6z/z9Z/VbC+WhoaGmzdvMg/DOvoBWjE9G/0QXTPxm+MTUj6fJ4TQiWRKpRLjaABs1NbWlkgkAoHA9l7OybcYJ2FYx9IDdNm7hzlVAQCAC+5p1AAAAI6GhAQAAFz48ve+9z3WMVhCUZSzZ882NDRUPPI1q+Ryubxv3z5zSyaEqKrq9/vNLZMG7Pf7n3jiCXNLJtYETCx+k607NxjSPwhVVbPZrIkft16yuR+KsTTTT9Hqwq0IW19jStjGks39EK17N0zn2hoS7aJRKBRML1mW5Wg0WiwWTS+5WCwmEgnTi915b5XNWBQwsfJNJlaeG6wYP4ixsTFRFK9evWp6yea+b8aP2PRT1Fi4uedSRWkmdn80lmzuh1j9VnN78rs2IWUymYmJiW03QKohmUz6fD6fz2d6yeFw2IpiN+utsnMWBUysfJOJlecGK8YPIhAIFAoFs47OWHIoFMpms2ZVK40fsemnqLFwc88lY2nmDg1T8YaY+CEaSzb3QzSdCxNSLBYrFAqFQiGZTJrbEJyWXCwWL1++nEwmTS/ZxAKNNuytwjkr3mSdFecGP8rl8uDgoBUV4tXV1Xg8Xi6XjR3+ts34EZt+ihoLN/dcMpZWKBQURTHl3SBVb4iJH6KxZHM/RNM5vh9SNdozSZKkoaEhcx9v0JITiYSmabIs0/+aWLJFdthbhYlCoWD6m6y7ffu26ecGP+jH/dxzz1lReDab9fl8pmQO40ds+ilqLNzcc8lYWjqdJuZ1f6x+Q8z6ECveARM/RNOhHxIAAHDBhbfsAADAiZCQAACAC0hIAADABSQkAADgAhISAABwAQkJANzG3C6rYBskJPZisdixY8eOHTsWi8Xof7dXCCFEkqTaG2yPsfCdlAOwDbEHJEkynoQ1TsUth/PBacwnF3aMdZx8Pk+HmQqHw3RNLpdbXV1NJpOKotDe2sPDw6qq6svZbJYQcvDgwampKULIn/zJn5RKJX2sB/oqn8+XTqfp0BIVXQ5pD3NCCO2tnc/nY7EY7UJIN04mk5IkNTU1ra6uiqJoLNy4CxoM7QdOX2XtOwVepfdplSRJUZRSqfS3f/u39LQMBoMV14jxbJckKZlM0utLVVX99KZ/pac9/RenNA9QQ+JRf38/7f6dy+WOHj0aCAQmJydzudzw8PDAwEA2my2VSqIoTk1N0b9ev349GAxGIhH6qkKhMDw8LIqiqqrlcnl4eLhilEZFUdLpNJ1m18i4calUomGEQiFj4ZQxsMuXLw8MDIiiaPn7Al5Fa0iFQkE/IU+cOEFPy+prpKenx/haVVVnZ2fD4fCG14IOpzQPUEPikT6qjc/n8/v9oVDI7/ffvn2brqcDf4VCIeNfZ2dnK0qgYymurq5uOM6xpmnVY+dUbFxjcB3jrp977rlyuUyzoFvH4wG29BqSLMsVf6q+RoLBoP5XURQnJibo8mbXQnU5OKVZQULiWltbGx2HKp1O0+GtVFUdHh6mt8iNf/X7/fq1SrfUNC2TySwtLdGV//M//6MXGw6H6U0/Qojf75ckiV51+sYVl72x8OrAJiYmlpaWrBhNHGBD9ISk/1ZcI7FYTBAEPSeFQqFTp04NDw+TjU5v48mPU5oHGMuOC/TUt/+3GL17bvNOAR6KoiihUIh1FGAHJCRPw6UOAPxAQgIAAC6glR0AAHABCQkAALiAhAQAAFxAQgIAAC4gIQEAABeQkAAAgAtISAAAwAUkJAAA4AISEgAAcAEJCQAAuICEBAAAXEBCAgAALiAhAQAAF5CQAACAC0hIAADABSQkAADgAhISAABwAQkJAAC4gIQEAABcQEICAAAuICEBAAAXkJAAAIALSEgAAMAFJCQAAOACEhIAAHABCQkAALiAhOQk58+fN3EzAACuNKyvr7OOwYtaWlroQmNj49LS0mabdXd3z83NGV+1vLy8WYH6n2psZgo9+H379s3NzVXvzuoAAMCVUENiZnl5eXl5eW1t7datW+fPnw8EAoFA4NatW6+++mogEOju7ta3pH+la2gyoP8eP368Ystqeslf//rXCSGnT58+ffq0cXfd3d3d3d1Xrlwxlqbvsbu727hxRfDXrl2r3tHrr7++urp6/Phxk98vcJqVlZU//MM/fPLJJ1tbW6PR6OLi4oab9fb29vb2Vi9XbNPS0vLkk08GAoGXXnopkUikUqnNCnlYMzMzBw8ebG1t7e7uHh0drb8oGlVLS8vevXtTqVT1C3cSlTehhsRGS0vLvn377ty5s2fPng8++CAQCPT398/Ozra0tFy/fn1kZIQQ8vLLL9OqRiAQePfddy9cuEDXLy8v0/WiKL755pvBYHBtbW2zGpJe8sLCwo9+9KNUKjU3N2fc3dzc3JkzZ15++WVjaQcOHHjrrbfoHhsbG/WNZVkmD9Lh/fv3n332Wb2GRIOcnZ194403PvzwQ9SQIB6P37lz58SJE+3t7devX//000/z+fyNGze6urq++93vfv/736fLmqY98sgjXV1d9+/fb25u3r17d3t7++uvv76ysvLSSy+Njo4SQujX+szMTDabvXz58vLycnd396FDh+hm3d3d//zP//zFF1/80R/90U9/+tN79+699NJLhBB6d6G7u7u/vz+VSi0uLvb29r7yyiupVIruOp/PNzc3p1Kpjz/+eGBgYH5+fmlp6ezZs1988cWJEyfm5uboZnv37n3//fcJIZ999hkh5MiRIxVRtbe3E0Lov6Ojo7T8Z5555tq1a7/4xS9+8IMfxONxVp+CszzCOgDvoldLIBAghOzZs+epp57q6el56qmn/vRP//Q///M/z5w509PTQ7fcs2fP008/3dfXRxOSbnl5+Y033jCu6e7u/vDDD41r9JKffvpp/UelcXdzc3Mvv/xyRWk+n0/fo3Fj464rDocGSQj5/ve/v9O3BlxhZWWlsbGRfmuvrKykUqnHH3/83Llz+Xz+z//8zzs6Oujy+++///nnnz/66KMnTpygG/f19T3++OPf/va3Dxw4UFFmPB6fnp7+4IMPHnnkkcXFRX2zO3fuEEJCoVB7e/ulS5dmZmaam5vX19cbGhrm5+dXVlZaW1u//e1v/+///u/Q0NC9e/cOHTq0urqazWZTqdTJkyfPnj375ptvEkKOHj16+PBhQsjNmzf1zd56662vfOUrDQ0Njz32GC1Qj+fatWs+n6+xsfFb3/oWXT8yMvnambYAABYWSURBVKIf5s9//vOvfe1ryEb1Q0JibG1tjRCyd+/eM2fOrK2tXbx4MZVK3blzh66n9u7dS+9X0EpVd3d3Y2MjIeTOnTuzs7OEkFdffZVuSZMcvY1ACPm93/s9Y8kvvvhiNput2J2+l4rSRFHU9169cbW9e/cGAoHl5eVsNnvy5MlXX3319ddfN/GNAidaWVmhX9Nnzpz55JNP2tvb6a3jhoaGJ598ki5//vnnLS0tu3btmp6epgnp5MmTc3NzuVyuqamJVmL0Aufn5xsbG3fv3l2x2RdffPHII48sLy/fuHGjXC7Tl+zatUt/4aOPPtrW1ra0tOT3+xsaGp599llCSFdXFyFkZGTk/v37P/zhD1Op1Ntvv/3ZZ5/t3r27ubl5bW2NbvYv//IvxqKMDh48ODMzQ5dp8HTX+qGZ/pa6G27ZsXH+/HlaL7ly5Qoh5MiRI2zjqXD8+PFnn332v/7rv2ZnZ42tKig9eIAaFhcXjx8//tFHH+3ataurq+s3f/M3L168qKrqU0899cILL1y5coUu+3y+PXv2dHV10ftmvb29qVTq0qVLKysrnZ2dP/rRj5qbm3t7e+kDy127dv3Wb/2WfpdP36y5ufnu3buEkJWVldXV1cbGxl/7tV974oknaCQvvPDCP/zDP9y9e7ezs/N3f/d3L168ePfu3a997Wvnzp1rb2+/cOFCOp1eWlratWvXN7/5TULI0tISvc1INyOE0BSoo0lIv2VHV9L/0n3RQ2ttbb13796rr77a19dn29vuaEhIsIHz58+fPHlyz549b731Fr0RBwAVaOWPVrPAFEhIAADABTT7BgAALiAhAQAAF5CQAACAC0hIAADABSQkAADgAhISAABwAQkJAAC44JKEpCiKJEnFYpF1IABWoSe5oijV65nEA2A6l4xll8vlBgcHx8bGwuEwXdPR0UEXpqam2MUF/NLPEKfQT/JMJqOvlGW5UCjk83l9jX5cCwsLNkcIsEMuSUiZTEaSJDpytg4XJLhJIBAQBEHTNH3NZnUjnPngUC65ZVcoFJLJZKlUYh0IgFXK5TIhRBAEfU2hUFAURVVVdkEBmMklCen27dtDQ0N+v591IABWaWtrSyQSPp+PECJJEiEknU4nk0ljigJwNNcOrtrR0YEbF+BBOPPBuVxSQwIAAKdDQgIAAC4gIQEAABeQkAAAgAtISAAAwAUkJAAA4AISEgAAcAEJCQDcqeEB48rFxcUNl4EHSEgA4Gbr6+t6TpqZmenr66teBk4gIQGAmzU0/HI8mt7e3ubm5upl4AQSEgC4k1vHRXMxJCQAcDOkJQdBQgIAD5mZmdlwGXiAhAQAAFxAQgIAAC4gIQEAABeQkAAAgAtISAAAwAUkJNhU9cgrGHYFAKyDhAQbo/3bjX04MOwKAFjKJQlJURRJkmRZZh2IexhHAKMw7AoAWOoR1gGYQ5blaDSazWZFUdRXdnR0EEKmpqbYxeV4DQ0NN27cWFhYoP9dW1vbcNmJ6OkBAPxwSUJKJpOapvl8PuNKR39d8qB6zJXGxkb9e9y4DACwcw3uGOipWCyqqhqJRPQ1HR0dSEjgQTjzjYxDfQP/XPIMqVAoyLIci8VYBwIAANvkmFt2mqb5/f7N/prJZOwMBgAATOeAGpIkSYqinDp1KpFIsI4FgBnalFRRlIo1xWKRYVQAJnJAQiqXy4SQSCSiaRrrWACYyeVy0Wg0l8tVrCkUCgyjAjCRAxJSKBQqFAqCIPT09LCOBYCZQCAgCILxZ1kmk5mYmAgEAhVbdnR0oAEkOJEDniFFIhHafG5gYIB1LADM0FsFgiDoawqFQjKZrG7Lg1Z24FC8J6SKiy2fz7OJA4C1tra2RCJB60OSJCWTydu3bw8NDdVo7APgLE5qpF8oFIw9jWpDbwzwJpz5RuiH5Cy815AIIZIkXb16lRDi8/nqT0gAAOAsDkhI5XJ5eHi4XC5j7FQAABdzQCs7n88ny/Ls7CzrQMABas/hRAiZn5+3PSgAqIsDakjJZJJ1COAkxokzZmZmUqmUnoQuXLgwPj4+MzPDLrpaNE2bnJwMBoOCIBhb0wF4hANqSLEHMFIDPCzjvE38142y2ezq6qqmaUNDQ6xj8aiK6jV5UMOm61OplP5rBjVvKzggIeUfqJhdAqAarR5t2LBqfHz8vffe43nm9XK5HAqF/H4/WnIzoc+SbKxh9/X10fV//Md/fPr06fHxcVI1Y/KFCxdSqRSboN3FAbfs9K5IwWCQaSC/RM9XNCfl02afSzabJXz/khVFkY4DdPToUdaxeJexhmSsYf/TP/0TIaS9vZ04rebtIFwnJEVRxsbGbNudfiLWzjT6D3B0cXCEiidG3D5AIoQEg8HBwUHyq8MxgG30ulHFdW1c39vbW/Gq8fHx5uZmnmveDsL1LbtQKJTP5+lVmk6n6dAplqqosAPYaWxsTNM0TdOOHTvGOhaP2uwnJl0fj8f1ipEum82Ojo7SmhPsENc1JKpcLguC4Pf7bUhI9aSizX5GAeyQfqrjcSk/jFVqete3en31f2F7HJCQRFGkN+6svrFef4JBKgIrRCIReqpjRBLwJtc+BcGIXuBNOPONtvGgd7OX4JmxDXivIUmSRAgplUr0vxjtG1xMkiSc6uBlvCckfZgGTdPQOQPcrampKZPJ4DwHz+I9IRFCJEkSRZHO3JzJZDbbTFVVNJYFR1tdXdWHI0ENCTzIAQmJNq6LRCLGJi4VisXi2NjYxYsXjSvpLM5TU1NWRwhOxOEk3/r9ANo9FsBrHJCQQqFQoVCIx+M9PT2bbRMOh6tbyuLRLjgLpv4Cj3NAQjp69Gi5XFZVFeOpgLth6i/wOAckJHqnTtO0bDaLG+vgYnTqL9ZRADDjgIRULpfD4fCWQyAjV4HTtbW10Tt1eIYE3uSAhIQhkMELYrFYqVSiNSRBEPAMCTzIAX2P9fbcmqbRwb7qeRX6q4PjKIoSCoV22OXOONAilknV4JNbLleM+r/lemctc/6Fz3sNKZFIaJomCEK5XNY0LZ1Os44IwCq0erRll7va9u/fr/8UM377eHO5ehxkLPOM6+knCCGapuXzeVVVM5lMPp+3uutr9QTGO9kM4KHoXe40TWMdCwADvCckOxln3ttyM0ybBKajXe4EQajR5Q7AxXi/pajPX04NDg6GQqF6XriNZ0h1TgWLGWOBZ3h6aoTRvp2F92dINjfmrmfmPUzQB1bAwPYAXCckRVHofGVUMBjUB/uyAiboA1YURdFTESEkGAxuuI0sy6Io6jcJ6JpQKCSKok2BAliJ62dIoVAon88Hg8HBwcF0Om3DFOYATNRzqudyuWg0StvgUbIsR6PRYrFoY6QAFuK6hkTRvkd+vx8JCdyt9qkeCAQEQTA2wEsmk5qmVQ8rTAcyx5MkcBwHJCRRFOmNO4zUAO5W+1SnWcrY86FYLKqqWn0fG6kIHMoBCSkYDJZKpWAwiPn3wN3C4XA4HN7sr21tbYlEIhAIEEIkSUomk4VCQdM0WZbRAgLcwQENGWnrI9pFo/4LD41fwZtw5huh2bezcN2ogSqXy6FQaMvRvgHcAcM0gGc54JYdRvsGj5AkSRTFHY5lB+BcDqghZbPZcDicyWQwSQy4G8ayA49zQEKijN0vAFwJY9mBxzkgIQWDwUgkQtvasY4FwCqxWKxcLqfTaUEQBgYGWIcDwIBr242grRE4jj41X7FYrNH+uzac+UZoZecsvDdqwIiT4B3ZbLanp2d2dpYQsu2EBOBcvN+ySyaTbW1t+Xw+n8/XGEFSURRJkhRF2eHuMEGfkUMP04kxU8lkkg6Waukgws5VcULOz8+nUqn5+XmGIYG5eE9IsVhsbGwsFovFYrEa+aZ63MltqHPmPY9M0OfQw6RhOytmip7kqqoWCoWKacCAbHRCjo+P9/f3j4yMsA0MTMT7Lbt8Pq8oypaT8lWPO3nz5s1tfCvpL9kyJ9WzmdM59DBptFvGzNvzgHw+r6rqxMREuVyORCKsw+FUxce6srKCGpKb8J6QCCGyLOuzIm32DKl63Mn9+/dv49FunTPveWSCPv3id9CROjFm3djY2MDAgN/vHxoawuPSCtUTY/b39xNC+vr6WIYFpnJAQmpqaspkMrXHDTKOO7kTdX6LOfHLbhuceJhOjFnn8/no9BMYJWtDFTdjp6enFxcXDx06xDAkMJcDGjJKkrSNVnZo/AqOUywWZVkmhITD4W1PAuvuM/9h216j2bezOKCGFI/HJycnMf0EuJ4gCD6fj07TxzoWAAZ4b2VHCMlms6urq5qmDQ0NsY4FwEJjY2OiKMbjcf2hKYCnOCAhYfoJ8Aj6DIk+RmIdCwADDrhlh+knwCNqT2EO4Hpc15AURYnFYtlsVtM0TdN2PhADALdUVSWERCIRn8+HcYPAm7hOSKFQKJ/PB4PBwcHBdDpNOxsBuNLQ0FAwGAyFQqIonjp1inU4AAw44JYdbXTk9/uRkMDF9KekW45LAuBWDkhIuLEOXkBvSodCIdyaBs9ybVcvd3cPBPdRFIWODuzz+eLx+La7Irn7zEfHWHdzQA1JH/nY7/dnMhmmsQBYJRQK4fQGj3NAQtKHC0LHWAAAF3NAQtJrSMFg0Op9YbRvAABWHJCQRFGk08PQ7rHW0e8R175ZXOdmLoC8CwB24j0hxWKxUqlEh0AWBAETl9mMDviPnAQANuA9IekzxtYzb+wOVc8AtpPNnA55iDeKosiyLIqi8UJQVRVDg4NrcD1SAyFEkqSrV6/mcrlCoWBDo4b19fV6voXr3MzRaN511uTl7pbL5aLRKG0aThWLxUQiUb1lR0dHR0eHjaEBmIP3hLS0tDQwMODz+dLpNB3sC2yz/gDrQIAQQgKBgCAImqbpa8LhsM/nq95yYWHBxV2RwMV4T0iaphWLxXK5rGkaxuQHL6NDZ+EGHbgY7wkpEonMzs729/dns1kMHQRe1tbWlkgkaJVIkiTW4QCYz7UPrt09gArAZtx95mPoIHfjvYZUPzxhAgBwNJckpM2aGwEAgFPw3g+pThs2N6ItX6emplhEBLxDw2gA3jg+IcViMX1soQouvpMOAOA+jk9I+ljgAACc03uaW9c+woZdWMclz5AIMhMAOIQ+/Jijd2EF9yQkAABHsCFVOLSROhISAJhgGyMfbvYSFw+iqCcJ67KFE/OQDgnpV9R5Jbj4ggHYtoe6TVS7/6lD7zjVyYac4dC0hIT0S3VeCV64YAAeCr0otlFD2vBVuLI8CwkJ6jU/P59Kpebn51kHAtzZ9u+z6nqSR6Zjhg0hIf2SPgNQPRP0eeGCaXiA/nd8fLy/v39kZIRtVMAnejnUf1HU3t71FxdsyPH9kMxV52XghatFz7jGX74rKyuoIQGARVBDgo1Vzxjb399PCOnr62MX1NYqanWEkMXFRbpAbznOzMwwCg0AtoCEBJuqmDF2enp6fHy8paWFbVRbMtZfZ2Zm9Aw6MjLyyiuvjI+PM4oLeFH9qwU4gVt2UK/R0VHWITy03t7e5uZmunzhwoVUKtXe3s42JOAEvQfghdvvDoIaErhKjW+ZbDY7OjqKZ2BAIRtxCAkJ3Gazb5mf/exn8XhcrzCBZz1sg0CwjWt/I7h7ImeAzTA8882aX7xGOQx3YVaNyoaamXMrf6ghAQAAF5CQAACAC2hlBwBgsm3Mkkdf8rDbW7oL+6GGBOB+m/W82UaPnO1t781OP/UP8bfhwCjm0rMRz58FEhKA13H1vekmzm1cwAoSEoDX4XvTdA/7fuq/CcwanXaz7Tn/rJGQALyi+puIw+9NN3moo7bnLeL8g0BCMkFHR4dbd4d9ucOWM6qYVRQrDD9Qr51LlnJJKztFUWRZDoVCoiiyjgXAJvS0F0UxFAqZW/I2WnAxRKNlcjPKhl0767PYIa7vJ9ZPkqRoNJrNZtPpNF2Dny1QmwsG8kgkEoODg2NjY5lMRl+J5gZQA+df+C6pISWTSU3TfD6fvsYFXzcAtQUCAUEQNE0zrjRxeBv+u61QerRMakhW75rh0dnP8QkpFouJoigIgqqqyWSSdTgA9imXy4QQQRBML9lB2Yhi+H2NselM5JLjTCQS9HdiPp9nHQuATXK5nKIogUAAP8XAHVySkAAAwOnQ7BsAALjw5e9973usYzCNoihnz571+/1PPPGEvlJVVb/fb9teVFUtlUqaphn/avoe6Zpyubxv3z579mLncTU0NJj7XKT2jqw4NCY2PDPtweo91K87+4/d0ku+NuOFyfBDt4Kraki5XC4ajeZyOX1NsVhMJBJ27mViYsLc3W24R1mWo9FosVi0bS/2HBddUygU7NyRFYfGRPVh2obJe2i87mw+dqsv+dqMFybDD90KrkpI1a1gw+GwsS24PXuRZfnq1auW7jGZTPp8PnMPbcu92HBcmUxmYmIiEAiYuJd6dmT6oTGxYStw29j/HhqvO5uP3epLvjbjhcn2QzedqxKSda1g699LNBpNp9Ozs7OW7rFYLF6+fNnctlW192LPcRUKhWQyWSqVTNzLljuy4tCYsOf83xDz99BTx268MBkeuBVclZDa2toSiQT94SBJEpO9TExMDA0NmXt+VO+xUCjIshyLxWzbiz3Hdfv27aGhIXOf+W25IysOjQnjYdqM+XvoqWM3XpgMD9wKaPYNAABccFUNCQAAnAsJCQAAuICEBAAAXEBCAgAALiAhAQAAF5CQHExRFNYhADBg7pmP64gfSEgONjY2xjoEAMtpmpZIJGKxWCKRoMnD3DMf1xE/HD9BH1CKotDxrIaHh7PZbDKZpENd+Xw+43pCyJ49ez744AM6iQ59FR1BRxRFfUvTu6YCbNvY2FgkEgmHw6qq6vO1S5JUKpXS6bSiKJcvX66eFEqSpKamptXV1WQyGYvF8vk8/VdfTzejr9JL0zSt4nqJRqMu6DTtFKghuUQulzt69GggEJicnCSEqKo6OzsbDoeN60ulkiiKX/7yl4eHh+nQW4VCYXh4uKenp1QqVZQAwIlyuRwOhwkhxkHbRFEcHBycmJi4fPnywMCAKIoVryqVSv39/dXDUOnr4/G4/le9tOrrBdnITkhILuHz+fx+fygUCgaDoijqIxAb1xNCQqHQ6uqqPkw4/Sv9U8WWAJzw+Xz0Tp2maXrdXRAEQRDK5fLAwEC5XM7lctUDjG5W0afrjX/VS6u+Xqw4ItgMbtk5Gx3PShTFtra2QqGgaVo6nRYE4dSpU8PDw4QQ43r6kqWlJbogyzIhRB+4vnpLAB7E43H6mEdVVXpWB4PBU6dOaZo2ODgoy/LS0hLNRpIkbTjisN/vlyRps/xkLO3q1au4ChjCWHaeVigU6G0QRVHMHTscAOBhoYbkacFgMJfL+Xy+wcFB1rEAgNehhgQAAFxAowYAAOACEhIAAHABCQkAALiAhAQAAFxAQgIAAC4gIQEAABeQkAAAgAtISAAAwAUkJAAA4AISEgAAcAEJCQAAuICEBAAAXEBCAgAALiAhAQAAF5CQAACAC0hIAADABSQkAADgAhISAABwAQkJAAC4gIQEAABcQEICAAAu/D8wADEjQXX0YgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Load data from comma-separated-value (CSV) file with header\n",
    "# Formatting: each %d, %s and %f corresponds to an integer, a string and a\n",
    "# number (double precision) respectively\n",
    "[degree, gender, salary] = textread ('data2.csv', '%d %s %f', 'Delimiter', ',', 'HeaderLines', 1);\n",
    "\n",
    "# ANOVA using wild bootstrap\n",
    "# https://gnu-octave.github.io/statistics-resampling/function/bootlm.html#7\n",
    "[STATS, BOOTSTAT, AOVSTAT] = bootlm (salary, {degree, gender}, 'model', ...\n",
    "                             'full', 'display', 'off', 'varnames', ...\n",
    "                             {'degree', 'gender'}, 'seed', 1);\n",
    "\n",
    "# Report ANOVA results\n",
    "fprintf ('ANOVA SUMMARY with effect of gender averaged over levels of degree\\n')\n",
    "for i = 1:numel(AOVSTAT.MODEL)\n",
    "  fprintf ('F(%u,%u) = %.2f, p = %.3g for the model: %s\\n', ...\n",
    "           AOVSTAT.DF(i), AOVSTAT.DFE, AOVSTAT.F(i), ...\n",
    "           AOVSTAT.PVAL(i), AOVSTAT.MODEL{i});\n",
    "end\n",
    "\n",
    "# Calculate regression coefficients and 95% credible intervals (using bayesian bootstrap\n",
    "# with an uninformative prior)\n",
    "MDL = bootlm (salary, {degree, gender}, 'model', 'full', ...\n",
    "                        'display', 'on', 'varnames', ...\n",
    "                        {'degree', 'gender'}, 'method', 'bayesian', 'prior', 1);\n",
    "\n",
    "# Calculate estimated marginal means and 95% credible intervals (using bayesian bootstrap\n",
    "# with an uninformative prior)\n",
    "EMMEANS = bootlm (salary, {degree, gender}, 'model', 'full', ...\n",
    "                        'display', 'on', 'varnames', ...\n",
    "                        {'degree', 'gender'}, 'dim', [1, 2], ...\n",
    "                        'method', 'bayesian', 'prior', 'auto');\n",
    "\n",
    "# Test whether salaries averaged over levels of college differ between men and \n",
    "# women (using wild bootstrap)\n",
    "POSTTESTS = bootlm (salary, {degree, gender}, 'model', 'full', ...\n",
    "                        'display', 'on', 'varnames', ...\n",
    "                        {'degree', 'gender'}, 'dim', 2, ...\n",
    "                        'posthoc', [1, 2]);\n",
    "\n",
    "# Save the figure of the diagnostic plots\n",
    "#print ('diagnostic_plots.svg','-dsvg')\n",
    "#print ('diagnostic_plots.png','-dpng','-r300')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example 3** (below): Bootstrap two-way ANOVA (Type I) using `bootlm` with data (in long format) as above but loaded from **data.xlsx**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ANOVA SUMMARY with effect of gender averaged over levels of degree\n",
      "F(1,18) = 87.20, p = 0.0001 for the model: salary ~ 1 + degree\n",
      "F(1,18) = 10.97, p = 0.00498 for the model: salary ~ 1 + degree + gender\n",
      "F(1,18) = 0.42, p = 0.572 for the model: salary ~ 1 + degree + gender + degree:gender\n"
     ]
    }
   ],
   "source": [
    "# Load data from Microsoft Excel or Open Office document spreadsheet\n",
    "# First output is for numeric data, second output is for text data\n",
    "# https://octave.sourceforge.io/io/function/xlsread.html\n",
    "[degree, gender] = xlsread ('data.xlsx', 'data2', 'A2:B23');\n",
    "salary = xlsread ('data.xlsx', 'data2', 'C2:C23');\n",
    "\n",
    "# As above, ANOVA using wild bootstrap\n",
    "[STATS, BOOTSTAT, AOVSTAT] = bootlm (salary, {degree, gender}, 'model', ...\n",
    "                             'full', 'display', 'off', 'varnames', ...\n",
    "                             {'degree', 'gender'}, 'seed', 1);\n",
    "\n",
    "# Report ANOVA results\n",
    "fprintf ('ANOVA SUMMARY with effect of gender averaged over levels of degree\\n')\n",
    "for i = 1:numel(AOVSTAT.MODEL)\n",
    "  fprintf ('F(%u,%u) = %.2f, p = %.3g for the model: %s\\n', ...\n",
    "           AOVSTAT.DF(i), AOVSTAT.DFE, AOVSTAT.F(i), ...\n",
    "           AOVSTAT.PVAL(i), AOVSTAT.MODEL{i});\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "**Example 4** (below): Bootstrap confidence intervals using `bootci` or `bootknife` with data specified in the code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------------------------------------------------------------------------------\n",
      "|          | norm     | perc     | basic    | bca      | stud     | cal      |\n",
      "------------------------------------------------------------------------------\n",
      "| CI_lower | 108.4481 | 95.83105 | 106.7427 | 115.6061 | 111.1151 | 113.3279 |\n",
      "| CI_upper | 247.9999 | 236.3254 | 247.237  | 262.1344 | 299.1946 | 281.9322 |\n",
      "------------------------------------------------------------------------------\n",
      "\n",
      "Summary of nonparametric bootstrap estimates of bias and precision\n",
      "******************************************************************************\n",
      "\n",
      "Bootstrap settings: \n",
      " Function: var\n",
      " Resampling method: Balanced, bootknife resampling \n",
      " Number of resamples (outer): 4999 \n",
      " Number of resamples (inner): 0 \n",
      " Confidence interval (CI) type: Bias-corrected and accelerated (BCa) \n",
      " Nominal coverage (and the percentiles used): 90% (11.9%, 98.7%)\n",
      "\n",
      "Bootstrap Statistics: \n",
      " original     bias         std_error    CI_lower     CI_upper  \n",
      " +171.5       -6.773       +42.49       +114.7       +262.4     \n",
      "\n",
      "\n",
      "Summary of nonparametric bootstrap estimates of bias and precision\n",
      "******************************************************************************\n",
      "\n",
      "Bootstrap settings: \n",
      " Function: var\n",
      " Resampling method: Iterated, balanced, bootknife resampling \n",
      " Number of resamples (outer): 4999 \n",
      " Number of resamples (inner): 499 \n",
      " Confidence interval (CI) type: Calibrated percentile\n",
      " Nominal coverage (and the percentiles used): 90% (11.0%, 99.8%)\n",
      "\n",
      "Bootstrap Statistics: \n",
      " original     bias         std_error    CI_lower     CI_upper  \n",
      " +171.5       -7.047       +45.72       +111.2       +298.6     \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Spatial Test Data from Table 14.1 of Efron and Tibshirani (1993)\n",
    "# An Introduction to the Bootstrap in Monographs on Statistics and Applied \n",
    "# Probability 57 (Springer)\n",
    "# AIM:\n",
    "# To construct 90% nonparametric bootstrap confidence intervals for var(A,1)\n",
    "# var(A,1) = 171.5\n",
    "# Exact intervals based on Normal theory are [118.4, 305.2]\n",
    "A = [48 36 20 29 42 42 20 42 22 41 45 14 6 ...\n",
    "    0 33 28 34 4 32 24 47 41 24 26 30 41].';\n",
    "ci1 = bootci (4999,{{@var,1},A},'alpha',0.1,'type','norm','seed',1);\n",
    "ci2 = bootci (4999,{{@var,1},A},'alpha',0.1,'type','per','seed',1);\n",
    "ci3 = bootci (4999,{{@var,1},A},'alpha',0.1,'type','basic','seed',1);\n",
    "ci4 = bootci (4999,{{@var,1},A},'alpha',0.1,'type','bca','seed',1);\n",
    "ci5 = bootci (4999,{{@var,1},A},'alpha',0.1,'type','stud','nbootstd',100,'seed',1);\n",
    "ci6 = bootci (4999,{{@var,1},A},'alpha',0.1,'type','cal','nbootcal',499,'seed',1);\n",
    "\n",
    "# Create table of confidence intervals using bootci\n",
    "header = {'','norm', 'perc', 'basic', 'bca', 'stud', 'cal'};\n",
    "tbl = table ({'CI_lower'; 'CI_upper'}, ci1, ci2, ci3, ci4, ci5, ci6);\n",
    "tbl = setVariableNames (tbl, (1:7), header);\n",
    "prettyprint (tbl(:,:))\n",
    "\n",
    "# Bias-corrected and accelerated bootstrap confidence intervals after balanced bootknife resampling \n",
    "bootknife (A,4999,{@var,1},[0.05,0.95]);\n",
    "\n",
    "# Calibrated bootstrap confidence intervals after iterated balanced bootknife resampling \n",
    "bootknife (A,[4999,499],{@var,1},[0.05,0.95]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "8.3.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}