Mixed Model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-3-7469fad3f73d>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m    117\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    118\u001b[0m \u001b[0mbiogeme\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmodelName\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'17lognormalMixture'\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 119\u001b[1;33m \u001b[0mresults\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbiogeme\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mestimate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    120\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresults\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\biogeme\\biogeme.py\u001b[0m in \u001b[0;36mestimate\u001b[1;34m(self, bootstrap)\u001b[0m\n\u001b[0;32m    362\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    363\u001b[0m         \u001b[0mopts\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m{\u001b[0m\u001b[1;34m'gtol'\u001b[0m \u001b[1;33m:\u001b[0m \u001b[1;36m1e-7\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'ftol'\u001b[0m \u001b[1;33m:\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfinfo\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfloat\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0meps\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 364\u001b[1;33m         \u001b[0mresults\u001b[0m \u001b[1;33m=\u001b[0m  \u001b[0mminimize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf_and_grad\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbetaInitValues\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mbounds\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbounds\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mjac\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0moptions\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mopts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    365\u001b[0m         \u001b[1;31m#        yep.stop()\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    366\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\_minimize.py\u001b[0m in \u001b[0;36mminimize\u001b[1;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[0;32m    601\u001b[0m     \u001b[1;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m'l-bfgs-b'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    602\u001b[0m         return _minimize_lbfgsb(fun, x0, args, jac, bounds,\n\u001b[1;32m--> 603\u001b[1;33m                                 callback=callback, **options)\n\u001b[0m\u001b[0;32m    604\u001b[0m     \u001b[1;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m'tnc'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    605\u001b[0m         return _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\lbfgsb.py\u001b[0m in \u001b[0;36m_minimize_lbfgsb\u001b[1;34m(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, **unknown_options)\u001b[0m\n\u001b[0;32m    333\u001b[0m             \u001b[1;31m# until the completion of the current minimization iteration.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    334\u001b[0m             \u001b[1;31m# Overwrite f and g:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 335\u001b[1;33m             \u001b[0mf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfunc_and_grad\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    336\u001b[0m         \u001b[1;32melif\u001b[0m \u001b[0mtask_str\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstartswith\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34mb'NEW_X'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    337\u001b[0m             \u001b[1;31m# new iteration\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\lbfgsb.py\u001b[0m in \u001b[0;36mfunc_and_grad\u001b[1;34m(x)\u001b[0m\n\u001b[0;32m    283\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    284\u001b[0m         \u001b[1;32mdef\u001b[0m \u001b[0mfunc_and_grad\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 285\u001b[1;33m             \u001b[0mf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    286\u001b[0m             \u001b[0mg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mjac\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    287\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mg\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\optimize.py\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[1;34m(*wrapper_args)\u001b[0m\n\u001b[0;32m    291\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    292\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 293\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    294\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    295\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\optimize.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, x, *args)\u001b[0m\n\u001b[0;32m     61\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     62\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnumpy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 63\u001b[1;33m         \u001b[0mfg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     64\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mjac\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     65\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mfg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\biogeme\\biogeme.py\u001b[0m in \u001b[0;36mf_and_grad\u001b[1;34m(x)\u001b[0m\n\u001b[0;32m    354\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    355\u001b[0m         \u001b[1;32mdef\u001b[0m \u001b[0mf_and_grad\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 356\u001b[1;33m             \u001b[0mf\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mg\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mh\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcalculateLikelihoodAndDerivatives\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mhessian\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mbhhh\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    357\u001b[0m             \u001b[0mscale\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minitLogLike\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[1;36m100.0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    358\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[1;33m-\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mscale\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mg\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mscale\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\biogeme\\biogeme.py\u001b[0m in \u001b[0;36mcalculateLikelihoodAndDerivatives\u001b[1;34m(self, x, hessian, bhhh)\u001b[0m\n\u001b[0;32m    279\u001b[0m                                                                \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbetaIds\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    280\u001b[0m                                                                \u001b[0mhessian\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 281\u001b[1;33m                                                                bhhh)\n\u001b[0m\u001b[0;32m    282\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    283\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mh\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbh\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import biogeme.database as db\n",
    "import biogeme.biogeme as bio\n",
    "import biogeme.draws as draws\n",
    "import biogeme.models as models\n",
    "\n",
    "pandas = pd.read_table(\"tripdata.dat\")\n",
    "database = db.Database(\"tripdata\",pandas)\n",
    "\n",
    "# The Pandas data structure is available as database.data. Use all the\n",
    "# Pandas functions to invesigate the database\n",
    "#print(database.data.describe())\n",
    "\n",
    "from headers import *\n",
    "\n",
    "# Removing some observations can be done directly using pandas.\n",
    "#remove = (((database.data.PURPOSE != 1) & (database.data.PURPOSE != 3)) | (database.data.CHOICE == 0))\n",
    "#database.data.drop(database.data[remove].index,inplace=True)\n",
    "\n",
    "# Here we use the \"biogeme\" way for backward compatibility\n",
    "#exclude = (( PURPOSE != 1 ) * (  PURPOSE   !=  3  ) +  ( CHOICE == 0 )) > 0\n",
    "#database.remove(exclude)\n",
    "\n",
    "# Alternative specific constants\n",
    "ASC_AUTO_DRIVE = Beta('ASC_AUTO_DRIVE',0,None,None,0)\n",
    "ASC_AUTO_PASS = Beta('ASC_AUTO_PASS',0,None ,None,0)\n",
    "ASC_METRO = Beta('ASC_METRO',0,None,None,0)\n",
    "ASC_TRAIN = Beta('ASC_TRAIN',0,None,None,0)\n",
    "ASC_WALK = Beta('ASC_WALK',0,None,None,1)\n",
    "\n",
    "# Sociodemo variables\n",
    "B_SEX = Beta('B_SEX',0,None,None,0) # 1 if female; 0 if otherwise\n",
    "B_AGE = Beta('B_AGE',0,None,None,0) # age\n",
    "B_EMP = Beta('B_EMP',0,None,None,0) # parttime or fulltime\n",
    "B_WKWK = Beta('B_WKWK',0,None,None,0) # work week\n",
    "B_INC = Beta('B_INC',0,None,None,0) # income\n",
    "B_NVEH = Beta('B_NVEH',0,None,None,0) # number of household vehicles\n",
    "\n",
    "# Generalized travel variables\n",
    "B_IVTT = Beta('B_IVTT',0,None,None,0) # in vehicle travel time\n",
    "B_DIST = Beta('B_DIST',0,None,None,0) # distance\n",
    "\n",
    "# Car-related variables\n",
    "B_DLIC = Beta('B_DLIC',0,None,None,0) # drivers license 1 or 0\n",
    "\n",
    "# Train-related variables\n",
    "B_NTERM = Beta('B_NTERM',0,None,None,0) # commuter train access station\n",
    "B_GFARE = Beta('B_GFARE',0,None,None,0) # commuter train fare\n",
    "B_GWAIT = Beta('B_GWAIT',0,None,None,0) # commuter train wait time\n",
    "B_GACCA = Beta('B_GACCA',0,None,None,0) # commuter train wait time\n",
    "B_EGRESS = Beta('B_GACCA',0,None,None,0) # commuter train wait time\n",
    "\n",
    "# Metro-related variables\n",
    "B_MSUB = Beta('B_MSUB',0,None,None,0) # closest metro station to destination\n",
    "B_MDIS = Beta('B_MDIS',0,None,None,0) # distance from metro to destination\n",
    "B_TFARE = Beta('B_TFARE',0,None,None,0) # transit fare\n",
    "B_TWAIT = Beta('B_TWAIT',0,None,None,0) # transit wait time\n",
    "\n",
    "B_IVTT_S = Beta('B_IVTT_S',0,None,None,0)\n",
    "B_DIST_S = Beta('B_DIST_S',0,None,None,0)\n",
    "B_NTERM_S = Beta('B_NTERM_S',0,None,None,0)\n",
    "\n",
    "# Define a random parameter, log normally distributed, designed to be used\n",
    "# for Monte-Carlo simulation\n",
    "B_IVTT_RND = (B_IVTT + B_IVTT_S * bioDraws('B_IVTT_RND','NORMAL'))\n",
    "B_DIST_RND = (B_DIST + B_DIST_S * bioDraws('B_DIST_RND','NORMAL'))\n",
    "B_NTERM_RND = (B_NTERM + B_NTERM_S * bioDraws('B_NTERM_RND','NORMAL'))\n",
    "\n",
    "# Utility functions\n",
    "\n",
    "# Updating variables to account for 'missing' values \n",
    "nveh_fix =  DefineVariable('nveh_fix', nveh * (nveh != 9), database)\n",
    "gfare_fix = DefineVariable('gfare_fix', gfare * (gfare != 9999), database)\n",
    "tfare_fix = DefineVariable('tfare_fix', tfare * (tfare != 9999), database)\n",
    "gwait_fix = DefineVariable('gwait_fix', gwait * (gwait != 9999), database)\n",
    "twait_fix = DefineVariable('twait_fix', twait * (twait != 9999), database)\n",
    "\n",
    "V1 = ASC_AUTO_DRIVE + \\\n",
    "     B_DIST_RND * trdist + \\\n",
    "     B_DLIC * dlic + \\\n",
    "     B_IVTT_RND * ivtt1\n",
    "V2 = ASC_AUTO_PASS + \\\n",
    "     B_NVEH * nveh + \\\n",
    "     B_IVTT_RND * ivtt2\n",
    "V3 = ASC_METRO + \\\n",
    "     B_TFARE * tfare_fix + \\\n",
    "     B_IVTT_RND * ivtt3\n",
    "V4 = ASC_TRAIN + \\\n",
    "     B_DIST * trdist + \\\n",
    "     B_GFARE * gfare_fix + \\\n",
    "     B_NTERM_RND * nterm + \\\n",
    "     B_IVTT_RND * ivtt4\n",
    "V5 = ASC_WALK + \\\n",
    "     B_SEX * female + \\\n",
    "     B_IVTT_RND * ivtt5\n",
    "\n",
    "# Associate utility functions with the numbering of alternatives\n",
    "V = {1: V1,\n",
    "     2: V2,\n",
    "     3: V3,\n",
    "     4: V4,\n",
    "     5: V5}\n",
    "\n",
    "# Associate the availability conditions with the alternatives\n",
    "\n",
    "av = {1: avail1,\n",
    "      2: avail2,\n",
    "      3: avail3,\n",
    "      4: avail4,\n",
    "      5: avail5}\n",
    "\n",
    "# The choice model is a logit, with availability conditions\n",
    "prob = models.logit(V,av,choice)\n",
    "logprob = log(MonteCarlo(prob))\n",
    "\n",
    "biogeme = bio.BIOGEME(database,logprob,numberOfDraws=1000)\n",
    "\n",
    "biogeme.modelName = '17lognormalMixture'\n",
    "results = biogeme.estimate()\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}