{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Using TensorFlow backend.\n" ] } ], "source": [ "import sys\n", "import itertools\n", "from keras.layers import Input, Dense, Reshape, Flatten\n", "from keras import layers, initializers\n", "from keras.models import Model, load_model\n", "import keras.backend as K\n", "import tensorflow as tf\n", "import numpy as np\n", "from seqtools import SequenceTools as ST\n", "from gfp_gp import SequenceGP\n", "from util import AA, AA_IDX\n", "from util import build_vae\n", "from sklearn.model_selection import train_test_split, ShuffleSplit\n", "from keras.callbacks import EarlyStopping\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "from sklearn.gaussian_process import GaussianProcessRegressor\n", "from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C\n", "import scipy.stats\n", "from scipy.stats import norm\n", "from scipy.optimize import minimize\n", "from keras.utils.generic_utils import get_custom_objects\n", "from util import one_hot_encode_aa, partition_data, get_balaji_predictions, get_samples, get_argmax\n", "from util import convert_idx_array_to_aas, build_pred_vae_model, get_experimental_X_y\n", "from util import get_gfp_X_y_aa\n", "from losses import neg_log_likelihood\n", "import json\n", "\n", "import isolearn.io as isoio\n", "import isolearn.keras as isol\n", "\n", "import warnings\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "from keras.backend.tensorflow_backend import set_session\n", "\n", "def contain_tf_gpu_mem_usage() :\n", " config = tf.ConfigProto()\n", " config.gpu_options.allow_growth = True\n", " sess = tf.Session(config=config)\n", " set_session(sess)\n", "\n", "contain_tf_gpu_mem_usage()\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "\n", "def get_z_sample_numpy(z_mean, z_log_var, n_samples=1) :\n", " \n", " n = z_mean.shape[0]\n", " m = z_mean.shape[2]\n", " \n", " epsilon = np.random.normal(loc=0., scale=1., size=(n, n_samples, m))\n", " \n", " return z_mean + np.exp(0.5 * z_log_var) * epsilon\n", "\n", "#Evaluate VAE Likelihood (ELBO) on supplied data\n", "def evaluate_elbo(vae_encoder_model, vae_decoder_model, sequence_one_hots, pwm_start=0, pwm_end=-1, n_samples=1, decoded_pwm_eps=1e-6) :\n", " _epsilon = 10**-6\n", " \n", " if pwm_end == -1 :\n", " pwm_end = sequence_one_hots.shape[2]\n", " \n", " #Get sequence VAE encodings\n", " z_mean, z_log_var = vae_encoder_model.predict(x=sequence_one_hots, batch_size=32, verbose=False)\n", "\n", " z_mean = np.tile(np.expand_dims(z_mean, axis=1), (1, n_samples, 1))\n", " z_log_var = np.tile(np.expand_dims(z_log_var, axis=1), (1, n_samples, 1))\n", " z = get_z_sample_numpy(z_mean, z_log_var, n_samples=n_samples)\n", " \n", " #Get re-decoded sequence PWMs\n", " decoded_pwms = np.zeros((sequence_one_hots.shape[0], n_samples) + sequence_one_hots.shape[1:])\n", "\n", " for sample_ix in range(n_samples) :\n", " decoded_pwms[:, sample_ix, :, :] = vae_decoder_model.predict(x=z[:, sample_ix, :], batch_size=32, verbose=False)\n", "\n", " decoded_pwms = np.clip(decoded_pwms, decoded_pwm_eps, 1. - decoded_pwm_eps)\n", " \n", " sequence_one_hots_expanded = np.tile(np.expand_dims(sequence_one_hots, axis=1), (1, n_samples, 1, 1))\n", " \n", " #Calculate reconstruction log prob\n", " log_p_x_given_z = np.sum(np.sum(sequence_one_hots_expanded[:, :, pwm_start:pwm_end, :] * np.log(np.clip(decoded_pwms[:, :, pwm_start:pwm_end, :], _epsilon, 1. - _epsilon)) / np.log(10.), axis=3), axis=2)\n", "\n", " #Calculate standard normal and importance log probs\n", " log_p_std_normal = np.sum(norm.logpdf(z, 0., 1.) / np.log(10.), axis=-1)\n", " log_p_importance = np.sum(norm.logpdf(z, z_mean, np.sqrt(np.exp(z_log_var))) / np.log(10.), axis=-1)\n", "\n", " #Calculate per-sample ELBO\n", " log_p_vae = log_p_x_given_z + log_p_std_normal - log_p_importance\n", " log_p_vae_div_n = log_p_vae - np.log(n_samples) / np.log(10.)\n", "\n", " #Calculate mean ELBO across samples (log-sum-exp trick)\n", " max_log_p_vae = np.max(log_p_vae_div_n, axis=-1)\n", " \n", " log_mean_p_vae = max_log_p_vae + np.log(np.sum(10**(log_p_vae_div_n - np.expand_dims(max_log_p_vae, axis=-1)), axis=-1)) / np.log(10.)\n", " mean_log_p_vae = np.mean(log_mean_p_vae)\n", " \n", " return log_mean_p_vae, mean_log_p_vae, log_p_vae\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "_5k_1\n", "WARNING:tensorflow:From /home/ubuntu/anaconda3/envs/tensorflow_p36/lib/python3.6/site-packages/tensorflow/python/framework/op_def_library.py:263: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n", "Instructions for updating:\n", "Colocations handled automatically by placer.\n", "mean log(likelihood) = -0.8908\n", "mode log(likelihood) = -0.1512\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcXFWd9/HPF5BNcBQIIUFjUMcRhEdIt4rIrmHYogwuqFGMD0kHUFwYZURwSITghqwaocMaDMqIjz4GQYOyyKaQjguRZUASUBJCIozKGsDf/HFuQ6VS1X073ffWre7v+/W6r6q62+/kvtL1q3PuuecoIjAzM6ua9VpdADMzs0acoMzMrJKcoMzMrJKcoMzMrJKcoMzMrJKcoMzMrJKcoMzMrJKcoMzMrJKcoMzMrJI2aHUBirLVVlvF+PHjW10MMytCT0967egoNMyynmUAjO0YW2ickaanp2dVRIzqb79hm6DGjx/PwoULW10MMytCd3d67eoqNExPd0qEHV3FJsKRRtIDefYbtgnKzIaxghNTLyem1vI9KDMzqyQnKDNrP93dLzbzFainu+eFZj4rn5v4zKz9TJ+eXgtu6rty+pWAm/paxTUoMzOrJCcoMzOrJCcoMzOrJCcoMzOrJCcoM7Makl5Y+lpnxXOCMjOzSnI3czNrPxHFh5gHMAOAkwBNLjyk1XENyszMKskJyszMKskJyszaT0dH4VNtAHSf0EX3CeUMTGtr8z0oM2s/ixaVEmb5Us8D1UquQZmZWSU5QZmZWSU5QZmZWSU5QZmZWSU5QZmZWSW5F5+ZtZ9p00oJM2Efz6bbSk5QZtZ+SpjuHWDS1PmlxLHG3MRnZmaV5ARlZu2npyctBVu2ZAzLlowpPI415iY+M2s/nZ3pteBRzeecOB2Ak+bNKDSONeYalJmZVZITlJmZVZITlJmZVVKpCUrSnpJ+LOkhSSFpSo5jdpJ0g6SnsuP+U5JKKK6ZmbVQ2TWozYDFwKeAp/rbWdLLgGuAFcCbs+M+BxxbYBnNzKwCSu3FFxFXAVcBSLo4xyGTgU2Bj0bEU8BiSW8AjpV0ekTBXXjMzKxlqt7N/G3AjVly6vUz4GRgPLCkFYUysxZbuLCUMNNOOa+UONZY1RPUNsCf69atqNm2RoKS1AV0AYwbN67wwplZi5Qw3TvA2O2WlxLHGhtWvfgiojsiOiOic9SoUa0ujpmZDULVE9TDwOi6daNrtpnZSNTVlZaCzT9/EvPPn1R4HGus6gnqVmAPSRvXrJsILAOWtqREZtZ6c+akpWCLrutg0XXlNCfa2sp+DmozSTtL2jmLPS77PC7b/mVJv6g55DLgSeBiSTtKOhT4POAefGZmw1zZNahO4DfZsgkwM3v/pWz7GOC1vTtHxF9JNaaxwELgW8A3gNPLK7KZmbVC2c9BXQ80HQUiIqY0WHcHsGdxpTIzsyqq+j0oMzMboZygzMyskqr+oK6Z2domTCglzJjxy0qJY405QZlZ+xnkdO95J0TomtU9qDg2OG7iMzOzSnINysxGrJi39jpNLr8c1phrUGbWfqS0FGzm5BnMnDyj8DjWmBOUmZlVkhOUmZlVkhOUmZlVkhOUmZlVkhOUmZlVkhOUmZlVkp+DMrP2c955pYQ5+Ij5pcSxxpygzKz9lDDdO0DHvoMbUskGx018ZmZWSU5QZtZ+urvTUrCeazvoubaj8DjWmJv4zKz9TJ+eXgtu6rvygkmAm/paxTUoMzOrJCcoMzOrJCcoMzOrJCcoMzOrJCcoMzOrJCcoMzOrJHczN7P2E1FKmJPmzSgljjXmGpSZmVVSrgQlaVTRBTEzM6uVtwb1kKQrJB0gSYMJKOloSUskPS2pR9Ie/ez/IUm/lfSkpIclfUfSNoMpg5m1uY6OtBSs+4Quuk8oZ2BaW1veBHUQsBr4AfCgpJMlvXagwSQdBpwFnArsAtwCXC1pXJP93w5cClwCvBE4BNgBmDfQ2GY2jCxalJaCLV86luVLxxYexxrLlaAi4pqI+BAwFvgKcADw35KulTRZ0sY54x0LXBwRcyLirog4BlgOHNVk/7cBf46IMyJiSUT8CjgHeGvOeGZm1qYG1EkiIv4nIr4VEZ3AJ4HdSDWcZZK+ImmzZsdK2hDoABbUbVqQnaeRm4ExkiYp2Qr4AHDVQMptZmbtZ0AJStIYSZ+XdDfwVeB7wF6kGtD+wI/6OHwrYH1gRd36FUDDe0oRcSspIc0jNTGuBAR8tEn5uiQtlLRw5cqVuf9dZmZWPXl78R0q6UrgAeD9wNnAthExJSJujIjLgUOBPYeycJJ2IDXpnUyqfe1PSmYN53uOiO6I6IyIzlGj3PHQzKyd5X1Q9yLgu8DbIqLZxCjLgVl9nGMV8Dwwum79aODhJsccD9wWEV/PPv9e0hPAjZK+EBF/zlV6MzNrO3kT1JiIeLKvHSLiKWBmH9tXS+oBJgLfr9k0kdQ7sJFNSUmtVu9nP2RsNlJNm1ZKmAn7eKLCVsqboP4uaUxEPFK7UtKWwCMRsX7O85wOXCrpNlIHiCNJPQPPzc43FyAiDs/2nw/MkXQU8DNgDHAmsCgiHswZ08yGmxKmeweYNHV+KXGssbwJqtnDuRuROi/kEhGXZ0ntRFKyWQwcGBEPZLuMq9v/YkmbA58AvgH8FbgW+I+8Mc3MrD31maAkHZu9DeBISY/XbF4f2AO4eyABI2I2MLvJtr0brDuH1FHCzCzpyZreCh5NYtmSMQCM3W55oXGssf5qUMdkrwKmsub9oNXAUlIznZlZeTo702vBo5rPOXE64FHNW6XPBBUR2wFIug44NCIeK6VUZmY24uW6BxUR+xRdEDMzs1pNE5Sks4HjI+KJ7H1TEfHJIS+ZmZmNaH3VoHYCXlLzvplyprY0M7MRpWmCqm3WcxOfmZmVbZ1HY5D0ugFMs2FmZjYguTpJSDoVuCciLslm1F0AvAP4q6QDsnmazMzKsXBhKWGmndJwXGorSd6RJCYDh2XvDwB2BnbN1n8ZcBOgmZWnhOnewQ/otlreBDUa6B05/EDgvyLiNkmPAuX8lDEzsxEl7z2ovwCvzt7vB/wie78BzcfpMzMrRldXWgo2//xJzD9/UuFxrLG8CeoHwGWSrgG2II0sDqmp774iCmZm1tScOWkp2KLrOlh0XTnNiba2vE18x5Jm0x0HHBcRT2TrxwDfLqJgZmY2suUd6ug50nQX9evPGPISmZmZkb8GhaRNSU16W7Nm02BExA+HumBmZjay5X0O6p3Ad4EtG2wO0txQZmZmQyZvJ4mzgJ8Ar4yI9eoWJyczMxtyeZv4xgPviohlBZbFzCyfCRNKCTNmvL/yWilvgroZ+BfgjwWWxcwsn94p3wvWNau7lDjWWN4EdS5wmqSxwB3As7UbI2LRUBfMzMxGtrwJ6orstdHPCXeSMDOzIZc3QW1XaCnMzAZC2QhrUex8qTMnzwDgpHkzCo1jjeV9UPeBogtiZmZWK/eEhZIOkHSlpDslvSpbN1XSO4ornpmZjVS5EpSkycB/AfeSmvtekm1aHziumKKZmdlIlrcGdRwwLSI+AzxXs/5XpOGPzMzMhlTeBPXPwK0N1j8OvGzoimNmZpbkTVDLgNc3WL8nA3x4V9LRkpZIelpSj6Q9+tl/Q0lfyo55RtKDkj45kJhmZtZ+8nYz7wbOljQ1+/yqLLF8DZiRN5ikw0jj+h0N3JS9Xi1ph4h4sMlh3wNeCXSR7oGNBjbJG9PMhqHzzislzMFHzC8ljjWWt5v51yT9E3ANsDFwHfAMcFpEfGsA8Y4FLo6I3qkwj5G0P3AUcHz9zpL2A94BvDYiVmWrlw4gnpkNRyVM9w7QsW85QypZY7m7mUfECcBWwFuAXYFREfHFvMdL2hDoABbUbVoA7NbksEOA24FjJf1Z0r2Szpa0Wd64ZmbWnnJPWJjZFFgSEX9Zh1hbkbqlr6hbvwJ4Z5NjXgPsTqqtvQd4OXAOMBZ4b/3OkrpITYGMGzduHYpoZm2hOxt1reCaVM+1HYBrUq3Sbw1K0taSLpL0GCmZPCLpMUnnS9q6hPIF8KGI+HVE/Az4BPAeSaPrd46I7ojojIjOUaNGFVw0M2uZ6dPTUrArL5jElRdMKjyONdZnDUrSS0mdGbYA5gJ3AgLeCHwQ2F1SR0Q8kSPWKuB5UieHWqOBh5scsxx4KCL+WrPurux1HGvXxszMbJjorwZ1DGnUiB0j4lMRcV5EnBsRxwA7ARuRajT9iojVQA8wsW7TROCWJofdDIytu+fU293d4wOamQ1j/SWoScCpEbFWDScilgNfBt41gHinA1OyMfy2l3QW6X7SuQCS5kqaW7P/ZcBfgIskvVHS20nd1K+IiEcGENfMzNpMf50k3kBq4mvmJlKSyiUiLpe0JXAiMAZYDBxYM1r6uLr9H5f0TlLHiNuBx4AfAZ/PG9PMzNpTfwnqZcCjfWx/lAEOdRQRs4HZTbbt3WDdPcB+A4lhZmbtr78mvvWAf/SxPXKcw8zMbMD6q0EJuEHSc022D/Q5KjOzwSt4Jt1enkm3tfpLMDNLKYWZmVmdPhNURDhBmZlZS/j+kZm1n46OtBSs+4Quuk8oZ2BaW5vvIZlZ+1m0qJQwy5eOLSWONeYalJmZVZITlJmZVZITlJmZVVLue1CSNiBNVjgO2LB2W0TMbXiQmVkLSWp1EWwQciUoSW8A5gPbkR7efT479lnSZIJOUGZmNqTy1qDOJE2VsTNp7qadgX8Cvk0a+NXMrDzTpg1o95i35mdNznfchH08k24r5U1Qbwb2iognJP0D2CAiFkk6jjTS+P8prIRmZvV6p3wv2KSp80uJY43l7SQh4Mns/Upg2+z9n4HXDXWhzMzM8tagFgNvAu4HbgP+Q9LzwDTgvoLKZmbWWE/W9FbwaBLLlowBYOx2ywuNY43lTVCzgJdm708EfgJcB6wCDiugXGZmzXV2pteCRzWfc+J0wKOat0quBBURP6t5fz+wvaQtgMciShr33szMRpRc96AkXShp89p1EfEosKmkCwspmZmZjWh5O0l8FNikwfpNgMOHrjhmZmZJn018WTOesuUVdTPrrg8cBKwornhmZjZS9XcPahUQ2XJng+0BnDTUhTIzM+svQe1Dqj1dC7wHeLRm22rggYhYVlDZzMxsBOtvyvcbACRtB/wpIv5RSqnMzPqycGEpYaadct5a6/oagNadmodW3m7mDwBIGkvj0cx/OfRFMzNrooTp3sEP6LZa3tHMxwKXAXuS7jspe+21/tAXzcyseuoHnoX8g8/awOTtZn4maYqNHUhj8u0BvA+4C9i/mKKZmTXR1ZWWgs0/fxLzz59UeBxrLO9QR3sBB0XE3ZICWBkRN0t6BjgZuKawEpqZ1ZszJ70WPKr5outSU6JHNW+NvDWoTUhdziH15Ns6e38nA5xqQ9LRkpZIelpSj6Q9ch63u6TnJC0eSDwzM2tPeRPU3cAbsve/BY6U9Grg48BDeYNJOgw4CzgV2AW4Bbha0rh+jnsFadbeX+SNZWZm7S1vgjoL2CZ7/yVgP9LUG0cDXxhAvGOBiyNiTkTcFRHHAMuBo/o57gLgEuDWAcQyM7M2lreb+bya94skjSfVqB6MiFXNjqslaUOgAzitbtMCYLc+jjsaGA2cAnwxTywzM2t/eWtQa4iIJyNiUd7klNmK1B29fuy+FbxYO1uDpJ1IQyl9OCKe7y+ApC5JCyUtXLly5QCKZmZmVdNvDUrSJsBxpKGOXkN6/ul+4PvANyLiqSIKJmkj4HLgsxGxJM8xEdENdAN0dnb6kW6z4WrChFLCjBnvkdxaqb/RzDcgjcM3AfgpaSZdkZ6H+k/gAEl7RcRzzc/yglWkZ6lG160fDTzcYP8xwPbARZIuytatl4ql54ADI2JBjrhmNtz0TvlesK5ZxXZjt771V4PqAl4HTIiIP9RukLQjadr3acC3+wsUEasl9QATSbWvXhOBHzQ45CFgp7p1R2f7/xuwtL+YZmbWvvpLUO8FZtUnJ4CIWCzpy6QRJfpNUJnTgUsl3QbcDBwJjAXOBZA0Nzv34RHxLLDGM0+SHgGeiQg/C2VmNsz1l6DeCHy6j+0/Bz6fN1hEXC5pS+BEUhPeYlJT3QPZLn0+D2VmBkDviOIFjx4+c/IMAE6aN6PQONZYfwnqFUBf3eFWAi8fSMCImA3MbrJt736OnQHMGEg8MzNrT/11M18f6KsDxD/wSOZmZlaA/mpQAr6TDQrbyEZDXB4zMzOg/wR1SY5zzB2KgpiZmdXqb8r3j5VVEDMzs1rrNNSRmZlZ0fJOWGhmVh3nnVdKmIOP8ESFreQEZWbtp4Tp3gE69i1nSCVrzAnKzNqaeh/atWHH96DMrP10d6elYD3XdtBzbUfhcawx16DMrP1Mn55ea5r6XpxW9UWaPLgwV14wCXBTX6u4BmVmZpXkBGVmZpXkBGVmZpXkBGVmZpXkBGVmZpXkBGVmZpXkbuZm1n4Knkm3l2fSbS3XoMzMrJKcoMzMrJKcoMys/XR0pKVg3Sd00X1COQPT2tp8D8rM2s+iRaWEWb50bClxrDHXoMzMrJKcoMzMrJKcoMzMrJKcoMzMrJKcoMzMrJLci8/M2s+0aaWEmbCPJypspdITlKSjgc8BY4A/AJ+OiBub7HsocCSwC7AxcCcwKyJ+XFJxzayKSpjuHWDS1PmlxLHGSm3ik3QYcBZwKinp3AJcLWlck0P2Aq4FDsr2vwr4oaQ9SiiumZm1UNk1qGOBiyNiTvb5GEn7A0cBx9fvHBGfqls1U9JBwCFAw1qXmY0APVnTW8GjSSxbMgaAsdstLzSONVZagpK0IdABnFa3aQGw2wBOtTnw2FCVy8zag6QX3veOZa7Guw6ZOSdOBzyqeauU2cS3FbA+sKJu/QpgmzwnkPRx4JXApU22d0laKGnhypUrB1NWMzNrsbbpZi7pPcDXgQ9FxAON9omI7ojojIjOUaNGlVtAMytFzGv83oafMhPUKuB5YHTd+tHAw30dKOm9pFrT4RHhbjVmZiNAaQkqIlYDPcDEuk0TSb35GpL0flJymhIRVxRXQjMzq5Kye/GdDlwq6TbgZtIzTmOBcwEkzQWIiMOzzx8gJafPAr+U1HuvanVEPFpy2c3MrESlJqiIuFzSlsCJpAd1FwMH1txTqn8e6khSGc/Mll43AHsXW1ozM2ul0keSiIjZwOwm2/bu67OZGQCnlBNm2innlRPIGvJYfGZtrPbZoEYios/tbWu7csL4Ad3WcoIyG8b6S2CNDNukZm3HCcpsGKh/HkiTW1OO0pyfvU4tNsz88ycBHjS2VZygzIaxRg+y9iavKia13DW+67LXghPUouvSWH95E1Sz8rtWum7aZiQJMzMbWVyDMrM19FWLKasm0FfNr8qqWCttZ65BmZlZJbkGZWZraNfaiw0/rkGZmVkluQZlZu1nfDlhxoxfVk4ga8gJyszaz6xywnTN6i4nkDXkJj4zM6sk16DMKm5dhiuquuH4b7Kh5xqUmbWfydlSsJmTZzBz8oziA1lDrkGZtYnh2P3bD7ZaX5ygzCrCzV5ma3KCMiuRk5BZfk5QZhVT5WYvj9ZtZXKCMmuB4Xg/yZqrwgC87cgJysxyG0jtzs2ZNlhOUGY2aKUnoyPKCXPwEUMzk65rzOvGCcrMClXIl/O+gzw+p459e8oJZA05QZmtIzdhvcg1BCuCE5QZ+ZKNb2ZXyLXZa8E1qZ5rOwDXpFrFCcpskPqqPVS5y3hbuyB7LThBXXnBJMAJqlWcoMxquKnKyuZny5pzgjLLyfeczMpV+mjmko6WtETS05J6JO3Rz/57Zfs9Lel+SUeWVVYrh6R+l6E8n1mVxLw1l17+P1xyDUrSYcBZwNHATdnr1ZJ2iIgHG+y/HXAVcCHwYWB3YLaklRHxg/JKbvUG80eyLk0XQ/lHua7n8v0kq4qRMjJF2U18xwIXR8Sc7PMxkvYHjgKOb7D/kcCyiDgm+3yXpLcCnwWGdYIa6l5l63K+on6pNW1zH+L7P04o1s7W9e9hKP9uW53sVFYBJG0IPAl8MCK+X7P+W8COEbFXg2N+CdwRER+vWfc+4DJg04h4tlm8zs7OWLhw4WDLPKjjzawYvd9aRf+FzmDGGq+WDDZvSOqJiM7+9iuzBrUVsD6wom79CuCdTY7ZBvh5g/03yM63vHaDpC6gK/v4uKR76uKvGnixhzVfk7X5mqytctekrJ+OfSSmyl2TMjX58T6Qa/LqPDsNq158EdENdDfaJmlhnow9kviarM3XZG2+JmvzNVlbEdekzF58q4DngdF160cDDzc55uEm+z/HCP71YmY2EpSWoCJiNdADTKzbNBG4pclhtzbZf2Ff95/MzKz9lf0c1OnAFElTJW0v6SxgLHAugKS5kubW7H8usK2kM7P9pwJTgNPWIXbDpr8Rztdkbb4ma/M1WZuvydqG/JqU1ovvhYDS0cBxwBhgMfCZiPhltu16gIjYu2b/vYAzgDcCy4CvRsS5pRbazMxKV3qCMjMzy6P0oY7MzMzycIIyM7NKGnEJStI2ki6V9LCkJyX9TvIgOJLeIukaSY9L+rukWyRt1epytZqSqyWFpPe2ujytImkLSedIulvSU5L+JOnbkrZsddnKNtABr4czScdLul3S3yStlDRf0o5Ddf4Rl6CAucD2wLuBHbPPl0ras6WlaqFsfMMFwPXArkAHqaeku/LDvwP/aHUhKmAssC2pg9NOpMGb9wS+28pCla1mwOtTgV1Ij8hcLWlcSwvWOnsDs4HdSNNHPgf8XNIWQ3HyEddJQtLjwDERcVHNugeAcyJiXbqvtz1JtwDXRcQJrS5LlUh6M/D/SAl7BfC+iLiitaWqDkkHAlcCL4+Iv7W6PGWQ9Gvg9xExrWbdvcAVEdFowOsRRdJmwF+BQyJi/mDPNxJrUDcB75e0paT1JL0bGMXaY/6NCJK2Bt4GLJd0k6RHJN0o6R2tLlsrSdqcNChxV0Q80uryVNTLgGdIg0APe9mA1x2k1oZaC0g1CIPNSXnlsaE42UhMUO8nDYa8ivTHNY80wvpvW1qq1nlN9jqTNO/WvwI3Aj+T9KaWlar1zgV+GhFXt7ogVSTp5cDJwJyIeK7V5SlJXwNeb1N+cSrpLOC3pFGABm1YJChJp2Q3sfta9s52P4X0H+2dQCfwdWDucPsyHsA16f0/cF5EXBgRv4mILwC3k+bjGjbyXhNJHwHeBHyu1WUu2gD/dnqP2QyYDzxEuidlhqTTSZPKvicinh+Scw6He1BZb7P+epw9SBq94j5g54j4Xc3xPweWRsTU4kpZrgFck9HA/cBHIuI7NcdfAGwTEQcVV8pyDeCazAYOZ83OEetnn2+NiN2LKWH58l6TiHgy238z0izXAg6IiMcLLmJlaB3mtBspJJ0BfADYJyLuHqrzDovpNiJiFTlGN5e0afa2Prs/zzCpTfYawDVZShpC6l/qNr0euGPoS9Y6A7gmJ7D2eI93kGZy/v8FFK1l8l4TeOG+3NWk5LT/SEpOkAa8ltQ74PX3azZNZJjP8N0XpTFVD2OIkxMMkwQ1AHeTalCzJX0W+AtwCOk/2LtbWbBWiYiQ9HVgpqTfA78h3afbFfhESwvXIhHxEKn56gVKE7T9KSLub0mhWixLTgtIHSMOAV4q6aXZ5kez2QpGgtNJj6XcBtxMagZ/YcDrkSarPX6E9H/iMUm99+IeH4ofMCMqQUXEs1nX2K+Q2tA3IyWsjw1Fl8h2FRFnStoI+AawJfAHUvPN7/o+0kaQDtKPFoD/rtu2D+kZumEvIi7PHk4+kRcHvD4wIh5obcla5ujs9Rd162dC8+mI8xoW96DMzGz4GVb3XczMbPhwgjIzs0pygjIzs0pygjIzs0pygjIzs0pygjIzs0pygjJrU5Kul/TNms9LswfQB3PONSZmrP0saXz2uXMwMdaxXJ1Z7PFlx7bWcYKySpP0Y0n1DwH2bts++9Lar2792ZKelzStwTFT+hgUdeMmcVryxSxphqTFfexyKFD0HERjSA+1m5XOCcqq7gJgnya/nI8AHqBmLq9sRIzJpNFCmg3++yTpi3eNJSKeHrJSlyAiHo2Ivxcc4+GIeKbIGGbNOEFZ1f2ENN/Ox2pXSnoJaQywCyOidtTxQ4GlwCxgB0k7NjhnZF+8ayyDKaSk6ZLuk7Q6e51Wt/31km6Q9LSkeyQdKOlxSVMGEXONJr4G2z8s6W+S3pV9lqTjJP1R0lOS7pD04X5irNHkl3m1pGskPSnpTkkT647ZU9Kvs3/rCklnZCOB927fSNKZ2banJf1K0u5159hf0t3Z9htJgxfbCOMEZZWWTYZ3CTBFUu3/10mkaSIuqjtkKvCdbHqIH9C8FjVkJP0b8E3gTGBH0qRtsyVNyravB/wQeI40nt0U4CRgowLL9CngHODgiPhxtvoUUq3z48AOwJeB8yQNdEqVWcDZpDmzbge+l03DgaRtSSOe/wbYJYv3wSxWr6+RRr/+v9k+dwA/lTQmO8ergB8B1wA7Z/+Orw2wjDYcRIQXL5VegH8mzYK8X826nwBX1+23HbCaNI8VwL6kqSQ2qtlnSnaux+uWW/qIPz47prPJ9ptJNbnadRcDN2Xv/5WUnLat2b5bds4pfcSdASzuY/v1wDdrPi8lTQlyMqnWuUvNtpcCTwF71J3jTOCqms8BvLfR55rrML1m+7bZut2zz7OAe4H16q75M8CmWTlWA4fXbF8f+CNwSvb5VNKAtKrZ58QszvhW/3/0Ut4yokYzt/YUEfdKuoH0i3uBpLGkL/0P1O16BPCLeLG57nrS/aZDgMtr9nuS9Mu81mDus2wPXFi37ibgXdn7NwDLIk3j0et21pwQcah8CtgceHNE3FuzfgdgY1JNpXaE6JeQEttA/L7m/bLsdevsdXvgV7Fms+tNwIbA62pi3ty7MSKel3RrVsbac9SWc0imELf24gRl7eICYI6kLUi/yB+lZvJASetn68dKeq7muPVIzXy1CSoi4r6iC0z6xV+2m4D9Sc1qX6pZ39s8Ook0a3CtZwcY44X9IyKyubLy3C4I0mSHfW03e4HvQVm7uAJ4GvgwqSY1NyJqv1j3J81l1UmqHfUuBwPvKPhchLq6AAACD0lEQVT5mbuAt9et2x24M3t/Nylxjq3Z3kkxf389wH7AsZK+WLP+TlIt8dURcV/dMpRzGd0F7Fp3v3B3UrPeH7NlNTXXK/tx8TZevF53AW9Vlvkyu2IjjmtQ1hYi4ilJl5Huy7yCVKOqNZV0T2pR3frFku4hJbX/zNapZubPWisj4vk+ivH6utoZpOTzdeD7StOBLyAly8mkHoWQbvbfA1ySPUi7CWlm1ufov9awsaT65sgnI6J+0sAXRMTt2bNhCyRFRJwSEX+XdBpwWvbF/0vShJ27Av+IiO5+ypHXbODTpE4iZwGvIXX5/2akjitI+jbwVUmrgCXAZ4DR2bGQZqf9d+BMSbOBnUgz19pI0+qbYF685F2ACaQv9Jvr1o8mNTt9qMlxXwL+RKqxTMnO0Wh5XZPjx/dxzI7ZPkeSZmd+NnudVneO15OSwjOkZHUwqSZxWB//3hlNYi7Mtl9Pg04SNZ/fAvwPcGL2WcAxvFibWklKnhNrjsnTSaKzrpz1x+wJ/DqLsQI4gzU7qmxE6pyxItvnV2SdLGr2OSi7Tk+T7ldNxp0kRtziGXXNWkDSm4Dfkr7se1pdHrMqcoIyK0H2rNQTpC7Y40lNfCJ1BfcfoVkDvgdlVo7Nga8CrwIeIzXPfcbJyaw516DMzKyS3M3czMwqyQnKzMwqyQnKzMwqyQnKzMwqyQnKzMwq6X8B5y9zGOG6z+EAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "_5k_2\n", "mean log(likelihood) = -0.8669\n", "mode log(likelihood) = -0.1757\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XucHFWd9/HPFxAEwVUhhgSNQV1WEB8xM7qKXMWwCERZb6hRjJJMMIoXFlkRHkmU4B25aIAJCASDsuqja7hoUC5yU8jECxFhRZOgJIREWBW5BOLv+ePUQE9nursm6equnv6+X696dXfVqVO/9GsyvzmnTp2jiMDMzKxstmh3AGZmZsNxgjIzs1JygjIzs1JygjIzs1JygjIzs1JygjIzs1JygjIzs1JygjIzs1JygjIzs1Laqt0BFGWnnXaKiRMntjsMMyvawEB67elpetWrBlYBML5nfNPr7mYDAwPrImJMo3KjNkFNnDiRJUuWtDsMMytaf3967etretUD/Sn59fQ1P/l1M0kr85QbtQnKzLpEAYlpkBNTe/kelJmZlZITlJl1tv7+p7r5mmygf+DJbj5rPXfxmVlnmzkzvRbQ1Xf5zMsBd/W1i1tQZmZWSk5QZmZWSk5QZmZWSk5QZmZWSk5QZmYVJD251dtnxXOCMjOzUvIwczPrbBHFVLsQYDYApwCaWshlrI6WtqAk7SfpB5LulRSSpjUof4Ck/5a0WtLDkn4t6f0tCtfMzNqo1V182wPLgI8Aj+QovzdwO/BWYE/gHKBf0rsKi9DMzEqhpV18EXElcCWApItylD+tatc5kg4E3gJc2vQAzazzDC6zMdD8KYn6T0qzU/TNLWYqJauvE+9BPRP4U7uDMLOSWLq0sKpXr/A6UO3UUQlK0uHAQcBraxzvA/oAJkyY0MLIzMys2TpmmLmk15K69T4cEbcOVyYi+iOiNyJ6x4xpuFijmZmVWEckKEn7AFcBn4qIc9odj5mZFa/0CUrSfqTkNDsizmh3PGZm1hotvQclaXvgxdnHLYAJkvYCHoiIeyR9FnhVRByUlT8AuAKYB1wqaefs3A0RsbaVsZuZWWu1epBEL3Btxec52XYxMA0YB7yo4vg0YDvg+GwbtBKYWFyYZtYxZsworOpJB3o13XZSFDRNSLv19vbGkiVL2h2GmXWYwQlh01RHFfuzqY5G6+/MVpI0EBG9jcqV/h6UmZl1JycoM+tsAwOFzCIBsGr5OFYtH1dI3dZYRz2oa2a2kd6sp6iArrf5J88E4JSFs5tetzXmFpSZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSE5SZmZWSh5mbWWcrcMaYGaeeV1jd1pgTlJl1tsEl3wswftfVhdVtjbmLz8zMSskJysw6W19f2gqw6PwpLDp/SiF1W2NOUGbW2ebPT1sBll7bw9Jri+tCtPqcoMzMrJScoMzMrJScoMzMrJScoMzMrJScoMzMrJT8oK6ZdbZJkwqretzEVYXVbY05QZlZZytouXeAvrn9hdVtjbmLz8zMSskJyszMSskJysw6m5S2AsyZOps5U2cXUrc11tIEJWk/ST+QdK+kkDQtxzkvk3S9pEey8z4lFfTTaGZmpdHqFtT2wDLgI8AjjQpLeiZwNbAGeGV23seB4wqM0czMSqClo/gi4krgSgBJF+U4ZSqwHfDeiHgEWCbpJcBxkk6PiCgsWDMza6uy34N6DXBDlpwG/QgYD0ysLiypT9ISSUvWrl3bohDNzKwIZU9QO5O69yqtqTg2RET0R0RvRPSOGTOm8ODMzKw4ZU9QZmbWpco+k8R9wNiqfWMrjplZtzvvvMKqPvzoRYXVbY2VPUHdAnxe0tMj4tFs32RgFbCibVGZWXkUtNw7QM/riptGyRpr9XNQ20vaS9Je2bUnZJ8nZMc/K+knFadcCjwMXCRpT0lvBj4BeASfmdko1+p7UL3AL7JtW2BO9v7T2fFxwIsGC0fEX0gtpvHAEuBrwJeB01sXspmVWn9/2gowcE0PA9f0FFK3Ndbq56CuA2rOAhER04bZdzuwX3FRmVlHmzkzvRbQ1Xf5BVMAd/W1i0fxmZlZKTlBmZlZKTlBmZlZKTlBmZlZKTlBmZlZKTlBmZlZKZV9Jgkzs/oKfGb/lIWzC6vbGnOCMrOu40W5O0OuLj5JXrvCzMxaKm8L6l5JPwAuAH7oefDMrDR6sqmIBkY+20Ms3Hifpj71vv+kNDtF39xiplKy+vImqMOA9wHfBf6cLdd+UUT8vqjAzMxyWbq0sKpXrxhfWN3WWK4uvoi4OiLeRZq09XPAG4D/kXSNpKmSnl5kkGZm1n1GNMw8Iv43Ir4WEb3Ah4G9gUuAVZI+J2n7IoI0M7PuM6JRfJLGAe8FpgHPA75Fui81HjiRtJzG65sbopmZdaNcCSpbKPD9wMHAMuAsYGG2XtNgmduAO4sI0szMuk/eFtSFwDeB10REraEyq4G5TYnKzMy6Xt4ENS4iHq5XICIeIa2Qa2bWOjNmFFb1pAO9UGE75U1Qf5M0LiLur9wpaUfg/ojYsvmhmZnlUNBy7wBTpi8qrG5rLO8ovlrzgmwDrG9SLGZmZk+q24KSdFz2NoBjJD1UcXhLYF88MMLM2mlwBonBGSWaaNXycQCM33V10+u2xhp18R2bvQqYDmyoOLYeWAEc0/ywzMxy6u1NrwXMwDb/5JmAZzVvl7oJKiJ2BZB0LfDmiHiwJVGZmVnXyzVIIiIOLDoQMzOzSjUTlKSzgBMj4u/Z+5oi4sN5LyhpFvBxYBzwG+CjEXFDnfLvAk4AdgP+CvwYOD4i7st7TTMz6zz1WlAvA55W8b6W3B2/ko4EzgRmATdmr1dJ2iMi7hmm/GtJc/0dD3wfGAvMAxYCB+W9rpmZdZ6aCaqyW6+JXXzHkZbpmJ99PlbSIcAHSHP5VXsN8KeI+Er2ebmks4GzmxSPmZmV1IhmM68k6cUjWWZD0tZAD7C46tBi0qzow7kJGCdpipKdgHcAV25KzGZm1jnyThZ7GnBXRFwsSaSkchDwF0lviIif5ahmJ9KzU2uq9q+hxgzoEXGLpHeQuvS2zeK9mjSj+nBx9gF9ABMmTMgRkpl1vCVLCqt6xqnnFVa3NZa3BTUVuCt7/wZgL+DVwALgswXEBYCkPUjdeZ8htb4OAXYGhv2piYj+iOiNiN4xY8YUFZaZlUlPTyEP6UJ6QNcP6bZP3rn4xgJ/yt4fCvxXRNwq6QEg758v60gP+o4dpu5aI/JOBG6NiC9mn38t6e/ADZI+GRF/qnGemZl1uLwtqD8DL8jeHwz8JHu/FbXn6RsiItYDA8DkqkOTgZtrnLYdQ2evoOLzJt8/M7NRpK8vbQVYdP4UFp0/pZC6rbG8v+S/C1wq6WrgOcCPsv17AXeP4HqnA9MkTZe0u6QzSavxngsgaYGkBRXlFwFvkvQBSS/Mhp2fBSwdbli6mXWh+fPTVoCl1/aw9Npiug+tsbxdfMcBK4EJwAkR8fds/zjgnLwXi4jLsiU6Ts7OXQYcGhErsyITqspfJGkH4EPAl4G/ANcA/5n3mmZm1pnyTnX0BClBVO//yjDFG9U1j/Sw7XDHDhhmn597MjPrQnlbUEjajtSl91yGdg1GRHyv2YGZmVl3y/sc1OuBbwI7DnM4SM83mZmZNU3eQRJnAlcAz4uILao2JyczM2u6vF18E4E3RsSqAmMxMxu5SZMKq3rcRP/Ka6e8Ceom4F+A3xcYi5nZyA0u+V6Avrn9hdVtjeVNUOcCX5I0HrgdeLzyYEQsbXZgZmbW3fImqO9kr8P9OeFBEmZm1nR5E9SuhUZhZraplM22FrnXTs1tztTZAJyycHbT67bG8j6ou7JxKTMzs+bJPeGqpDdIulzSHZKen+2bLslLr5uZWdPlSlCSpgL/BfyO1N33tOzQlsAJxYRmZmbdLG8L6gRgRkR8DHiiYv/PSNMfmZmZNVXeBPXPwC3D7H8IeGbzwjEzM0vyJqhVwG7D7N8PP7xrZmYFyDvMvB84S9L07PPzJe0LfAGYXURgZma5nHdeYVUffvSiwuq2xvIOM/+CpH8CrgaeDlwLPAZ8KSK+VmB8Zmb1FbTcO0DP64qbRskay70eVEScJGkusAepa/COiHiosMjMzKyr5U5Qme2A5RHx5yKCMTMbsf5sBrYCWlID1/QAbkm1S8NBEpKeK+lCSQ8Ca4D7JT0o6XxJzy0+RDOzOmbOTFsBLr9gCpdfMKWQuq2xui0oSc8AbgSeAywA7gAEvBR4J7CPpJ6I+HvRgZqZWXdp1MV3LGnWiD0j4r7KA5JOA24GPgR8vpjwzMysWzXq4psCnFadnAAiYjXwWeCNRQRmZmbdrVGCegmpi6+WG7MyZmZmTdUoQT0TeKDO8QfwVEdmZlaARglqC+AfdY5HjjqGkDRL0nJJj0oayGakqFd+a0mfzs55TNI9kj48kmuamVnnaTRIQsD1kp6ocXxEz1FJOhI4E5hF6h6cBVwlaY+IuKfGad8Cngf0kZb7GAtsO5LrmtkoVsBKuoO8km57NUowc5p8veOAiyJifvb5WEmHAB8ATqwuLOlg4CDgRRGxLtu9oskxmZlZCdVNUBHRtAQlaWugB/hS1aHFwN41TjsCuA04TtJRwCPAVcAnh5tmSVIfqaXFhAkTmhS5mZm1w4juH22mnUgr8K6p2r8G2LnGOS8E9gFeDryF9MzVIcBFwxWOiP6I6I2I3jFjxjQjZjMru56etBWg/6Q++k8qbjJaq2+kc/G12hakgRjvioi/AEj6EPAjSWMjojrZmVm3Wbq0sKpXrxhfWN3WWCtbUOuADaRBDpXGAhs9CJxZDdw7mJwyv81e3YdnZjaKtSxBRcR6YACYXHVoMmnKpOHcBIyXtH3FvsGVfVc2N0IzMyuTVragAE4HpkmaLml3SWcC44FzASQtkLSgovylwJ+BCyW9VNJrScPUvxMR97c4djMza6Hc96AkbQW8itS1tnXlsYhYMOxJVSLiMkk7AicD44BlwKERMdgamlBV/iFJrwfOJo3mexD4PvCJvHGbmVlnypWgJL0EWATsSnp4d0N27uOkpd9zJSiAiJgHzKtx7IBh9t0FHJy3fjMzGx3ytqDOIN0/2os0oGEv4J+Ac0itITOz9pgxo7CqJx3olXTbKW+CeiWwf0T8XdI/gK0iYqmkE0jdb/+nsAjNzOoZXPK9AFOmLyqsbmss7yAJAQ9n79cCu2Tv/wS8uNlBmZmZ5W1BLSPN5vAH4FbgPyVtAGYAdxcUm5lZYwNZN1wBs0msWj4OgPG7rm563dZY3gQ1F3hG9v5k4ArgWtLDt0cWEJeZWT69vem1gFnN5588E/Cs5u2SK0FFxI8q3v8B2F3Sc4AHIwqc697MzLpWrntQkr4uaYfKfRHxALCdpK8XEpmZmXW1vIMk3svwiwRuCxzVvHDMzMySul18WTeesu3ZVSvrbgkcxsbLZ5iZmW22Rveg1pGWuwjgjmGOB3BKs4MyMzNrlKAOJLWeriEtGPhAxbH1wMqIWFVQbGZm1sUaLfl+PYCkXYE/RsQ/WhKVmVleS5YUVvWMU8/baJ+kmuU9qLm58g4zXwkgaTzDz2b+0+aHZmaWQ0HLvYMf0G23vLOZjyetzbQf6b6TstdBWzY/NDOzzVOvtbOpYuEw15na9MsY+YeZn0FaYmMP0px8+wJvIy2/fkgxoZmZ5dDXl7YCLDp/CovOn1JI3dZY3qmO9gcOi4g7JQWwNiJukvQY8Bng6sIiNDOrZ/789FpnVvPqVk/eFs/Sa1P3oWc1b4+8LahtSUPOIY3ke272/g681IaZmRUgb4K6E3hJ9v6XwDGSXgB8ELi3iMDMzKy75e3iOxPYOXv/aeCHwDtJy72/t4C4zMysy+UdZr6w4v1SSRNJLap7ImJdrfPMzMw2Vd4W1BAR8TCwtMmxmJmZPalhgpK0LXACaaqjF5Kef/oD8G3gyxHxSKERmpnVM2lSYVWPm+iZ3Nqp0WzmW5Hm4ZtEuu90Bekh3T2ATwFvkLR/RDxRuxYzswINLvlegL65tYeuW/EajeLrA14MTIqIN0XEiRHxiYh4Iylp7QbMGMkFJc2StFzSo5IGJO2b87x9JD0hadlIrmdmZp2pUYJ6KzA3In5TfSAilgGfJc0okYukI0kjAk8DXgHcDFwlaUKD854NLAB+kvdaZmbW2RolqJeSuvhq+TGw5wiudxxwUUTMj4jfRsSxwGrgAw3OuwC4GLhlBNcys24gpa0Ac6bOZs7U2YXUbY01SlDPBtbWOb4WeFaeC0naGugBFlcdWgzsXee8WcBY4NQ81zEzs9GhUYLaEqg3AOIf5J/JfKesbPUS8Wt46iHgISS9jLRi77sjYkOjC0jqk7RE0pK1a+vlVTMzK7tGw8wFfCObFHY42zQ5nqcuLG0DXAYcHxHL85wTEf1AP0Bvb69XDjMz62CNEtTFOepYkPNa60hLdoyt2j8WuG+Y8uOA3YELJV2Y7dsCkKQngEMjorq70MzMRolGS76/r1kXioj1kgaAyaSHfAdNBr47zCn3Ai+r2jcrK//vwIpmxWZmZuWzSVMdbYbTgUsk3QrcBBwDjAfOBZC0ACAijoqIx4EhzzxJuh94LBvibmZmo1hLE1REXCZpR+BkUhfeMlJX3cqsSN3noczMNnLeeYVVffjRXqiwnVrdgiIi5gHzahw7oMG5s4HZTQ/KzDpXQcu9A/S8rrhplKyxvAsWmpmZtZQTlJl1tv7+tBVg4JoeBq7pKaRua6zlXXxmZk01c2Z6LaCr7/ILpgDu6msXt6DMzKyUnKDMzKyUnKDMzKyUnKDMzKyUnKDMzKyUnKDMzKyUPMzczDpbFLeyzikLZxdWtzXmFpSZmZWSE5SZmZWSE5SZdbaenrQVoP+kPvpPKm4yWqvP96DMrLMtXVpY1atXjC+sbmvMLSgzMyslt6DMbFSQ1O4QrMncgjIzs1JyC8rMRoVYuPE+TW19HNY8bkGZmVkpuQVlZh1tcC3dIgaDTzrQCxW2kxOUmXW0bD3dQhLUlOmLCqjV8nIXn5mZlZJbUGbW0SYVWPeq5eMAGL/r6gKvYrU4QZlZRyvyLtH8k1MHomc1b4+Wd/FJmiVpuaRHJQ1I2rdO2TdLWixpraS/Sfq5pDe2Ml4zM2uPliYoSUcCZwKnAa8AbgaukjShxin7A9cAh2XlrwS+Vy+pmZnZ6NDqLr7jgIsiYn72+VhJhwAfAE6sLhwRH6naNUfSYcARwA2FRmpmZm3VshaUpK2BHmBx1aHFwN4jqGoH4MEa1+iTtETSkrVr125aoGZmVgqt7OLbCdgSWFO1fw2wc54KJH0QeB5wyXDHI6I/InojonfMmDGbE6uZlYykYTcbvTpmFJ+ktwBfBI6MiJXtjsfMzIrVygS1DtgAjK3aPxa4r96Jkt4KLACOigg/2m3Wxaonhe3JJoQtYrj5jFPPK6BWy6tlCSoi1ksaACYD3644NBn4bq3zJL0duBh4b0R8p9gozazTFLeerh/QbbdWd/GdDlwi6VbgJuAYYDxwLoCkBQARcVT2+R2k+03HAz+VNHivan1EPNDi2M3MrIVamqAi4jJJOwInA+OAZcChFfeUqp+HOoYU4xnZNuh64IBiozWzTlBkJ9yi86cAnjS2XVo+SCIi5gHzahw7oN5nM7NqRcxiPmjptT2AE1S7dMwoPjOzsqs17D0iWhzJ6ODlNszMrJTcgjKzhvI8EOtWwsZD4DW1PXGMFk5QZgbkS0JmreQEZWa5VbcQwK0EK44TlJkN0WlJaHAGiZ4C6h43cVUBtVpeTlBmHczdctCbvRZxB6xvbn8BtVpeHsVnZmal5BaU2ShQa/RYve46jzizsnOCMrPS2JQuyyIHt8+ZOhuAUxbOLvAqVosTlFnJ+T6TdSsnKDMrnU4bSWjFcIIy6xBl/6Xteeis2ZygzEpitHbl1ft3OXlZPU5QZtYUmzIqcLQmZWsOJyizkhltw7/L3jXZCm5FbhonKDNrm2Yk48EFC4uY8+Hwo71QYTs5QZkVwF1XrTM/ey0iQfW8bqBxoRzcitw0nurIzMxKyS0oswJtyhRENjIzCqx74Jo0R3qzWlI2Mk5QZtbRipxv/PILpgDFJigPoKjNCcqsAd9PMmsPJygznISsfdzdW1vLE5SkWcDHgXHAb4CPRsQNdcrvD5wOvBRYBXwhIs5tRaxWXnkSynDdI5uTiHw/yay1WpqgJB0JnAnMAm7MXq+StEdE3DNM+V2BK4GvA+8G9gHmSVobEd9tXeRWrdEv+k1JDs1OKJu0dIOTjZVIt89v2OoW1HHARREx+OjCsZIOAT4AnDhM+WOAVRFxbPb5t5L+FTge6NoEtVmtgBH8YLc6ObQyoYy22RqsuzS7S7r698Km9lA0m1qViSVtDTwMvDMivl2x/2vAnhGx/zDn/BS4PSI+WLHvbcClwHYR8Xit6/X29saSJUs2N+bNOt/Mijf4G6yI/62zmT3k1Z6yOblD0kBE9DYq18oW1E7AlsCaqv1rgNfXOGdn4MfDlN8qq2915QFJfTw188lDku5qEM+6xmF3DX8fG/N3MlQpv48i/4xskJhK+X20yjB/wI/k+3hBnkKjahRfRPST87EISUvyZPBu4e9jY/5OhvL3MZS/j6GK+D5aOdXROmADMLZq/1jgvhrn3Fej/BN08V8uZmbdoGUJKiLWAwPA5KpDk4Gba5x2S43yS+rdfzIzs87X6sliTwemSZouaXdJZwLjgXMBJC2QtKCi/LnALpLOyMpPB6YBX2pCLEXOkNKJ/H1szN/JUP4+hvL3MVTTv4+WjeJ78oLpQd0TSA/qLgM+FhE/zY5dBxARB1SU3x/4Ck89qPt5P6hrZjb6tTxBmZmZ5eH1oMzMrJScoMzMrJS6PkFJ2lnSJZLuk/SwpF9JnvhG0qskXS3pIUl/k3SzpJ3aHVc7KblKUkh6a7vjaQdJz5F0tqQ7JT0i6Y+SzpG0Y7tjayVJsyQtl/SopAFJ+7Y7pnaQdKKk2yT9VdJaSYsk7dms+rs+QQELgN2BNwF7Zp8vkbRfW6Nqo2y+w8XAdcCrgR7SyMluH9r/H8A/2h1Em40HdiENdHoZaRLn/YBvtjOoVqqY9Po04BWkx2SukjShrYG1xwHAPGBv4HWkZ1R/LOk5zai86wdJSHoIODYiLqzYtxI4OyKaMZy940i6Gbg2Ik5qdyxlIemVwP8jJes1wNsi4jvtjaocJB0KXA48KyL+2u54iibp58CvI2JGxb7fAd+JiOEmve4akrYH/gIcERGLNrc+t6DSsh9vl7SjpC0kvQkYw8ZzAHYFSc8FXgOslnSjpPsl3SDpoHbH1i6SdiBNUNwXEfe3O54SeibwGGky6FEtm/S6h9TDUGkxqRXR7XYg5ZUHm1GZExS8nTQh8jrSf7KFpBnXf9nWqNrnhdnrHNI6XP8G3AD8SNLL2xZVe50L/DAirmp3IGUj6VnAZ4D5EfFEu+NpgXqTXu/c+nBK50zgl6RZgDbbqExQkk7NbmTX2w7Iip9K+qF7PdALfBFYMNp+GY/gOxn8mTgvIr4eEb+IiE8Ct5HW5xoV8n4fkt4DvJy0CvSoNcL/M4PnbA8sAu4l3ZOyLibpdNKism+JiA1NqXM03oPKRps1GnF2D2k2i7uBvSLiVxXn/xhYERHTi4uytUbwnYwF/gC8JyK+UXH+BcDOEXFYcVG2zgi+j3nAUQwdHLFl9vmWiNinmAhbK+/3EREPZ+W3J612LeANEfFQwSGWgjZhXbtuIOkrwDuAAyPizmbVO6qW2xgUEevIMdu5pO2yt9XZfgOjrHU5gu9kBWlKqX+pOrQbcHvzI2uPEXwfJ7Hx3I+3k1Z1/u8CQmuLvN8HPHlP7ipScjqkW5ITpEmvJQ1Oev3tikOT6dJVvpXmVD2SJicnGKUJagTuJLWg5kk6HvgzcATph+1N7QysXSIiJH0RmCPp18AvSPfpXg18qK3BtUFE3EvqwnqS0kJtf4yIP7QlqDbKktNi0sCII4BnSHpGdviBbNWC0e500qMotwI3kbq+n5z0uptkLcf3kH4WHpQ0eB/uoWb84dLVCSoiHs+GyH6O1Je+PSlhva8ZQyQ7VUScIWkb4MvAjsBvSN04v6p/pnWBHtIfKwD/U3XsQNKzc6NaRFyWPZh8Mk9Nen1oRKxsb2RtMSt7/UnV/jlQfzniPEblPSgzM+t8o+o+i5mZjR5OUGZmVkpOUGZmVkpOUGZmVkpOUGZmVkpOUGZmVkpOUGYdStJ1kr5a8XlF9sD55tQ5ZDHGys+SJmafezfnGpsYV2927Ymtvra1jxOUlZqkH0iqfghw8Nju2S+tg6v2nyVpg6QZw5wzrc5kqE+vcZ22/GKWNFvSsjpF3gwUvf7QONJD7GYt5wRlZXcBcGCNv5yPBlZSsXZXNgPGVNLsILUm+32Y9It3yBYRjzYt6haIiAci4m8FX+O+iHisyGuY1eIEZWV3BWmtnfdV7pT0NNIcYF+PiMqZxt8MrADmAntI2nOYOiP7xTtk25wgJc2UdLek9dnrjKrju0m6XtKjku6SdKikhyRN24xrDuniG+b4uyX9VdIbs8+SdIKk30t6RNLtkt7d4BpDuvwyL5B0taSHJd0haXLVOftJ+nn2b10j6SvZLOCDx7eRdEZ27FFJP5O0T1Udh0i6Mzt+A2myYusyTlBWatkieBcD0yRV/rxOIS0PcWHVKdOBb2TLQnyX2q2oppH078BXgTOAPUmLts2TNCU7vgXwPeAJ0jx204BTgG0KjOkjwNnA4RHxg2z3qaRW5weBPYDPAudJGukSKnOBs0jrZN0GfCtbfgNJu5BmOv8F8Irseu/MrjXoC6TZr9+flbkd+KGkcVkdzwe+D1wN7JX9O74wwhhtNIgIb95KvQH/TFr1+OCKfVcAV1WV2xVYT1q3CuB1pCUktqkoMy2r66Gq7eY615+YndNb4/hNpJZc5b6LgBuz9/9GSk67VBzfO6tzWp3rzgaW1Tl+HfDVis8rSMuAfIbU6nxFxbFnAI8A+1bVcQZwZcXnAN463OeK72FmxfFdsn37ZJ/nAr8Dtqj6zh8DtsviWA8cVXF8S+D3wKnZ59NIE9GtdIwcAAADJklEQVSqoszJ2XUmtvvn0Vvrtq6ezdw6Q0T8TtL1pL+4F0saT/ql/46qokcDP4mnuuuuI91vOgK4rKLcw6S/zCttzn2W3YGvV+27EXhj9v4lwKpIS3cMuo2hiyA2y0eAHYBXRsTvKvbvATyd1FKpnCH6aaTENhK/rni/Knt9bva6O/CzGNrteiOwNfDiimveNHgwIjZIuiWLsbKOyjibsoS4dRYnKOsUFwDzJT2H9Bf5A1QsGChpy2z/eElPVJy3BambrzJBRUTcXXTApL/4W+1G4BBSt9qnK/YPdo9OIa0UXOnxEV7jyfIREdn6WHluFwRpkcN6x82e5HtQ1im+AzwKvJvUkloQEZW/WA8hrV3VS2odDW6HAwcV/PzMb4HXVu3bB7gje38nKXGOrzjeSzH//waAg4HjJP3fiv13kFqJL4iIu6u2Zq5j9Fvg1VX3C/chdev9PtvWU/F9ZX9cvIanvq/fAv+qLPNlXo11HbegrCNExCOSLiXdl3k2qUVVaTrpntTSqv3LJN1FSmqfyvapYuXPSmsjYkOdMHarap1BSj5fBL6ttBT4YlKynEoaUQjpZv9dwMXZg7TbklZlfYLGrYanS6rujnw4IqoXC3xSRNyWPRu2WFJExKkR8TdJXwK+lP3i/ylpgc5XA/+IiP4GceQ1D/goaZDImcALSUP+vxpp4AqSzgE+L2kdsBz4GDA2OxfSyrT/AZwhaR7wMtKqtdZt2n0TzJu3vBswifQL/aaq/WNJ3U7vqnHep4E/klos07I6htteXOP8iXXO2TMrcwxpNebHs9cZVXXsRkoKj5GS1eGklsSRdf69s2tcc0l2/DqGGSRR8flVwP8CJ2efBRzLU62ptaTkObninDyDJHqr4qw+Zz/g59k11gBfYehAlW1IgzPWZGV+RjbIoqLMYdn39CjpftVUPEii6zavqGvWBpJeDvyS9Mt+oN3xmJWRE5RZC2TPSv2dNAR7IqmLT6Sh4P5PaDYM34Mya40dgM8DzwceJHXPfczJyaw2t6DMzKyUPMzczMxKyQnKzMxKyQnKzMxKyQnKzMxKyQnKzMxK6f8DjmG7KK7HYjQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "_5k_3\n", "mean log(likelihood) = -0.8558\n", "mode log(likelihood) = -0.1705\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X2cXVV97/HPFxQE0aqAIaPGYK0VxCtkRqvIoxpEIL68aks1FqMkE0TRK7W0Fq4EJVAVEXzgYQKK0dBS8bY1CDVUHpQHhZnYKynIFQ1oSQhBqBUhROB3/1h74OQwZ2bP5Kxz9pz9fb9e+3XO2XvtvX4cJvObtfbaaykiMDMzq5ptuh2AmZnZWJygzMyskpygzMyskpygzMyskpygzMyskpygzMyskpygzMyskpygzMyskpygzMyskp7W7QBy2WWXXWL27NndDsPMumFkJL3292erYt3IOgD6+vuy1dGrRkZG7ouIXScq17MJavbs2QwPD3c7DDPrhqGh9Do4mK2KkaGUBPsH8yXBXiXprjLlejZBmVmNZUxMo5yY8vM9KDMzqyQnKDPrPUNDT3bzZTIyNPJEN5/l4S4+M+s9ixen14xdfZctvgxwV19ObkGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmZUgaYttrP3WXk5QZmZWSR5mbma9JyLfpVeMvlsCwMmA5merrtbcgjIzs0pygjIzs0rqaIKSdICkb0u6W1JIWjBB+YMk/Yuk9ZIekvQTSe/vULhmNl3192ddagNg6MRBhk7MPyltnXX6HtROwBpgebFNZF/gFuAzwHrgzcCQpE0RcXG2KM1selu9OnsV6+/0OlC5dTRBRcTlwOUAki4qUf60pl3nSjoYeAfgBGVm1sOm4z2oZwMPjHVA0qCkYUnDGzdu7HBYZmbWTtMqQUk6AngjMOY8+hExFBEDETGw664TriZsZmYVNm0SlKTXk7r1PhwRN3U7HjMzy2taJChJ+wFXAJ+IiHO7HY+ZmeVX+ZkkJB0AfAc4OSLO6nY8ZjYNLFqUvYo5B3s13dw6mqAk7QS8tPi4DTBL0t7A/RHxS0mnA6+JiDcW5Q8iJadzgIsl7Vac+1hEeBSEmY0t83LvAPMWrsxeR911uotvAPhxse0AnFK8/2RxfCbwhw3lFwA7Ah8jPQc1ut3cmXDNzKxbOv0c1DVAyznpI2LBGJ8XjFXWzKylkaL7LeNsEuvWzgSgb/f12eqou8rfgzIzm7SBgfSacVbzZSctBuDkFUuy1VF302IUn5mZ1Y8TlJmZVZITlJmZVZITlJmZVZITlJmZVZITlJmZVZKHmZtZ7xkezl7FolPPz15H3TlBmVnvybzcO/gB3U5wF5+ZmVWSE5SZ9Z7BwbRltPKCeay8YF7WOurOCcrMes+yZWnLaPXV/ay+On9XYp05QZmZWSU5QZmZWSU5QZmZWSU5QZmZWSU5QZmZWSX5QV0z6z1z5mSvYubsddnrqDsnKDPrPaNLvmc0uHQoex115y4+MzOrJCcoMzOrJCcoM+s9UtoyOmX+Ek6ZvyRrHXXnBGVmZpXU0QQl6QBJ35Z0t6SQtKDEOa+UdK2kh4vzPiFl/tPIzMy6rtMtqJ2ANcBHgIcnKizp2cCVwAbg1cV5fwUcnzFGMzOrgI4OM4+Iy4HLASRdVOKU+cCOwHsj4mFgjaSXA8dLOjMiIluwZmbWVVW/B/U64AdFchr1XaAPmN2ViMzMrCOqnqB2I3XvNdrQcGwLkgYlDUsa3rhxY/bgzMwsn56aSSIihoAhgIGBAXf/mdXV+ednr+KIo1dmr6Puqp6g7gFmNO2b0XDMzOypMi/3DtD/hvzTKdVd1bv4bgT2l/SMhn1zgXXAnV2JyMzMOqLTz0HtJGlvSXsXdc8qPs8qjp8u6XsNp1wMPARcJGkvSW8H/gbwCD4za21oKG0ZjVzVz8hV/VnrqLtOt6AGgB8X2w7AKcX7TxbHZwJ/OFo4In5DajH1AcPAl4HPAWd2LmQzm3YWL05bRpddOI/LLpyXtY666/RzUNcALWeBiIgFY+y7BTggX1RmZlZFVb8HZWZmNeUEZWZmleQEZWZmleQEZWZmleQEZWZmlVT1mSTMzCavA49JnrxiSfY66s4tKDMzq6RSCUrSrrkDMTMza1S2BXW3pEslvcXLrZtZ5fX3p20KJI25NRs6cZChE/NPSltnZe9BHQ68D/gW8OtiNdyLIuLnuQIzM5uy1auzV7H+zr7sddRdqQQVEVcCV0p6DmkZ9vcBfyvpWuBC4FsRsSlfmGZmnRUrtvys+d2Jo84mNUgiIv4rIr4cEQPAh4F9ga8D6yT9naSdcgRpZmb1M6lh5pJmAu8FFgAvBP6B1ILqAz5Omq38Te0N0czM6qhUgirWYXo/cAiwBvgCsKJYDmO0zM3AT3MEaWZm9VO2BfVV4O+B10VEq3WO1wNL2xKVmZnVXtkENTMiHhqvQEQ8TFqA0MysuxYtyl7FnINb/a1u7VI2Qf1W0syIuLdxp6SdgXsjYtv2h2ZmNkWZl3sHmLdwZfY66q7sKL5WD+duD2xuUyxmZmZPGLcFJen44m0Ax0h6sOHwtsD+eGCEmVXNSNH9NsXZJMpYt3YmAH27r89WR91N1MV3XPEqYCHwWMOxzcCdwDHtD8vMbCsMDKTXjLOaLztpMeBZzXMaN0FFxO4Akq4G3h4RD3QkKjMzq72yUx0dnDsQMzOzRi0TlKQvAB+PiN8V71uKiA+3PTIzM6u18UbxvRJ4esP7Vttek6lQ0rGS1kraJGlE0v4TlH+3pH+X9JCkeyR9Q9Juk6nTzMymn5YtqMZuvXZ18Uk6EjgbOBa4rni9QtKeEfHLMcq/njQZ7ceAfwZmAOcAK4A3tiMmMzOrpikv+S7ppZKeMcnTjietI7UsIm6LiONIUyR9oEX51wH/GRGfj4i1EfFD4IvAn0w1bjMzmx7KThZ7GnB7RHytWFF3FakF8xtJbykSx0TX2A7oB85oOrSKtGzHWK4HTpM0D7gM2Bn4c+DyMnGbWU0ND2evYtGp52evo+7KtqDmA7cX798C7A28FlgOnF7yGruQHu7d0LR/AzDmPaWIuJGUkFaQnrvaSHom671jlZc0KGlY0vDGjRtLhmVmPWcrlnwvq2/39X5IN7OyCWoG8J/F+8OAf4yIm0jdbfvkCAxA0p5FHZ8itb4OJSWzMf90iYihiBiIiIFdd901V1hmZtYBZRPUr4EXF+8PAb5XvH8arefpa3YfaSaKGU37ZwD3tDjn48BNEfHZiPhJRHyXNLDiLyS9sGS9ZlY3g4Npy2jlBfNYecG8rHXUXdkE9S3gYklXAs8Dvlvs3xu4o8wFImIzMALMbTo0F7ihxWk7suX0SjR8nvIADzPrccuWpS2j1Vf3s/rqvN2IdVd2uY3jgbuAWcAJEfG7Yv9M4NxJ1Hcm8HVJN5EGQBxDWi7+PABJywEi4qii/EpgmaQPkJLiTOAsYPVYw9LNzKx3lJ3q6FHgc2Ps//xkKouIS4o1pE4iJZs1wGERcVdRZFZT+YskPQv4UFH/b4CrgL+eTL1mZjb9lG1BIWlHUpfe89myey0i4p/KXiciziE9bDvWsYPG2PdF0kAJMzOrkbLPQb0J+HvSc0jNgjR83MzMrG3KDjQ4G/gO8MKI2KZpc3IyM7O2K9vFNxt4a0SsyxiLmVl7zJmTvYqZs/3rMLeyCep64I+Bn2eMxcysPUaXfM9ocOlQ9jrqrmyCOg84Q1IfcAvw+8aDEbG63YGZmVm9lU1QlxavY/3J4EESZmbWdmUT1O5ZozAzaycVM7BFZKvilPlLADh5xZJsddRd2Qd175q4lJmZWfuUns9O0lskXSbpVkkvKvYtlOSVbc3MrO1KJShJ84F/BH5G6u57enFoW+CEPKGZmVmdlW1BnQAsioiPAo827P8hafojMzOztiqboP4IuHGM/Q8Cz25fOGZmZknZBLUOeNkY+w/AD++amVkGZYeZDwFfkLSw+PwiSfsDnwGW5AjMzGzKzj8/exVHHL0yex11V3aY+Wck/QFwJfAM4GrgEeCMiPhyxvjMzCYv83LvAP1vyD+dUt2VXg8qIk6UtBTYk9Q1eGtEPJgtMjMzq7XSCaqwI7A2In6dIxgzs7YYKmZly9iSGrmqH3BLKqcJB0lIer6kr0p6ANgA3CvpAUkXSHp+/hDNzCZp8eK0ZXTZhfO47MJ5Weuou3FbUJKeCVwHPA9YDtwKCHgF8C5gP0n9EfG73IGamVm9TNTFdxxp1oi9IuKexgOSTgNuAD4EfDpPeGZmVlcTdfHNA05rTk4AEbEeOB14a47AzMys3iZKUC8ndfG1cl1RxszMrK0mSlDPBu4f5/j9eKojMzPLYKIEtQ3w+DjHo8Q1tiDpWElrJW2SNFLMSDFe+e0kfbI45xFJv5T04cnUaWZm089EgyQEXCvp0RbHJ/UclaQjgbOBY0ndg8cCV0jaMyJ+2eK0fwBeCAySlvuYAewwmXrNrGYyrqQ7yivp5jdRgjmlzfUdD1wUEcuKz8dJOhT4APDx5sKSDgHeCPxhRNxX7L6zzTGZmVkFjZugIqJtCUrSdkA/cEbToVXAvi1OextwM3C8pKOAh4ErgL/1NEtmZr1tslMdbY1dSCvwbmjavwF4U4tzXgLsR5qY9h3Ac4AvAn3AO5sLSxokdQUya9astgRtZtNQf5qGiJF80xANnZimURpcOvTEPkljlo0OdDn2ok4mqKnYhjQQ490R8RsASR8CvitpRkRskewiYoi0NAgDAwP+iTCrq9Wrs1ex/s6+7HXUXScT1H3AY6RBDo1mAE95ELiwHrh7NDkVbiteZ/HU1piZWdfEii0/a3534ugVkxoivjUiYjMwAsxtOjSXNGXSWK4H+iTt1LBvdGXfu9oboZmZVUnHElThTGCBpIWS9pB0Nul+0nkAkpZLWt5Q/mLg18BXJb1C0utJw9QvjYh7Oxy7mZl1UOkuPklPA15D6lrbrvFYRCwf86QmEXGJpJ2Bk4CZwBrgsIgYbQ3Nair/oKQ3kQZG3Aw8APwz8Ddl4zYzs+mpVIKS9HJgJbA76eHdx4pzf08aYVcqQQFExDnAOS2OHTTGvtuBQ8pe38zMekPZFtRZpPtHe5MGNOwN/AFwLqk1ZGZWHYsWZa9izsFeSTe3sgnq1cCBEfE7SY8DT4uI1ZJOIHW//Y9sEZqZTdbQ0MRlttK8hSuz11F3ZQdJCHioeL8ReEHx/j+Bl7Y7KDMzs7ItqDXAq4BfADcBfy3pMWARcEem2MzMpmZ0BonRGSUyWLd2JgB9u6/PVkfdlU1QS4FnFu9PAr4DXE16+PbIDHGZmU3dwEB6zTjF0LKTFgOe1TynUgkqIr7b8P4XwB6Sngc8EJ5kysymmVZz5lm1lLoHJekrkp7VuC8i7gd2lPSVLJGZmVmtle3iey/p4djfNu3fATgKeH87gzIz64TmufPA8+dVybgJqujGU7E9t2ll3W2Bw/GErWZmlsFELaj7SMtdBHDrGMcDOLndQZmZmU2UoA4mtZ6uIi0YeH/Dsc3AXRGxLlNsZmZWYxMt+X4tgKTdgV9FxOMdicrMbGsMD2evYtGp52evo+7KDjO/C0BSH2PPZv799odmZjZFGR/QHeUHdPMrO5t5H2ltpgNI951UvI7atv2hmZlZnZWdi+8s0hIbe5Lm5Nsf+FPS8uuH5gnNzGyKBgfTltHKC+ax8oJ5Weuou7LPQR0IHB4RP5UUwMaIuF7SI8CngCuzRWhmNlnLlqXXjLOar746dSN6VvN8yragdiANOYc0ku/5xftb8VIbZmaWQdkE9VPg5cX7fweOkfRi4IPA3TkCMzOzeivbxXc2sFvx/pPAvwLvIi33/t4McZmZWc2VHWa+ouH9akmzSS2qX0bEfa3OMzMzm6qyLagtRMRDwOo2x2JmZvaECROUpB2AE0hTHb2E9PzTL4BvAp+LiIezRmhmNllz5mSvYuZsz/KW20SzmT+NNA/fHNJ9p++QHtLdE/gE8BZJB0bEo62vYmbWYaNLvmc0uDTfEHZLJmpBDQIvBeZExH80HpC0F2nZ90XAuXnCMzOzuppomPk7gaXNyQkgItYAp5NmlChN0rGS1kraJGlE0v4lz9tP0qOS1kymPjMzm54mSlCvIHXxtfJvwF5lK5N0JGnI+mnAPsANwBWSZk1w3nOB5cD3ytZlZjUmpS2jU+Yv4ZT5S7LWUXcTJajnAhvHOb4ReM4k6jseuCgilkXEbRFxHLAe+MAE510IfA24cRJ1mZnZNDZRgtoWGG8AxOOUnMlc0nZAP7Cq6dAqYN9xzjsWmAGcWqKOQUnDkoY3bhwvr5qZWdVNNEhCwDeKSWHHsv0k6tqFlMw2NO3fALxpzMqlV5KWlH9tRDymCZrsETEEDAEMDAzEuIXNzKzSJkpQXytxjeXtCKSZpO2BS4CPRcTaHHWYmVl1TbTk+/vaWNd9pDWlZjTtnwHcM0b5mcAewFclfbXYtw0gSY8Ch0VEc3ehmZn1iLKzmW+1iNgMjABzmw7NJY3ma3Y38Epg74btPOCO4v1Y55iZWY+Y0lx8W+FM4OuSbgKuB44B+kiJB0nLASLiqIj4PbDFM0+S7gUeKZ7BMjMb2/nnZ6/iiKO9UGFuHU1QEXGJpJ2Bk0hdeGtIXXV3FUXGfR7KzKyUzMu9A/S/If90SnXX6RYUEXEOcE6LYwdNcO4SYEnbgzIzs8rp2D0oM7OOGRpKW0YjV/UzclV/1jrqruMtKDOz7BYvTq8Zu/ouu3Ae4K6+nNyCMjOzSnKCMjOzSnKCMjOzSnKCMjOzSnKCMjOzSnKCMjOzSvIwczPrPZF/tZ2TVyzJXkfdOUGZWc+aaA05qzZ38ZmZWSW5BWVmvad/yymIYsWWhzV/66sYOjHNUjG4NO+USnXmBGVmvWf16uxVrL+zL3sddecuPjMzqyQnKDMzqyQnKDMzqyQnKDMzqyQnKDMzqySP4jOz3rNoUXpdtixbFXMO9kKFuTlBmVnvGV3uPWOCmrdwZbZrW+IuPjMzqyS3oMys94zk735bt3YmAH27r89eV105QZlZ7xkYyF7FspMWA57VPKeOd/FJOlbSWkmbJI1I2n+csm+XtErSRkm/lfQjSW/tZLxmZtYdHU1Qko4EzgZOA/YBbgCukDSrxSkHAlcBhxflLwf+abykZmZmvaHTXXzHAxdFxOjQmuMkHQp8APh4c+GI+EjTrlMkHQ68DfhB1kjNzKyrOtaCkrQd0A+sajq0Cth3Epd6FvBAu+IyM7Nq6mQX3y7AtsCGpv0bgN3KXEDSB4EXAl9vcXxQ0rCk4Y0bN25NrGY2TUh6yma9Ydo8ByXpHcBngXdHxF1jlYmIoYgYiIiBXXfdtbMBmpm1MFYSdTKdWCfvQd0HPAbMaNo/A7hnvBMlvRNYDhwVEX5828yeYotVc9cWryflq2/Rqefnu7gBHUxQEbFZ0ggwF/hmw6G5wLdanSfpz4CvAe+NiEvzRmlmPWH3/FVM5gHd5iXnoT3Lzve6To/iOxP4uqSbgOuBY4A+4DwAScsBIuKo4vOfk+43fQz4vqTRe1WbI+L+DsduZmYd1NEEFRGXSNqZ1PCeCawBDmu4p9T8PNQxpBjPKrZR1wIH5Y3WzKatC/JXsfKCeYAnjc2p41MdRcQ5wDktjh003mczs1Kuzl/F6qv7ASeonDwXn1kPKjM6LCI6EInZ1DlBmdm04CHZ9eMEZdbDPHrMpjMnKLNprI6tiuak64Tbu6bNTBJmZlYvbkGZ9QC3KprMLl7vzFfFzNnr8l3cACcos9pq1T041ui+aTcqcGnxmjFRDy4dyndxA9zFZ2ZmFeUWlFlNTaVbcLxRgZNpkZmV4QRlVnF1HKm31TpwD+6U+UsAOHnFkvyV1ZQTlJm1hQdqWLs5QZlNE5166LabLTa3Fq2RB0mYmVkluQVlVhFVaT1UYXqkKsTQKeP9f6/7ABO3oMzMrJLcgjLroFIPvHqwQa3UqbU4WU5QZhlUpbuuto4uXi/MV8URR3uhwtycoMy6oO5/NWdP4G8oXjMmqP43jOS7uAFOUGZZubvOrUmbOicoM+uabAn8qjZdZxwjV/UDbknl5ARlNkVuGZTTle7MjF17oy67cB6QN0HVfX5DDzM3M7NKcgvKbCvVfcCD5VP3e5gdb0FJOlbSWkmbJI1I2n+C8gcW5TZJ+oWkYzoVq3WGpAm3dpxT9ryym1m31OXnsqMJStKRwNnAacA+wA3AFZJmtSi/O3B5UW4f4HTgi5Le0ZmI66Gdv4hz/WKfyjV6/R+v2Vh6KXF1uovveOCiiFhWfD5O0qHAB4CPj1H+GGBdRBxXfL5N0p8AHwO+lT3aaSjHD2K7r9nyxu8Uusqm2r3WqutkvOvVvbvFqmMqP/ft/nfciYEa6tRoEEnbAQ8B74qIbzbs/zKwV0QcOMY53wduiYgPNuz7U+BiYMeI+H2r+gYGBmJ4eHhrY96q882sO0Z/q+X8F7yEJVu81s3W5A5JIxExMFG5TragdgG2BTY07d8AvKnFObsB/zZG+acV11vfeEDSIDBYfHxQ0u1bE3CP2QW4r9tBVJy/o3Iq/z114k/LCRJT5b+jrbWVf8C/uEyhnhrFFxFDwFC346giScNl/mKpM39H5fh7mpi/o/bo5CCJ+4DHgBlN+2cA97Q4554W5R+lx/86MTOru44lqIjYDIwAc5sOzSWN0hvLjS3KD493/8nMzKa/Tj8HdSawQNJCSXtIOhvoA84DkLRc0vKG8ucBL5B0VlF+IbAAOKPDcfcCd31OzN9ROf6eJubvqA06NorviQqlY4ETgJnAGuCjEfH94tg1ABFxUEP5A4HPA68A1gGfjojzOhq0mZl1XMcTlJmZWRmeLNbMzCrJCcrMzCrJCaqGJL1G0pWSHpT0W0k3SNql23FVjZIrJIWkd3Y7niqR9DxJX5T0U0kPS/qVpHMl7dzt2LpJk5wM28bnBFUzxVyGq4BrgNcC/aRRkR62/1R/CTze7SAqqg94AWnA0yuB9wAHAH/fzaC6abKTYdvEPEiiZiTdAFwdESd2O5Yqk/Rq4P+QEvgG4E8j4tLuRlVtkg4DLgOeExH/3e14Ok3Sj4CfRMSihn0/Ay6NiLEmw7YJuAVVI5KeD7wOWC/pOkn3SvqBpDd2O7YqkfQs0oTEgxFxb7fjmUaeDTxCmhS6VorJsPtJvRONVgH7dj6i3uAEVS8vKV5PAb4CvBn4AfBdSa/qWlTVcx7wrxFxRbcDmS4kPQf4FLAsIh7tdjxdMN5k2Lt1Ppze4ATVAySdWtzIH287iCf/f58fEV+JiB9HxN8CN5PW3upZZb8jSX8BvAr4q27H3A2T+FlqPGcnYCVwN+melFlb9NRs5jV2FvCNCcr8kicn3r216ditQK/fyC37HS0A9iQt19J47BJJN0bEfnnCq4yy3xPwRHK6vPh4RERsyhVYxU1lMmybgBNUD4iI+ygxu7ukO0nTRf1x06GXAbe0P7LqmMR3dCJPnevxFtIqzv+SIbRKKfs9wRP36q4gLb90aEQ8mDO2KouIzZJGJ8P+ZsOhuXj17ylzgqqRiAhJnwVOkfQT4MfAn5GGm3+oq8FVRETcTeqqekLRkvpVRPyiK0FVUJGcVpEGRrwNeKakZxaH7y9WL6ibM4GvS7oJuJ7Ubf7EZNg2eU5QNRMRZ0naHvgcsDPwH8BbIuL/djcym2b6SX/YAPy/pmMHk56zq5WIuKR4UPkknpwM+7CIuKu7kU1ffg7KzMwqyaP4zMyskpygzMyskpygzMyskpygzMyskpygzMyskpygzMyskpygzKYpSddI+lLD5zslfWwrr7nF4oyNnyXNLj4PbE0dU4xroKh7dqfrtu5xgrJKk/RtSd9rcWyP4pfWIU37vyDpMUmLxjhnwTiToD6jRT1d+cUsaYmkNeMUeTuQe52hmaSJYM06zgnKqu5C4OAWfzkfDdwF/NvojmKWjPnA3wELW1zzIdIv3i226TbRaUTcHxG/zVzHPRHxSM46zFpxgrKq+w5pTZ33Ne6U9HTgL4CvRETjsuxvB+4ElgJ7StprjGtG8Yt3i21rgpS0WNIdkjYXr4uajr9M0rWSNkm6XdJhkh6UtGAr6tyii2+M4++R9N+S3lp8lqQTJP1c0sOSbpH0ngnq2KLLr/BiSVdKekjSrZLmNp1zgKQfFf+tGyR9XmlBv9Hj20s6qzi2SdIPJe3XdI1DJf20OP4D0oTGVjNOUFZpxeJ3XwMWSGr8eZ1HWiTuq02nLAS+EREPkWaRbtWKahtJ/xP4Emmpir2As4FzJM0rjm8D/BPwKGn+ugXAycD2GWP6CPBF0hIY3y52n0pqdX6QtKTI6cD5kg6f5OWXAl8grZt1M/APxbIbSHoBaYbzHwP7FPW9q6hr1GeAI4H3F2VuAf5V0sziGi8C/hm4Eti7+O/4zCRjtF4QEd68VXoD/ggI4JCGfd8BrmgqtzuwGdit+PwG0tIR2zeUWVBc68Gm7YZx6p9dnDPQ4vj1pJZc476LgOuK928mJacXNBzft7jmgnHqXQKsGef4NcCXGj7fSVoW5FOkVuc+DceeCTwM7N90jbOAyxs+B/DOsT43fA+LG46/oNi3X/F5KfAzYJum7/wRYMcijs3AUQ3HtwV+DpxafD6NNAGtGsqcVNQzu9s/j946t3k2c6u8iPiZpGtJf3GvktRH+qX/501Fjwa+F092111Dut/0NuCShnIPkf4yb7Q191n2AL7StO864K3F+5cD6yIt5THqZuBx2u8jwLOAV0fEzxr27wk8g9RSaZwh+umkxDYZP2l4v654fX7xugfww9iy2/U6YDvgpQ11Xj96MCIek3RjEWPjNRrjvHGSMVoPcIKy6eJCYJmk55H+Ir+fhgUEJW1b7O+T9GjDeduQuvkaE1RExB25Ayb9xd9p1wGHkrrVPtmwf7R7dB4NK+IWfj/JOp4oHxFRrJdV5nZBkBY3HO+42RN8D8qmi0uBTcB7SC2p5RHR+Iv1UNL6VgOk1tHodgTwxszPz9wGvL5p337ArcX7n5ISZ1/D8QHy/PsbAQ4Bjpf0vxv230pqJb44Iu5o2tq5XtFtwGub7hfuR+r5x4q0AAAB2ElEQVTW+3mxbabh+yr+uHgdT35ftwF/oiLzFV6L1Y5bUDYtRMTDki4m3Zd5LqlF1Wgh6Z7U6qb9ayTdTkpqnyj2SdJuY1SzMSIeGyeMlzW1ziAln88C31Ra8nsVKVnOJ40ohHSz/3bga8WDtDuQVl99lIlbDc+Q1Nwd+VBENC8S+ISIuLl4NmyVpIiIUyPit5LOAM4ofvF/H9iJ9Iv/8YgYmiCOss4B/hdpkMjZwEtIQ/6/FGngCpLOBT4t6T5gLfBRYEZxLqQVaP8SOEvSOcArSavTWt10+yaYN29lN2AO6Rf69U37Z5C6nd7d4rxPAr8itVgWFNcYa3tpi/Nnj3POXkWZY4A7ijjuABY1XeNlpKTwCClZHUFqSRw5zn/vkhZ1DhfHr2GMQRINn18D/BdwUvFZwHE82ZraSEqecxvOKTNIYqApzuZzDgB+VNSxAfg8Ww5U2Z40OGNDUeaHFIMsGsocXnxPm0j3q+bjQRK127yirlkXSHoV8O+kX/Yj3Y7HrIqcoMw6oHhW6nekIdizSV18Ig0F9z9CszH4HpRZZzwL+DTwIuABUvfcR52czFpzC8rMzCrJw8zNzKySnKDMzKySnKDMzKySnKDMzKySnKDMzKyS/j9k8RyC72/oTgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Evaluate ELBO distribution on GFP dataset, decoder epsilon = 1e-6\n", "\n", "n_z_samples = 128\n", "\n", "for it in range(3) :\n", "\n", " TRAIN_SIZE = 5000\n", " train_size_str = \"%ik\" % (TRAIN_SIZE/1000)\n", " num_models = [1, 5, 20][it]\n", " RANDOM_STATE = it + 1\n", "\n", " X_train, y_train, gt_train = get_experimental_X_y(random_state=RANDOM_STATE, train_size=TRAIN_SIZE)\n", "\n", " L = X_train.shape[1]\n", "\n", " vae_suffix = '_%s_%i' % (train_size_str, RANDOM_STATE)\n", " \n", " print(vae_suffix)\n", "\n", " vae_0 = build_vae(latent_dim=20, n_tokens=20, seq_length=L, enc1_units=50)\n", "\n", " vae_0.encoder_.load_weights(\"models/vae_0_encoder_weights%s.h5\" % vae_suffix)\n", " vae_0.decoder_.load_weights(\"models/vae_0_decoder_weights%s.h5\"% vae_suffix)\n", " vae_0.vae_.load_weights(\"models/vae_0_vae_weights%s.h5\"% vae_suffix)\n", "\n", " #Compute multi-sample ELBO on test set\n", " log_mean_p_vae_test, mean_log_p_vae_test, log_p_vae_test = evaluate_elbo(vae_0.encoder_, vae_0.decoder_, X_train, n_samples=n_z_samples)\n", "\n", " #Log Likelihood Plot\n", " plot_min_val = None\n", " plot_max_val = None\n", "\n", " f = plt.figure(figsize=(6, 4))\n", "\n", " log_p_vae_test_hist, log_p_vae_test_edges = np.histogram(log_mean_p_vae_test, bins=50, density=True)\n", " bin_width_test = log_p_vae_test_edges[1] - log_p_vae_test_edges[0]\n", "\n", " mean_log_p_vae_test = np.mean(log_mean_p_vae_test)\n", " mode_log_p_vae_test = log_p_vae_test_edges[np.argmax(log_p_vae_test_hist)] + bin_width_test / 2.\n", "\n", " print(\"mean log(likelihood) = \" + str(round(mean_log_p_vae_test, 4)))\n", " print(\"mode log(likelihood) = \" + str(round(mode_log_p_vae_test, 4)))\n", "\n", "\n", " plt.bar(log_p_vae_test_edges[1:] - bin_width_test/2., log_p_vae_test_hist, width=bin_width_test, linewidth=2, edgecolor='black', color='orange')\n", "\n", " plt.xticks(fontsize=14)\n", " plt.yticks(fontsize=14)\n", "\n", " if plot_min_val is not None and plot_max_val is not None :\n", " plt.xlim(plot_min_val, plot_max_val)\n", "\n", " plt.xlabel(\"VAE Log Likelihood\", fontsize=14)\n", " plt.ylabel(\"Data Density\", fontsize=14)\n", "\n", " plt.axvline(x=mean_log_p_vae_test, linewidth=2, color='red', linestyle=\"--\")\n", " plt.axvline(x=mode_log_p_vae_test, linewidth=2, color='purple', linestyle=\"--\")\n", "\n", " plt.tight_layout()\n", " plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:tensorflow]", "language": "python", "name": "conda-env-tensorflow-py" }, "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 }