ch14densityratio.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np               \n",
    "import matplotlib.pyplot as plt \n",
    "import statsmodels.api as sm\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 密度比の推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# データ設定 \n",
    "n, m = 100, 200\n",
    "nu_mean, nu_sd = -0.5, 1\n",
    "a_mean,  a_sd = 1, 0.8\n",
    "newdat = np.linspace(-4,4,500).reshape(500,1)  # 予測点\n",
    "tnu =  norm.pdf(newdat, nu_mean, nu_sd)            # 確率密度の計算\n",
    "tde = (norm.pdf(newdat, a_mean, a_sd)+tnu)/2\n",
    "tw = tnu/tde                                   # 予測点上での真の密度比"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# データ生成\n",
    "nu = np.random.normal(loc=nu_mean, scale=nu_sd, size=n).reshape(n,1)\n",
    "ma = np.random.binomial(m,0.5); mb = m-ma\n",
    "de = np.r_[np.random.normal(loc=nu_mean,scale=nu_sd,size=ma).reshape(ma,1), np.random.normal(loc=a_mean, scale=a_sd, size=mb).reshape(mb,1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzsnXd4VFXawH9nJj2k0EvoXQRBiCBI\nEESKZcWGFBULwmLZVVRUdAXrWtfPsqwKCsrasKCirBRpIr0TQEDASEISWkggPZk53x9vJpmESTKT\nTPr5Pc997sy9557zTjLz3nPf8xaltcZgMBgMdQdLVQtgMBgMhsrFKH6DwWCoYxjFbzAYDHUMo/gN\nBoOhjmEUv8FgMNQxjOI3GAyGOoZR/AaDwVDHMIrfYDAY6hhG8RsMBkMdw6eqBXBFo0aNdNu2bata\nDIPBYKgxbNu27ZTWurE7baul4m/bti1bt26tajEMBoOhxqCU+tPdtsbUYzAYDHUMo/gNBoOhjmEU\nv8FgMNQxjOI3GAyGOoZR/AaDwVDHKFXxK6VaKaVWKaV+U0rtVUo96KKNUkq9rZQ6pJTarZTq7XTu\nDqXU73nbHd7+AAaDwWDwDHfcOXOBR7TW25VSIcA2pdRyrfU+pzZXAZ3ytn7Au0A/pVQDYCYQCei8\naxdprc949VMYDAaDwW1KVfxa6wQgIe/1OaXUb0AE4Kz4RwHztdRx3KiUCldKNQcGA8u11kkASqnl\nwEjgc69+CuDQiVSeXBgtb1ShHSr/vSr8vpjj5F+niunn/PMF16pi2ronw/myFx7Dx2LBx6LwsSp8\nrRasjtcWee1rVfhY89pYFH4+VoL8rAT6yT7Iz0qgrw/B/o5jPgT7WfPHMRgMtR+PAriUUm2Bi4FN\nRU5FALFO7+PyjhV33FXfk4HJAK1bt/ZELADSsnLZHJPk8XUG8LEowoP8qB/kS/1g2TcI9qNJSAAR\n4YG0CA+kRXgALcIDCfC1VrW4BoOhnLit+JVS9YBvgIe01meLnnZxiS7h+PkHtZ4NzAaIjIz0uAJ8\nhyb1WDD50vzOHTXkNfkvCg1e9Lw+77wuLGxx7XXx1+giF58/hnsyOLBrTY5NY7Nrcm32/Nc5dju5\nNk1u3nHZa3LtdrJy7KTn2MjIziU920Z6to2MbBvpOblkZNtIy7KRkWPjVGoWp1KzKI1G9fzp1KQe\nnZvWo2PTELo0DeHCFqEE+1fLIHCDweACt36tSilfROl/qrVe6KJJHNDK6X1LID7v+OAix1eXRdDS\nqOfvQ7/2DSui61pPVq6N5PQcktKyOZOeTXJ6DqfTsjmekkl8cgbHkjOIT8kgITkz/wax4cjp/Ost\nCro0C6V363Ai29ZnYMfGNA7xr8JPZDAYSkIVnVWe10CMvx8DSVrrh4ppcw3wAHA1srj7tta6b97i\n7jbA4eWzHejjsPkXR2RkpDa5eqofNrsmPjmDQydS+f3EOQ4eT2V/4ll+SziHzV74e9QjIozLOzdm\nWLemXNQyzKwhGAwVjFJqm9Y60q22bij+gcBaIBqw5x1+EmgNoLV+L+/m8G9k4TYduEtrvTXv+rvz\n2gO8qLWeV5pQRvHXLDKybUQfS2H70TNsOHyajUdOk5Vrzz/ftmEQ1/WK4KbeEbRpGFyFkhoMtRev\nKv6qwCj+YtAa1q2D3r0hKEiObd8OublwySXnuyVVEZk5NjYeOc2q/SdYHJ2Yv3agFAzp0oQ7B7Ql\nqlMj8xRgMHgRTxS/idytKezeDf37Q1SUKHsHs2dDv37Qty+sWVN18jkR4GtlcJcmPDuqOxunX8En\nE/txY+8IfK0WVu4/wYS5mxn55lqW7Ek8bwHbYDBUPEbx1wQ+/RQiI2HTJmjUCLKzC8516iTHtm6F\nwYPhsccgJ6fKRC2Kj9XCwE6NeOOWXmx44gqmjehCs9AADhw/x5RPtjFq1jrW/n6yqsU0GOoUxtRT\n3fnwQ7jnHnk9eTK8/jqEhBRuk54Or74KL7wANps8FSxaBOHhlS+vG2Tl2liwJZZ3Vh7i5DkxA13d\noxlPX9uN5mGBVSydwVAzqbs2/vR0eOghGDpUtkaNvC9cZbJsGVx9tSjzV1+FadNKbr9uHYwZA8eO\nQa9e8OuvEFx9F1Mzsm3MW/8H/155iPRsG0F+Vp665gLG921t7P8Gg4fUXcW/dCmMHFnw/uKLYdgw\nuPJKGDgQAmvQbDInR8w4f/4J06fDP//p3nWxsXLTu+kmuaYGKND45Aye+2EfS/YmAjC0axNevuki\nEwtgMHhA3VX8f/4JX3wBP/8Ma9dCllMkqr+/KMXGbtUirh4cPAjvvw+vvQYWD5ZjUlIgNLRGKH1n\nftgVz1PfRnM2M5fGIf68d1sf+rSpX9ViGQw1grqr+J3JyBBTx88/w/LlYgbav7/g/M03y4z62mvh\n0kvBWktz0MTHw6xZYv+vATeChJQMHvxiJ5v/SMLXqnhuVHfG9fU8d5PBUNcwit8VmZkQECCvY2Kg\nXbuCc40aiS395pth+HB5Oqgqzp2DxYvFVl9eRW23i8//rl1i9pk+3TsyVjA5NjsvLv6Nj9bHADDl\n8g48PrKLsfsbDCVg/Phd4VD6ABER8iTw4IPQvj2cOgXz58N110GTJrB5c9XJ+c47MG6cyFZeLJaC\nmf5TT8H//lf+PisBX6uFZ667kFdvvggfi+K9NYd5/Jvd5NrspV9sMBhKpe4ofmd8fWUB9M034dAh\n2LsXXnxRPGFyc6F794K2778vrpFZpWeuLDdnz4q7JsCoUd7p89pr4bnnJOp3/Hg4fNg7/VYCt0S2\nYs4dkQT4Wvhyaxx/+3yHUf4GgxeoO6YedzlxQmb9IMq+SRNRyOHhYn65/XYYMKBi7OXPPw8zZsCg\nQbB6tffGsNvFjPXtt9Cnj7h9VqU5y0O2Hz3DnXM3czYzl7/0bMGbY3phtRizj8HgjDH1lAeH0gdx\nqZw+HXr2hORkmf0PHAgdO8LMmZCQ4L1xU1PhjTfk9XPPeffGYrHA3LnQti1s2wZPPOG9viuB3q3r\n8/Hdfann78MPu+J57Ovd2O3Vb8JiMNQUjOIviXr1REnu3AnR0ZIOISICjhwR5XzqVEHb8j45zZ8v\nN5cBA+Dyy8vXlyvCw8XV1d8fwsLKL28lc3Hr+sy76xICfa18sz2O15cdqGqRDIYaizH1eIrNBqtW\nyfbiiwXHhw2D1q0lrULfvp7N2LWGbt3E3XTBArjlFu/L7eDYMbl51VDWHDzJ3R9twWbXvHJTD8Zc\nYlw9DQYw7pyVz5Ej0KFDwfsePWDSJLjtNqjvRgBSVpYs6i5aJLEHvr4VJ6szJ06IfJU1npf4bNNR\nnvw2Gh+L4uO7+3JZxxqemsNg8AJetfErpeYqpU4opfYUc36aUmpn3rZHKWXLq7yFUipGKRWdd64G\naXIPad8eDhyQXDqNGolZ6O9/hxYtYMIEmWWXhL+/uFtu2lR5Snj1arlBPfVU5YznRcb3a81fL29P\nrl3zwGfbOZacUdUiGQw1Cnds/B8hlbVcorV+TWvdS2vdC5gOrClSWnFI3nm37kQ1ls6dJZHasWPw\n5ZeSHygzUzxpnLNpVpcnLB8fOH1a0kH8+GNVS+Mxj43oyqDOjTmTnsP9n24nO9e4eRoM7lKq4tda\n/wKUWCPXiXHA5+WSqKbj5wejR0uaiEOH4L//lbw5IDeCiy4Sj6D4eDn2ySfwyiuQmFi5cg4cWJD4\nbcIEyXNUg7BaFG+O6UVEeCA7Y5N5cfG+qhbJYKgxeM2rRykVhDwZfON0WAPLlFLblFKTvTVWjaFD\nB7j++oL3ixfDnj3iEdSmDYwdC88+K55DGzZUvnyPPgrXXANnzkiMgnOBlxpAg2A/Zt3aG1+r4uMN\nf/LzvuNVLZLBUCPwpjvnX4B1Rcw8l2mtewNXAfcrpQYVd7FSarJSaqtSauvJk7W0ItONN4pt/aab\nJKhqwQJ5KvDxES+gyjYDWSzw8cfijbRpEzz+eOWO7wV6tQrn8ZFdAXhi4e78+r4Gg6F4vKn4x1LE\nzKO1js/bnwC+BfoWd7HWerbWOlJrHdm4JqVO9gSlxEf/66/hjz8kihYkTcQ//lE1MjVsKDcgHx94\n+21ZpK5h3H1ZO/q3b8ip1GymL4w2dXwNhlLwiuJXSoUBlwPfOx0LVkqFOF4DwwGXnkF1kpYtC+z6\nf/tb4WjdmBhZdE1J8d54Nptk/kxMlPWFkyclFQVIWur//EcqfnXp4r0xKwmLRfH6LT0JCfBh+b7j\nfL0trqpFMhiqNe64c34ObAC6KKXilFITlVJTlFJTnJrdACzTWqc5HWsK/KqU2gVsBhZrrZd4U/ga\nzcaN4gHUujW89ZaYgRz83/9JlHCrVpIywh3Tl9YQFycZOB0KHcR8ExIiM/rQUGjeXAK4mjSBfv0K\n2k2aJDP+e+6Rxebvv/duSooKJiI8kGevuxCAFxb/Zkw+BkMJ+JTWQGs9zo02HyFun87HjgA9yypY\nrefrr2V/883nR/leey3s3i3rAS+/LAr53ntlMbZZM2ljt0u8wOrVsq1dK+6ZACtWwBVXyGurVfIA\nKSX1d4OCxLafk1M4uMxmgx9+OH+doWVLiUS+//6CPqspN1wcwbc7jrH291M8/+M+3hp7cVWLZDBU\nS0yunqriyitlkXfMmPPPDRsmKSE2bhSvm/R0+Ne/pHjMvHlirmneXNJIP/QQfPedKP2GDWHIkMLV\nxKZNg6SkAlPP8eMykz91CtavLzzu0qXytOEo0O7nJ08RCxfKzcPBpk2wZo3cfKoRSilevL4HAb4W\nvt8Zz+oDJ6paJIOhWmJSNlR3tJYEbq+/Lq6gW7dKVa3OnSXVw5AhMHiwbG3beier57p1MrvPzhZT\nUdeucpNyBKKNHSsLwq1aSdGY228vXMOginlvzWFe/mk/rRoEsnzq5QT41tKymgaDEyYtc23g3Dmp\nxtWzJ9x5pyj9VavEE0gpUc7du8vsfMQIeRrwVirnyy6DDz6Q16+8Ijcf5+jjXr0kDiE2VqKVe/SQ\nG8/CheKhVMVMHNiOLk1DiE3K4MNf/6hqcQyGaodR/FXBjBliT8/JOf9cfLwEdLVqJfl+oqMl/88j\nj0hxeAdnz4pp5t13JVDs0UcLp4kuL7ffLhXKQBZ8Fy4sOPfEE5KYbu1a+OtfJX31mjXyVPDKK96T\noYz4Wi3M+Es3AGatOkRiSmYVS2QwVDO01tVu69Onj661HDmiNWgdGqp1dnbhc2fPal2vnpwHrQcN\n0vqrr7TOynLd1969Wt90U0H7evW0njFD6+Rk78k7c6b03aCB1ufOuW6TkqL1229r3a2b1rGxBcf3\n7NE6NdV7snjI5PlbdJvHf9RTv9hRZTIYDJUFsFW7qWPNjL+yWbxY9iNGSCbOkycLPGlCQmSx96ab\nZGF3zRrx+vHzc91Xt27iHbRtG1x1lSzAPvccXHih92oEz5wpC8QLFsjM3hWhoRKLsGePeAGBmHxu\nuEFMUP/5T5WYgJ66uht+VgsLdxxjx9EzlT6+wVBdMYq/snFkwrzySsnT07ateOU4mD1blLmzj31p\n9O4t/vtr10q93vHjC2rq2u3lU7pKiR3/yisLjhWX0M15jeH4cXEXPXlSXEF79oQllRvG0bphEBOj\n2gHwypL9JqLXYMjDKP7KJDVVFmiVkpn5M8+Iq+batQVtLOX4lwwcKD79zpXBPvlEMoL++KN3cgEt\nWSIeRc8/X3J/ERHy1PLNN1KvYN8+eSq55hpJV1FJTLm8A2GBvmw8ksTa3724BmIw1GCM4q9M5s8X\nF0mtJWq3Z08x5ziKrHsDpQoXc5k3D377Df7yF3HR3LatfP3HxMii9IwZcOutkFFCERSlJCJ53z5J\nQREaKk8mV1xRaaafsEBfplwu1dFeW3rAFGk3GDCKv/JYsULs4ACBgeKNs22bmGYqkqVLxTunQQN5\nGoiMlJKQZc2/P2WKlIisVw8+/1ziCEpL7eDvL15HBw+K3/+LL0oKiUrizgFtaRLiT/SxFH7aU8l1\nDwyGaogJ4KosUlMlL4/W4ho5ZEjljn/mDLz0kuQFys4WZbxihfjsl4XoaHmK+PNPaNq04CbgKa+9\nJusQjz5aOOLYy3yy8U/+8d0e2jcOZtlDg/CxmjmPoXZhAriqAzYbzJolNnyQGfKRI6KAK1vpgyy0\nvvqqpF0eP15uQpdcUvb+evSAzZvlsxw/XrrZxxVxcVLz94knpJ+YmLLLUwpjLmlFm4ZBHDmZxsLt\npdRANhhqOUbxVwSxsWLHfuABePjhguPh4VUnk4O2beHTT2H79gI30VOnRF7nRWZ3aNJESkw+9ZSU\nmAwM9Oz6li0lC2izZjL2RRfBV1951oeb+FotPDhUAuD+s/oQubbqlWfIYKhMjOL3Nl9/LQrsl1/E\nBOJIt7x4Mfz+e/Uptu7sk//66+JtNGiQ5OGJjXW/H6sVXngBhg4tODZrFmzZ4t71V10lZqMbbpA0\nFbfcAlOnuo5qLifX9WxBm4ZBxJxOZ3F0zUk5bTB4G6P4vUVmpqQvGD0akpPFbXH3bhg+XJTY2LHi\nBlkdc9zPmCGupQEBEqjVpYu4a3pqugHJIfS3v4lr6Zw57t3oGjUSt8+335ZF3zffhAcf9HzsUvCx\nWrhvsHj4/HvlIePhY6izGMXvDc6ehagoCb7y95fkaj/8IKYQkDTGqalwwQXQokXVyuqKoCCJ0N2/\nX25cGRlyM+jWTRaAPSEyUjx/srNh8mSYONG9G4hScsP45RdJAvfEE2X7LKVww8UtaREWwO8nUlm2\nz3j4GOom7lTgmquUOqGUclk2USk1WCmVopTambfNcDo3Uil1QCl1SClVMb/k6kBIiMzm27aFDRvE\ntu8cxbp8ueyHDasS8dymTRv48ksx+1x0kSy2OiKA3cXfX1I0zJ8vNv9588RzyN2grf79Zf2hdWt5\nr7XI4yX8fCxMyZv1v7PykInmNdRJ3JnxfwSMLKXNWq11r7ztOQCllBWYBVwFdAPGKaW6lUfYaoXW\nBRWvlJLZ/rZtcLGLqk8//yx757QH1ZnBg+WzLF0qJhsHCxe6b/65/Xa5CXboADt2SAoKd4PHnG+a\nL74oC8+PPiqeUl7glshWNA7xZ2/8WVYfcKOspcFQyyhV8WutfwGSytB3X+CQ1vqI1job+AIYVYZ+\nqh+ZmRIEdcUVBe6awcESJFWUlBQx9VitcPnllStnefDxkfUJB460yz16uJ9zp2dPKRwzbJj05zB9\neULz5nLtv/4lC+XOlcDKSICvlclR7QF4/5fD5e7PYKhpeMvG318ptUsp9ZNS6sK8YxGAs3tIXN6x\nms2JE+LB8tln4pe/d2/J7deskZnqpZdKyoKaSkCAZP08fFg8cUaPlrQTpREeLh5Na9dKjQFPmTgR\nli2TOIRFi2QtJbH8tvmxfVsR4u/DxiNJ7I5LLnd/BkNNwhuKfzvQRmvdE3gHcKSadFUOqliDqlJq\nslJqq1Jq68mT1fTxe+9eMVmsXy9KbN260oOg4uLE1l3NC5WXSr9+YrJ57TV5uvn6aynJ+MYbpefd\n8fUVk4+DV1+VIvLuMmSIJHzr2BF27pQ1g8Plm6mHBPgyrp+sI8xZa6p0GeoY7iTtB9oCe9xsGwM0\nAvoDS52OTwemu9NHtSzEsmSJFE8Brfv21Tohwf1rs7KkWElt4ehRrW+4oaAAzNNPu3/t3r1aKyXX\nzZyptd3u/rXHj2vdp49cO2KEx2IX5diZdN1h+mLdfvpiHZuUVu7+DIaqhMosxKKUaqaUrMYppfoi\nTxGngS1AJ6VUO6WUHzAWWFTe8aqEnTvFL//sWTFxrF4t0abu4udXs808RWnVShZ6f/xRTFgPPeT+\ntd26icePxSL1CGbMKP0aB02aiIfP3XfDRx95LHZRWoQHcu1FzbHZNfPWxZS7P4OhpuCOO+fnwAag\ni1IqTik1USk1RSk1Ja/JzcAepdQu4G1gbN4NKBd4AFgK/AZ8qbUuxSBeTenZU+rOPvEEfPGF+6kJ\nUlIqJAK12nDNNeK541jUzsoS+39pi7+33SZ/R0fUrydpqUNC4MMPC268Wkvkbxm5J2+R94vNR0nJ\nqMX/K4PBGXcfDSpzqxamHptN61OnCt57YpJw8OijWgcHa/3hh96TqzrzzjsF5p/bbtP65MmS28+f\nX9D+4489H89u1/qhh7T29dX6u+/KJrPWetzsDbrN4z/q99ccKnMfBkNVg6m5W05ycuCOOyR3TVKe\nJ6tytVZdCmvWQFpaQR3a2s6UKbJwGxAglb8uuEASwhUXJHX77fB//yevZ8woW4oIpeT/dfPNsuBc\nBiYOlPKM8zf8ic2kcTDUAYziL0pGBlx/vSiuP/+U6lVl4dw5iUC1WmHAAO/KWF3x8ZHC7NHR4sV0\n6pSYda6+uvjCLw89JLl51qzxPLunUuLf/9hj4lk0dqyYkDxkcJcmtG4QRNyZDFbtP+Hx9QZDTcMo\nfmfS0uDaa6U8YKNGsohb1kIl69aJ/35kZOFMmHWBjh0lWnnuXPG/X7JEcvAUx4MPSroIB9nZ7o+l\nlLiG/uMf8ve+9VaJsfAAq0Vx+6Uy/scbYjy61mCoiRjF7+DsWRgxAlaulGjRNWtEaZeVNWtkX5Oi\ndb2JUnDXXVJv95//lJm/A0e0c1G0liL0UVGemX2Ukmyizz4r1bwmTJCbtwfcEtmKAF8La38/xaET\n5Y8ONhiqM0bxgyiiYcNklt6qlSjtbuVMK1TXFb+DZs1g+vSCNZIDB2R2/5//iJJ2JiVF3DQ3bxaX\nTU8TqM2YIUVhLrgA+vTx6NKwIF9uuFgCy/+7IcazcQ2GGoZR/CC25QEDoF07MUl06lS+/tLSpBCJ\nxVI4yZlB8v2fOgX33y9J65zLLYaHSzrrkBCx1b/5puf9P/+8uJg2berxpRP6twXg621xnMs0rp2G\n2otR/CCz0TfeEGXdtm35+wsKkgXOzz+vXYFb3mDGDCmv2KiRBGN17w7vvlsw+7/wwoLgrGnTSl4b\ncIVSBWsqWsOTT0q6Bze4oHkofds1IC3bxrc7TF1eQ+2l7ir+kyelzN/x4/JeKWjY0Dt9KyV5bG65\nxTv91TZuvlls/6NHy9PRffeJqc3h+XPjjeKpY7PJ3zA+vmzjzJ8PL70EI0dK9LUb3DmgLQAfr48x\nufoNtZa6qfhPnxYzw1dficnBUPk0bixFX778Umb/69dL5K+DF1+U5GzHj0tlrrJw662SSjolRSKK\n3SgGM6xbU5qFBnD4ZBobDp8u27gGQzWn7in+M2dkdrl7t8zKZ83ybv+5ubKgO3Wq1wqH1GpGj5as\np199JVXMQEw0Z8+KnX/MGPj3v8vWt4+PBJBdcYWkch4xQp70SsDXamFsX0kf/enmo2Ub12Co5tQt\nxZ+SIj/+HTvE13zFijItApbI7t1il/7hBwneMpROkyYSP+Hg44/lJvDrr6L8mzcve9/+/vDtt5Jv\n6fffJb9QKcVcxlzSCouCZXsTOZWaVWJbg6EmUncU/7lz8ri/ZYt476xcWTGFz9etk31ZA78MUrjl\n9Gkx09x1l8z+bTZ45hm3bfWFCA2Fn36ShfstW2T9oASahwVyRdcm5Ng0X2+LK9NHMBiqM3VH8c+f\nL25+rVuL0i9LNSh3MIq//CxYAO+8Izl/PvpICr9PnSoBWjfeWJA/yROaN5cawtdcI/2Uwvi8Ii1f\nbD6K3eTvMdQy6o7iv+8+WTBcudI7LpvFYRR/+bFY4IEHxCTXp494+7zzjgSD/fEHjBtXtvWTzp2l\nhkDjxqU2vbxzE1qEBRBzOp0NR8wir6F2UbsVf26uLOaCuFg++WThEoDe5uhRKbUYHi7Ro4by0bWr\nPKU9/bTcDEJCxANo2TI5Vh5sNnmK+OQTl6etFsWYS2TW/9kms8hrqF3UXsWvNUyeLJGz7hQF9wbr\n18u+f39RVIby4+sr+Xt+/VU8f778UhbNX3pJ3peVxYslMnjiRCkE7wLHIu/SvYmcPGcWeQ21B3cq\ncM1VSp1QSu0p5vytSqndedt6pVRPp3MxSqlopdROpdRWbwpeIlrDo4/CvHliGjhaSTO2jh0lLmD0\n6MoZry7Rv7945gwZIgXfQUw+q1eXrb/rrpP4gOxsuOEGOHTovCbNwgK4omtTcu1mkddQu3BnWvoR\nMLKE838Al2utLwKeB2YXOT9Ea91La12OVJce8s9/SgoGX19x5evfv3LGjYwUn/O77qqc8eoqY8eK\nm6bNJou9P/xQtn7eeENqBZw+LYu+DrOgE7c6Fnm3mEVeQ+2hVMWvtf4FKNaNQmu9Xmvt+MVsBKq2\n3NSsWZKbXSmx344YUaXiGCqA5s0l6GvoUFHW110ni8GeVvDy8ZE4gR494OBBcR8tUgtgUOfGRIQH\n8ufpdNabSF5DLcHbhuiJwE9O7zWwTCm1TSk1uaQLlVKTlVJblVJbT5YSXVksn34qCgDg/fcrN1fO\nvn0SeFRcpSmDd+nQQRZ5X39dFPisWfLEtXevZ/2EhIinT7NmkjSuiKunLPKK6+9nm53+t1qL19HM\nmfLU0L69BAN27SoxIwZDdcadwrxAW2BPKW2GAL8BDZ2OtcjbNwF2AYPcGa/MxdYfe0wKd7/6atmu\nLw/PPitjP/hg5Y9d17n66oKi7ZMmla2PLVu0HjVK69OnzzuVkJyh209frDtMX6xPnM3UeskSrS+8\nsGBM561r18IXL1+udW5u2WQyGDyAyi62rpS6CPgAGKW1zn8e1lrH5+1PAN8Cfb0xXrG8/LKkYZg2\nrUKHcYnx3686nn9eSjwCJCR4VrrRQWQkfPcdNGhw3ilZ5G1Crl3zzfY4WVvYu1dcS++9VzyN9u+X\nLKI//1xw4aZNkhfKjRxBBkNlUm7Fr5RqDSwEbtdaH3Q6HqyUCnG8BoYDLj2DvIZSkpCrsrHZxN8c\njOKvCnr3FoVbv76Yba66CmJj5buwapXn/eXkwBNPFJjtfvmFx3d/D8CCLbHokSNh4UJR9P/5j3hx\ndekiaw8REQX9pKdLHqIVKyrXrdhgKI3SHgmAz4EEIAeIQ+z4U4Apeec/AM4AO/O2rXnH2yPmnV3A\nXuApdx9DymzqqSp27pTH/LZtq1qSus22bVo3bSr/i0aNZK+U1tOna52d7X4/06fLtV26aD1smNag\n7RaLHvvgh7rN4z/q9YdOud8N7d7lAAAgAElEQVRXXJzWPXtKf507uzQlGQzeAG+aerTW47TWzbXW\nvlrrllrrD7XW72mt38s7f4/Wur4Wl818t02t9RGtdc+87UKt9YtevF9VL4yZp3rQu7eYV3r0kPKO\nvXvLU+BLL0mqbHfjOUaMgOBgqQ+8fDmEhKBmzOCywb0A+NyTdM0REZIm5KKLCjyHckxZR0PVYsJL\nvYFD8Q8YULVyGKSQ+/r1YqpZvVq2li3FFNezZ8k+/1pLdbDBg6UymIOpU2HmTK6/vCtKwZI9iZxJ\n82AdoUEDMUE1by7yPPNMmT6aweAtjOL3BjYb+PmZGX91oV49meWHhEBUlKRi9veH5GTx+b/3Xnjv\nPXH5/cc/Cmz5SkFgoGxTp8KHH8rx55+HJUtoWT+IQZ0ak22zs9DTmrytWkkN5latJP7AYKhClK6G\ndUUjIyP11q2Vl+HBK2RlSaSwydFT/YiNlZn85s2uz3//vdwQQG4CQUEFGTyfeUZ8++vXh23bWJIW\nwJRPttOpST2WTR2EUsozWbKy5CZkMHgZpdQ27WaGBJ+KFqbOYH7M1ZdWrcT2f+iQmIGioyXI6sQJ\nuVl37FjQtk2bwtfOmAHbtsnTQHg4Q1uH0aieP7+fSGX70TP0aXO++2eJOH9PkpJcuo8aDBWNUfzl\nJTFRZoemzGL1p2PHAiWfni65/vfvl4jb115zffO2WCStQ2AgWCz4AqMjW/Lu6sN8vjnWc8UPBSmh\nP/oIfvutsAuowVAJGLtEebnuOpm1bd9e1ZIYPCEwEKZMkRn/O+/I+szhw67bBgcXmPCys5mgEgD4\ncXc8ZzPL4KFjtYpP/7lzUhzIYKhkjOIvD+npkq8lNRU6dapqaQyeoBQ8+KB4ZLVrJ+ac3r1LzvGf\nlgaXX07z669mrF8SmTl2vvd0kdfBCy/IzWTOHEkdbjBUIkbxl4ctW6TK10UXiQeJoeZxySXytHbT\nTVLU/ZZb4JFHXLcNCpLKapmZPP3RDEIzU/l8cyxlcpC44AK49Vb5/jz3XPk+g8HgIUbxlwdHxS3j\nv1+zCQ+Xmf6sWeKW27On63ZKSZuePQmOjeHNpW+zLz6F6GMpZRv3mWcks+j8+WbWb6hUjOIvDw7F\nb/z3az5KwX33SXTthAkFxxMSCrcLDIRvvoGwMK7Yv57Jmxfy+ebYso3Zvr1UEbPb4a23yi67weAh\nRvGXFbvdzPhrI87unNHR4gX05JNiknHQoYPM0oHH13xM/Pc/kZaVS5l45BG56ZjsnYZKxCj+snLw\noPhhN29+vu+3oXaweTNkZkoU8PDhcPx4wbnrroMnnsCq7YzcuYIfd8eXbYyePSVo7NNPvSOzweAG\nxo+/rHTqBDt3iinA0+hNQ81g4kSZ3Y8dK+mdL74YFiyQNBAAzz/P5notmH62Hb02xzLmktZlG6dV\nK+/JbDC4gZnxlxWrVWZrI0uqQ2+o8QweLC67gwbJTX7IECn3qDX4+NBj2n2EBPqyMzaZ3+LLuMgL\n0t+aNVDTUpUYaiRG8RsMpdG8uRRTeewxibqdORNiYgAI9LNyw8URtEw5TtCVVxSs+3jK7NlykzGZ\nOw2VgFuKXyk1Vyl1QinlsoKWEt5WSh1SSu1WSvV2OneHUur3vO0ObwlepSQlyeP+zJlVLYmhsvDx\ngVdekfKMc+dK0FceYy9pzdhdS2nz23b0uHHy/fCUG2+UKOKffnK/boDBUEbcnfF/BJRk07gK6JS3\nTQbeBVBKNQBmAv2QerszlVL1yypstWHjRvj1VymwYahbjBoFY8YUvJ8/n24/fcWqW6aws3ln1NGj\ncNddYrrxhMaNRfnb7fDBB96V2WAogluKX2v9C1DSNGYUMD+vAthGIFwp1RwYASzXWidprc8Ayyn5\nBlIzMIVXDABxcfDXv8I99/Daind5+OqppAXWg0WL4O23Pe9v8mTZ//e/cgMwGCoIb9n4IwDnKJa4\nvGPFHa/ZmMAtA0hlr/feg8BA2v/4Fe8s/hcvRuUFf02b5vlC7eDB0mdMTMHkwmCoALyl+F35M+oS\njp/fgVKTlVJblVJbT1bnYJacnIKCHv37V60shqrnjjukrGP79lyYeIhpv/yXA/2GyPdkzBjJwOku\nFovk7wH45JOKkddgwHuKPw5wdkZuCcSXcPw8tNaztdaRWuvIxo7qR9WR3bslK2enTgVVmgx1m549\nYetWzl4xjPqZ5+i0aRX2tm0lBURwsGd93XabeBG1bFkhohoM4D3FvwiYkOfdcymQorVOAJYCw5VS\n9fMWdYfnHau5mDQNBlfUr0/Isp/4ZMSdpPkF8surH0g6Bk9LcXbvLqUin366YuQ0GHAzclcp9Tkw\nGGiklIpDPHV8AbTW7wH/A64GDgHpwF1555KUUs8DW/K6ek5rXQZft2rERRfJItyIEVUtiaGaoaxW\n7DNnMrjTlXRNCWaw48SWLVKsp0MH9zoy1dwMFYwptm4weJGUjBz6/fNnMnPsrJk2mDavvwAvvyx+\n/7/95n5t5lOnJD3E+PFS6N1gKAVPiq2byF2DwYuEBfpydY/mACzYEisunyD59gcOLJzlsyTGj4cH\nHhDXUIPByxjF7wnbt8Pnn0N8GTMxGuoE4/pKsravtsWRM+8jKawO4t4ZGQknTpTeyY03yr6kUpAG\nQxkxit8TPv5YZmIffljVkhiqMZFt6tOxST1Onsti5YGT8MYbcM89cnLXLujVCzZtKrmTG26QheFl\nyyClHMnfDAYXGMXvCcajx+AGSinGXiJezJ9vzsu78957BXEfCQkS8VtSdG7TppIRNCcHfvihgiU2\n1DWM4neXtDRJz2uxQN++VS2NoZpzY++W+FktrDl4kmPJGeKp8+WXslBrscA//lG6q+fNN8vemHsM\nXsYofnfZulVS8vbsCSEhVS2NoZrTINiPEd2boTUscMz6W7aEzz6TSF+HUtca/vlP1xk5b7xRivws\nXQpnz1ae8IZaj1H87vLrr7I3aRoMbnJbP1nk/WzzUbJybXJw5MjCT4zz58NTT8mi75o1hTto3lzy\n90RFubcgbDC4iVH87rJ2rewdZfcMhlLo264BFzQP5VRqNot3J5zf4PvvxVNs+HAptj50qGT1dI6t\nWb5cto4dK09wQ63HKH53sdulGIdR/AY3UUpx14C2AHy0PoZCwZLx8ZLE7e23YdIkePxxMSU++KDk\n88/MlHYmitdQARjF7y4Ot7qImp9V2lB5XNerBfWDfNkdl8KO2OSCEy1aSEUvkBQg998PX3wBQUHi\nNhwVJd4/IDeE9esl8tdg8AJG8XtCUFBVS2CoYQT4WhmbF9D10bqYwif//ne4+mo4cwZuv10WfDds\nkPQOWVkQGirtXnpJaj+8807lCm+otRjF7w5xcaYikqHM3HZpGywK/hedwPGzmQUnlIJ588Rnf80a\nePVVSQK4ZYv47jtSOl95pewXLfK8pKPB4AKj+EtDa+jdW3LvJyZWtTSGGkhEeCAjLmxGrl3z6aYi\nbptNmohpB2DGDCny07AhtGkjx7SWmX5gIBw7JovBBkM5MYq/NA4eFI8LX1+ZmRkMZeDOvEXezzb9\nWeDa6WDECMnnExFx/pPlwYPw7beQkSHv//vfihfWUOsxir80fvlF9lFR8mhuMJSBvu0a0LVZSPGu\nnS+9BDt3wqWXFj7epYvY/Zs0kfezZhV8Jw2GMmIUf2k4/PcHDapaOQw1GqUUd13WFoC56/7gvDoY\n/v4QHl7w3jlgq2dPSRditUpa5yuugH//29j7DWXGLcWvlBqplDqglDqklHrCxfn/U0rtzNsOKqWS\nnc7ZnM7VvOTiJnDL4CVG9YqgYbAfe46dZcPh064b5eaKL3+XLhATU3C8RQsYNUoWfG02+NvfTJZY\nQ5kpVfErpazALOAqoBswTinVzbmN1nqq1rqX1roX8A6w0Ol0huOc1vo6L8pe8cTGyo8vNBR69Khq\naQw1nABfa76t/71fjrhuZLXCn39CcrIUXncu3PLBB5Kz59NPxb1z/PiKF9pQK3Fnxt8XOKS1PqK1\nzga+AEaV0H4c8Lk3hKtyHPl5LrvMRFAavMJtl7Yh0NfKLwdP8luCi8RrSomCb9EC1q2DF18sOOfI\n7Dl+vDyJOuJKUlNh48bK+QCGWoE7ij8CiHV6H5d37DyUUm2AdsBKp8MBSqmtSqmNSqnryyxpVTB6\ntPhUP/NMVUtiqCXUD/ZjTF6u/tnFzfobNZLkbUrBc8/JDcCZ06cLsnlqDXffLabIWbOM3d/gFu4o\nfleuLMV9u8YCX2utnf3VWucVAB4PvKmU6uByEKUm590gtp48edINsSoBHx/Jmmjy7xu8yMSB7bBa\nFD/sipdc/a4YOhQee0zcO2+9taAK12efiYfPP/4h7+12aN1aTEIPPAATJxbk+TEYisEdxR8HtHJ6\n3xIorujsWIqYebTW8Xn7I8Bq4GJXF2qtZ2utI7XWkY0bN3ZDLIOhZtKqQRDX9GhOrl0z99c/im/4\n3HPQp4/Y/B2Kvk8fUfb/+58oe6sVXn9d7P6BgRIJfPnlBUXeDQYXuKP4twCdlFLtlFJ+iHI/zztH\nKdUFqA9scDpWXynln/e6EXAZsM8bglc4ixbJD+iTT6paEkMtZPKg9gB8sfkoKek5rhv5+ckM/9Zb\n4dln5ViXLtCpEyQliX+/g/HjJZFbmzYS/RsZWbBGZTAUoVTFr7XOBR4AlgK/AV9qrfcqpZ5TSjl7\n6YwDvtCFHZQvALYqpXYBq4CXtdY1Q/EvXSqBMs4udQaDl+geEcbAjo1Iy7bx8YaY4ht27iyTjwYN\nCo795S+y//HHwm179ZJKcUOHwvHjsk5gMLhAnRdIUg2IjIzUW7durVohLrgA9u+XhTVTXN1QAaw/\nfIrxczYRFujLuieuoJ6/T8kXZGfDu+9Ct25SvOWCC2Cfi3lUbq4EeE2ZAgEBFSO8odqhlNqWt55a\nKiZy1xXx8aL0g4PhkkuqWhpDLaV/+4Zc0rY+KRk5fLw+pvQLxoyBhx6ClSshLEzy8x8+fH47Hx9p\n51D6Z8/ChAmS5M1gwCh+16xaJftBgyQ5m8FQASil+PvQTgB8sPYIaVm5JV/w8MPix//KK7LIq5TY\n80vjscckuVtk5PmuoYY6iVH8rnAo/iFDqlYOQ61nYMdG9G4dzpn0HD7Z+GfJjaOi4MknxVd/3z55\nKh03rvRBnn9eirYnJsp3+v33vSK7oeZiFL8rVubFn11xRdXKYaj1OM/6Z/9yhPTsUmb9M2ZIBs/E\nRJg+3b2ArcaNpXTogw9CTo7Y/idPlipfhjqJUfxFsdvlx3XPPeIlYTBUMJd3bkzPlmGcTsvmvxtK\nmfX7+orPfkgILFxYkL+nNHx94c03peiLvz/MmSOz//R073wIQ43CKP6iWCxw553ywzD5eQyVgFKK\nqcM6A/Cf1YdJySjGr99B+/bi3QMyi7//fvcHmzBB/PtbtpR0z6aOdJ3EKH6DoRpweefG9GvXgJSM\nHOYUl8PHmVtvlclJRoZE8dpspV/jIDJSSji+9VbBseTk4tsbah1G8TujNTz+uDxCe/JDMhjKiVKK\nx0Z2BeDDX//gxDk38u1MnAgdO0oU7/r1ng3YuLFEBoMo/UsugXvvlVgBQ63HKH5n9u+HV1+VH4DF\n/GkMlUufNvUZ1q0pGTk23llxqPQLlCqI4h0zBn7+uWwDb94stSfee08cGhITy9aPocZgtJszS5fK\nftgwU1/XUCVMG9EFi4LPNx/lz9NppV9w7bWyT0gQ+/2pU54POny45PePiBA//8hI2LTJ834MNQaj\n+J1xKP4RI6pWDkOdpXPTEG7s3ZJcu+bln/aXfkFUlETxgij/iRPLlpP/kktg2zYYOFAifAcNgrlz\nPe/HUCMwit9BZiasWSOvhw+vWlkMdZpHhncm0NfKT3sS2XikmNq8Dnx9YeRIeR0QIFll33uvbAM3\nbQorVsB994mtf+JEM/OvpRjF72DtWvGQ6NVLfgAGQxXRPCyQKZdLvaJnf9iHzV7KDP6xx+T7O2+e\nvH/4Ydi7t2yD+/lJJa8PP5R++/UrWz+Gao1R/A6MmcdQjZg8qD0R4YH8lnCWBVtiS27cu7eYaMaO\nhbvukqfXcePKV4nr7rslJ5CD3bulDKmhVmAUv4PWrSWgxfHYbDBUIYF+VqZfLe6dry87UHpQl4O3\n3oLu3UVxO9w1y0tSEowaJesJH33knT4NVYpR/A7+/nfYuVOSWRkM1YBrejSnb9sGJKVl83/LD5bc\nePt2WZuaOhV27JC0zN5ySa5XTyZEWVnyRPH3v0vOH0ONxa1vhlJqpFLqgFLqkFLqCRfn71RKnVRK\n7czb7nE6d4dS6ve87Q5vCm8w1GaUUsy8rhtWi+LjDTHsii0hujYwEJYvh+++K+yKfORI+fPw+/lJ\niojZs2Ux+Z13xOX5xIny9WuoMkpV/EopKzALuAroBoxTSnVz0XSB1rpX3vZB3rUNgJlAP6AvMFMp\nVd9r0nuLRYvgkBsBMwZDJXNhizAmDmyH1jB9YTS5Nrvrhl27QocOcPp0QS3eX34R+//YsVKVq7xM\nmiSeb82byz4yUlxADTUOd2b8fYFDWusjWuts4AtglJv9jwCWa62TtNZngOVA9TKiZ2VJ3pNOnSAu\nrqqlMRjO46ErO9GyfiD7Es4yd90frhs5R/H+8IPsu3aVKnK//ioZZ71B//5S17d/f4n2/f137/Rr\nqFTcUfwRgLNbQVzesaLcpJTarZT6WinVysNrUUpNVkptVUptPXnypBtieYnVqyE1VRZ2W7asvHEN\nBjcJ8vPhheu7A/DG8oPEJhWTSrmo4m/SBD7/XGz9L70ES5Z4R6AWLaRY0TffyNOEocbhjuJ3lbug\nqGPxD0BbrfVFwM/Axx5cKwe1nq21jtRaRzZu3NgNsbzEokWyv+66yhvTYPCQwV2a8JeeLcjMsfP4\nN7uxu/Ltd0TxOtfiHTQInntOXt9+u/fq7vr7w403FrzfuhWuvhoqc9JmKDPuKP44oJXT+5ZAvHMD\nrfVprbWjnM8coI+711YpWhcofsdsyWCopsz8SzcaBvux/vBpPnJVnN05ivfHHwuOT58ui7GnTol/\nvzfs/c5oLTUBfvqpIOWzoVrjjuLfAnRSSrVTSvkBY4FFzg2UUs2d3l4H/Jb3eikwXClVP29Rd3je\nserBzp1i12/eXIpXGwzVmEb1/Hnpxh4AvLJkP4dOnDu/0ZQpkrLhllsKjlks8Mkn8j3fuNH7aRiU\nklTm/frB0aNw2WVSJcxQbSlV8Wutc4EHEIX9G/Cl1nqvUuo5pZTDPvJ3pdRepdQu4O/AnXnXJgHP\nIzePLcBzeceqBw5b6F/+YtIwG2oEwy9sxug+LcnKtTN1wS5yinr5DB4Mf/2rKHlnmjSBr76Shd7L\nLvO+YBER4ulz990SMXzbbfDII95/ujB4BaXLksmvgomMjNRbt26t+IGeeEJ8khcsKEhvazBUc85l\n5nDVW2uJO5PB/UM6MG1E16oWqQCtxef/wQdF6V9/PXz7bVVLVSdQSm3TWke607ZuT3NfflnsniY/\nj6EGERLgy79G98SiYNaqw6w6UCSQ6vhxSdR2553Fd7JwoTg0eHtGrpRk91y5UpIdTpjg3f4NXqFu\nz/gNhhrMv1f+zuvLDhIe5Mviv0cRER4oJ5KSROlqLdW0GjUqfGFaGnTuDPHxkoHTORmbN0lNlXQP\nDmJioG3bihnLYGb8brFmjQRvGQw1lPsGd2Rwl8Ykp+dw/6fbyc7Ns/c3aABXXil1o12ZWYKDxb/f\napVSowsXVoyAzkp/40bo0gWmTTN2/2pA3VT8MTGyCNa+vSmqbqixWCyK/7ulFy3CAtgZm8wLi/cV\nnLz5Ztl//bXriwcNEqUPYhLa70a1r/Kwfz/Y7fD665JMzuT5qVLqpuJ3/BiiomTWYzDUUOoH+zHr\n1t74WS3M3/An8zfEyInrr5fv9ooVkr/HFVOnitvnuXMSjHXOhXuot7jzzgK7/6pVkkNo48aKG89Q\nInVT8X/5pexHj65aOQwGL3Bx6/q8fJP49z/7wz7WHDwJDRvC0KHyRPvdd64vVEoqbXXrJtG+995b\nsYJGRUlw12WXFdT1fffdstUINpSLuqf4Y2KkklBwMFx1VVVLYzB4hRt7t+SBIR2x2TUPfLqd34+f\nK5jYOCY6rqhXT2z83btLDv+KpkULmfk7cvpPny5eSIZKxaeqBah0vvhC9tdeC0FBVSuLweBFHh7W\nmSOnUvlfdCJ3zN3MN2NG0vyvfy3dpbJLF9i1q/KCGP38pFJYv34QGgrNmlXOuIZ86taMX2v4OC9/\n3K23Vq0sBoOXsVgU/xrdiz5t6hOfksmtCw9x6vW3YMAAdy4ueP3NN5WTonz8+MKBk3PmwOLFFT+u\noY4p/lOn5AvetKmprWuolQT6WZl75yVc0DyUI6fSmPDhZvfr9QJ88IF4BI0eXbnuznv2SODXtdfC\nzJnG266CqVuKv3Fj+YJt3y6ZDA2GWkhYoC/z7+5Lu0bB7Es4y3t/e5WcoVeKOac0Ro2CVq3E4+a+\n+ypv4bVbN0kfbbHI/tprJRDNUCHULcUP4snQokVVS2EwVCiNQ/z578S+RIQH0mL7enxXriBj7kdu\nXNhYvIACA2HuXLHFVwYWiyz0LlkiHklLlkjG3B07Kmf8OkbdUfzR0QXFKSoYu93O8ePHKS0dhqft\nbDZbfnvna4ued25X0pjujJ+bm8uePXuw2+2l9uEsR2xsLCtXriQ+Pj6/f+e+HK9zc3Pzr8/OzmbV\nqlXnXbN7926OHTtGfHw8cXFx7Nq1i2PHjuW/TkhIyG9vt9tJSEggPj4+fyvaztEmMTGRrKwsVq9e\nTVZWFrt37yYhIQGbzZbfR2JiIlrrYuVwtD9+/Dg5OTns3LmTnTt3FjruPGZR+aOjo7Hb7aSnp/PB\nBx8QFxdX4v80JycnX87i/m+OayPCA/lySn82DLgagPR58zmRlFrs/zX/2ogI9Ny5cuCRR2CpdzOp\nl/i9GzZM6vhGRooHXv/+xbujGsqO1rrabX369NFeZ+RIrUHrL7/0ft9O2Gw2PWjQIO3j46MHDRqk\nbTabV9pZrVYdFhamfXx8dFRUVKFro6KitI+Pjw4LC9MWiyW/XdF+nccs2oer8XNycnRYWJgGdFhY\nmM7JyTlPbsfYzv2FhIRopNKaBvTAgQN1VlZWfl+hoaE6NDRUA9pqtWqr1aoHDhyorVZrsdeUtg0c\nOFDn5OToQYMGldguKipKDxw4UANaKeWyTb169Qqdu+yyy/LldbWFhoZqpVQh+QEdEhKirVarjoqK\n0lFRUSWO5/z+sssuK/Q/dfzvi44xcODA8/5vrr5Xx5PTdWyjllqDfvyOF/SB+OTz/q+urrU/9ZT8\nZsLCtN6/v6w/iTJ973VGhtaTJmldv77Whw97ZezaDrBVu6ljq1zJu9q8rviPHNFaKa39/LQ+edK7\nfRchMTFR+/j4aED7+PjoxMREr7VzbFartdC1RRWOYyvar3NfRftwNX50dHSh/qKjo8+T2zG2c39F\nN4vFoletWlWiQrZYLB5fU7R9dHR0sTI4/+2KjlXaVpzCdndz3Nw8+SzutLdYLOf934r7XqU+/0+t\nQa9s30df+PT/tH+r7oX+ry6vjY/X+oYbtG7bVus9e8r0eyjr9z6f2NiC13a71sePe0WO2ohR/EWZ\nNk0+6u23e7dfF9jt9sKzJrvdK+3KOuN37td5zKJ9uBrfZrMVmhnabLbz5HZnxh8VFaVzc3M9nvE7\nX1Pa5pj9VtcZ/6BBg9ye8TueENyZ8Rf9vxX7vTp5Utv9/bVNKT3wrx/o1o9+p4O7D83/vxZ77blz\nXp0sufu9d8m//qV1w4ZaL1niNXlqE15X/MBI4ABwCHjCxfmHgX3AbmAF0MbpnA3Ymbctcmc8ryr+\n9HStGzSQj7pxo/f6LQGbzaYTExNL/VJ72i43Nze/vfO1Rc87tytpTHfGz8nJ0dHR0eeZjEqT4+jR\no3rFihX62LFj+f079+V4nZOTk399VlaWXrly5XnX7Nq1S8fFxeljx47p2NhYvXPnTh0XF5f/Oj4+\nPr+9zWbT8fHx+tixY/lb0XaONgkJCTozM1OvWrVKZ2Zm6l27dun4+Hidm5ub30dCQoK22+3FyuFo\nn5iYqLOzs/WOHTv0jh07Ch13HrOo/Lt379Y2m02npaXpOXPm6NjY2BL/p9nZ2flyFvd/K/b/OmGC\n1qB/vWmibvP4j7rN4z/qJxfu1pk5ue5/J73wG3L3e18Iu13ra6+V37FSWj/9tNa5uaVfV4fwquIH\nrMBhoD3gB+wCuhVpMwQIynt9L7DA6Vyqu8I4Nq8q/rlz5WNGRsqXx2Coq2zfrvVbb2l95oz+bNOf\nutOT/9NtHv9RX/fOWh2blFb69Y89Jr+l99+veFldYbNp/dxzovhB66FDtS7NVFSH8ETxu+PV0xc4\npLU+orXOBr4ARjk30Fqv0lqn573dCLR0o9+KR+sCd7T77xdXToOhrnLxxZIjJzyccX1b8829A2hZ\nP5BdcSlc9dZavttxzDFZc02XLrK/7z746afKkdkZiwWefhqWL5cawitWyGf65ZfKl6WG447ijwBi\nnd7H5R0rjomA87ciQCm1VSm1USl1fXEXKaUm57XbevLkSTfEcoP0dPlitG4NY8d6p0+DoTaQnk6P\niFB+/NtArrygKecyc3lowU4e+GwHZ9KyXV9z993w1FMSVTt6dNX52A8dKmNHRUFCgtTOLumGZTgP\ndxS/q2myy7+yUuo2IBJ4zelway3lwMYDbyqlOri6Vms9W2sdqbWObNy4sRtiuUFwMMybBwcOQECA\nd/o0GGo677wjk6HlywkP8mPOhD68clMPgv2sLI5OYNj/reHbHXGuZ//PPy95rtLS4Jpr4OjRypcf\nCrJ8zpgBn31mnuY9xB3FHwe0cnrfEogv2kgpdSXwFHCd1jo/yYfWOj5vfwRYDVxcDnnLhlH6BkMB\n6elSnOWllwBQSjHmksUSr2gAABUvSURBVNYseWgQ/do14FRqNlMX7GLs7I2S3tkZRw7/wYNltn31\n1ZCcXPmfAcDHB559tqCOr9bwwAOSdt1QIu4o/i1AJ6VUO6WUHzAWWOTcQCl1MfA+ovRPOB2vr5Ty\nz3vdCLgM8f6peB55BL7/3jwCGgxFmTJF0iGvXl2oClarBkF8MflSXh/dkwbBfmz6I4mr3lrLsz/s\n5XSqU8I2f3/J4X/BBVI/NyWl8j+DK+bPh1mzpNDLv/9tfvsloEpczHE0Uupq4E3Ew2eu1vpFpdRz\nyCryIqXUz0APICHvkqNa6+uUUgOQG4Iducm8qbX+sLTxIiMj9datW8v2iUC+zP37i6knJgYaNSp7\nXwZDbeTJJ2XGf+WVslhahOT0bF5beoDPNh9Fa6jn78OkqPZMjGpHPf+8Mh5Hj8pvrGHDSha+GLKy\n4NFHRemDlJX84AMICalauSoJpdS2PLN66W3dUfyVTbkV/9ChYv978kl48UXvCWYw1BZOn4YOHWS2\nvmyZ5Mhxwd74FF5beoDVB8ThomGwH3cPbMdtl7YhLNApw63W0s/w4VVvb1+wAO65B1JToXNnqbHd\no0fVylQJeKL4a1+Stp9/FqUfFiZ3f4PBcD4NG4o3DMDjj0ORZG0OLmwRxkd39eWLyZfSu3U4p9Pk\nSWDASyt4cfE+ElIypOG0aVLj4p//rKQPUAJjxkiitx494OBBqfRlCrsXonbN+HNyxH1z7175Ak6f\n7n3hDIbaQno6dOokHj6LF0ODBiU211qz7tBp3l1ziHWHTgNgtSiu6NqEB5N2cOHDf0VpDf/5T8UX\nbneH9HRZ7N2zB379VUo+1mLqrqnn7bfhwQehfXtR/sabx2AomYMHoWNHj+vtRsel8N4vh1m6J5Fc\nu+iQew+u4PFv/w+tFOqzz6pP7ExqqhSVBzFx/fkn9O5dtTJVAHXT1GO3S81OgDffNErfYHCHzp3L\nVGS9R8swZo3vzfrpVzBtRBda1g/k3c5DeXXQBJTW5N56G0tem0tiSmYFCO0hDqWvtQShXXqp6Ihq\nOOmtLHyqWgCvYbGIHe/LLwsXcDYYDKVz9KikNbnjDqm56yZNQgK4f0hH7r28A2sPnWJhz+bMy07l\nro0LGTJ9ChN2nCBzwEAGd2nCkK5NuCgiDIulyOKv1vDbb7B2LezeDYcOiSu2Y/J2772ydtewodyo\nunWTQi0DBkBQkPufMTdXzFo5OTB1qqwFzptXfbySKpHaZeoxGAxlY9YssYeHhMjCaKdOZe4qIyuX\nE7dPxHf9r9x88/PEB4Tln2sY7MflnRvTv00Yg+L30mTJItT338OJE4U72btXFDzAVVdJKcai+PrC\nuHHw8ceeCfjttzLzT06Gli0l8jcqysNPWf2ouzb+CiQjI4OAgACUUmityczMJDAwsMxtPenPkzFK\n6jcpKYnw8HCSk5MJDQ3l7Nmz+e8b5C3sZWRkoLXO70trjcViISAg4DwZk5KS8Pf3JygoCK01f/zx\nB2FhYQQFBeHv709ycjKBgYFkZmZSv359MjMzSUtLIzMzE39/f7TWZGVlERAQwNmzZwkNDSU4OJj0\n9HSSk5Np3rx5/jgWiwV/f38CAgKwWCykp6eTlZVFixYtiI+Px2azkZOTg4+PPMRmZ2ejtebMmTM0\na9YMpRQpKSnY7XYaNmxIRoZ4ozRs2JDU1FSUUgQFBeHr60tqaiphYWEkJyeTmZlJ48aNyc7OJi0t\njdTUVJo2bYrWmhMnTuDn55fff0xMDOHh4WRnZ9OwYUNyc3Ox2Wzs2bOHDh06EBgYSHx8PFarlaCg\nIEJCQsjNzcXPz4/MzEwiIiJIT0/nzJkzBAYGEhAQwJkzZ/D39ycrKwu73Y7Kc5XUWpOTk0NISAjZ\n2dmkpKTg6+uLxWLBarWSlpZGaGho/udy/B8CAgKIiYkhIiKCY8eOER4eTmxsLOFhYTT529+o//PP\nZHfpwqEPPySkdWsaNGhAYmIiISEhhISEkJqaSlZWFg0bNsTPz4+UlBQCAgLIyMigQYMGZGRkkJSU\nBHY7zYODOWP15dcDx1kXk8L6P85yLFn+7j0SfueH+VPzv0vpjZqSeekAfC7pjd+FF+B/5ZUQGkpG\nRgbqzBkCU1Ph+HHYv18Wa9evh+3bxSMpLwKZ06fh++/JGDWKgAYNSvxtZR44gP+dd6I2bkRbLOQ+\n/TS+M2dWvStqOTCK38tkZGQwZ84cOnfuzIgRI1i6dCkHDx5k0qRJLhV6aW096c8TeW677TY++eQT\nl/1mZGQwa9YsGjRokK9IbTYbDRs2JCkpifvvv5/AwEDef/99MjMz6d69O9HR0dhsNurVq0fHjh05\ncuRIvoxJSUnMmjULpRR9+vThyJEjnDp1CgBfX198fX1JT0/Hx8cHu91Ojx49+OOPPzh79myJn83X\n15ecnJz89xaLpVBdWKVUvuwAHTt25NChQ+79I13gUA7uHq9IunXrxv79+/M/b9HPXh58fHyw2Wwo\npYrt0y8zk0lz5tDo9GliW7bkk9tvJycgIP/vUK9ePVJTpWZvSEgINpuN9PT0/P9Zz5492bt3L7m5\nuQD5Nxtts3Hjt9/S7qKLeO3iy0gPbcNZ/8ZMePavbGvamf91Hcj+xm1BKRSaBpZMujQJpHWIIiPh\nEM2C4KF7Xfw2kpMlYZzDVPPGG/DII2QFBpIwfDhtXnqJpX/+Wezvr0v79gxfuxb16qvsGDyYrv/7\nn0eTr+qGUfxeRmvN0qVL2bRpU/6xfv36MWLEiPwZmCdtPenPE3mGDx/OsmXLXPartWbOnDkkJiae\n11+zZs2YNGkSSimWLFnC5s2bXY7rLKPdbmf27NkcP368VHmbNm3qVjtDxeLO/yE0JYW75s4lPCWF\n2JYt+WLcONKDg8s0nm9WFhft3k3U2rWEnT2LBnY+/jiL8pSr1uDbsCV7TmRxwh7MaXswKToA7SIv\nZMNgPzo2qZe/tW0UTKv6gbSsH0SAr1Uaff89+uWXUU4++wc6dyZ50iT6PvwwKm8Ru+jvp3VMDC1u\nvJHh114rvz9nL6AaRN306qlAlFKMGDGi0LHilLQ7bT3pzxN5LBZLsf1aLBYmTZrksr9JkyZhsVhQ\nSjFy5Mhix3WW0WKxMHnyZLfkdbedoWJx5/9wNiyM+XfcQUpoKK3i4uiyf7/H4zQ8dYqRP/3Ew2+8\nwbWLF4vSDwhAAb3eeov2eU9oSsH0+++iq89JBvnFcEPAXqKfGc7Vfvvp53uUDtZTNFJpBPlZOZ2W\nzaY/kvh001Ge/WEfd83bwpVv/ELXp5cQ+cLPXD9rHQ+kteKVGfP4bu6PrLp4ANk+PnQ5eJB+06bJ\nwnUeRX8/R9u2LVD6SUnQvTs89hhkF5OeuhZgfeaZZ6pahvOYPXv2M9VJWThmCMeOHcs/lpmZSYcO\nHYqd8ZfU1pP+PJGnffv2LFu2zGW/jhm/41HdmYMHD9I7z695yZIlha53xllGx4w/LS2tVHkPHDjg\nVjtDxeLu/yEzMJA93bsDsO6y/2/v/oOkrus4jj9fu8fdIeLBdUDAAVLJTahwhhEOQ+AP8hKHxsYK\nyx/9GLUpzGak0pwJyBprKs3GZvphTDOClD9SEUrE0XJSI1AOBe0uwAQK2gPuPIFb8G7f/fH5Hrec\nd9yPuPsut+/HzA733f3sfl9z7L7v+/18P/v5zDzW7/3ePXs4NGTIu/vBzY7dV/DOO9x8zz2M37WL\ngpYWdo4bx9Nz5/LC9ddTOWoUiY0bOXvrVupGjmTfiBHU1tYe9558Y/s2OHyAEYlDTEg2UFGwj89/\nZAzfXnAhsytGcu7YEsYOG8ywwYOQ4PDRZg4eaWZvY5ra/x5k45v1PLkfHn//LB6YWkW6oIhJ+3Zy\ny/iL+eEbxuOb/s1f/7qFFX/ZQm1TManM6ezPnMarb9aRLhjCO0+spuzB5ej55zmyajUHps3g6LBS\npPBlte4enPVWc0sGiV7tZ+nSpXuWLFnyq+60HTjDOftQOp2mtrb2WFdHa9/57Nmz333RqBtte/J6\nPckzffr0Tl+3qamJVCpFWVnZu/r4U6nUsQuxNTU1FBUVddjHn52xoaGBuro6kslkl338dXV1TJ06\n1fv4u9DXffypVKrbr3lw6FDWRWd4Zsaw+npu/OUvOVJYyH/GjOFoSQkZYHBDA6UNDfxs4ULOOf98\ntm7dymuTJ9NcUMCrs2axs7QUM0P19ey94w5SO3ZQ+dxzfPrhh6n+2tdYnUgwbdo05s2bx5o1a9i0\naROFhYVUVlYiierqav5ZW8tFF15I+aQRzJ50/FodLRkj9Xaa3fVN7K4/zBupRv6yvprEaSUkJkzk\n/vd8jvtmXMHBgmLsQBO7DjRx3eqfMGfHS9x/3mUsP28edacPZ2Otsbx2EzCKD131Q366+seMf6Wa\noTNn8N2Lr2fl1EtRQhQXJClIioKEKEgmKEiIZOL47URUtDNmmIFhZI4tcxvuz1go8kdbMhxpzvBO\nS4ajzRkyBs8umsPEst51r3WX9/F3k4/q8VE9A3pUz7BhHDp0iPLyctLpNI2NjZSUlBwb1VNaXU3J\nokUk/vWvDt+Pb/3xj5xRVdU2qgcYPXo06XSapqYmEolEGPVz+DDFd96Jvvc9rLKS+iefZPjIkcfe\nr/X19RQXFx93IVbS//XZaGpqIpMsJNWYpuFgmvddOY9hL4frWC2DBvHq7MtZN/czbBv9fg4fbeHw\n0RZobOTGB+/iYy+HmUvXTbqAb1TdRMPgM7qdozckWPv1jzJpVM9nFPWLu865vrFnD2zeHGb1bG6G\n0aNh4sSwGEpPuid+/WuYPx9GjeqzqJ0yg+efD6OAHnus7Ru8l1wS7sueyXPlyrB+QWMjLU+sJj33\nUppbjOZMhpaM0Zyxd223RFNYJCQSCRAiEXXfZP87KJlgUDJBYTJBYUGCQclw1tBbXvidc6eOTCbM\nm/+FL4QvZfWnbdvgnntg2bIwqdv27WGuLwh/3EpKwpoeTzwBN93U9rys6xq54qSP6pFUJalG0jZJ\nt3bweJGk30ePr5d0ZtZjt0X310i6tP1znXN5bvFiuPHGsIZvF9eATroPfCCsQbxrFzz0UFvRNwtT\nQpx9dvijMGtW25nBCy+EqZ57MeIpV3R5xC8pCdQCcwnr724ArjKz17LafAWYYmZflrQAuMLMPiNp\nMrASmA6MAZ4GJplZy4n26Uf8zuWRDRvC/FqpFEyZEqaILi+PN9POnTB16vHrCU+YEBZ5Wr8+TCkx\neDDcfTfccENOHP2f1K4eSRcAS8zs0mj7NgAzuzOrzdqozYuSCoC9wAjg1uy22e1OtE8v/M7lmR07\nwsLtNTXhusGKFaHwZq8RUFfX+fOHDGmbsC2dhrff7rxtWVlboW5oCJO2dSSRCNNCPPJImN8nez6h\nOXPCmsUAH/4wXHttuDZQUQFNTSFDY2PYT+vSjy0t4cthJSVh3eKhQ2HfvnD/W2+FL41VVPR6ZuGe\nFP5oiFHnN+BK4L6s7WuAe9u12QKUZ21vB8qAe4Grs+7/DXBlV/ucNm2aOefyzP79ZrNnm4VOFbOK\niuMfP+20tsfa3+66q63d8uWdtwOzgwfb2s6a1Xm7z362rd2WLSd+zZN1W7Om178+whroXdZ0M+vW\nOP6OzmHanyZ01qY7zw0vIN0A3AAwfvz4bsRyzg0opaVh3d7Fi2H5chg+/PjHy8rCBdiOZA/3LCoK\nbbujpKTzttmLtBcWdtwukwlH61Io3aWlcOhQOOpv7U1pPbvI3i4qCkf4Bw60lX3ot4vb3tXjnHMD\nwMke1bMBOEvSREmFwAJgVbs2q4Drop+vBJ6JTj1WAQuiUT8TgbOAjmcAc8451y+67Ooxs2ZJC4G1\nQBJYZmZbJX2X0Ke0itB3f7+kbcABwh8HonYPAq8BzcBXrYsRPc455/qWf4HLOecGAJ+W2TnnXKe8\n8DvnXJ7xwu+cc3nGC79zzuUZL/zOOZdncnJUj6Q64M1ePr0M2HcS45wsnqtnPFfPeK6eGYi5JpjZ\niK6b5Wjh/39I2tjdIU39yXP1jOfqGc/VM/mey7t6nHMuz3jhd865PDMQC/+v4g7QCc/VM56rZzxX\nz+R1rgHXx++cc+7EBuIRv3POuRMY0IVf0iJJJqmbqzL0LUl3SHpFUrWkpySNiTsTgKQfSfpHlO1R\nScPizgQg6VOStkrKSIp1BIakKkk1krZJujXOLNkkLZOUkrQl7izZJI2T9Kyk16P/w5vjzgQgqVjS\n3yVtjnItjTtTK0lJSZskre7rfQ3Ywi9pHGGB+J1xZ8nyIzObYmaVwGrgO3EHiqwDzjGzKUAtcFvM\neVptAT4JPBdnCElJ4OfAx4HJwFWSJseZKctvgaq4Q3SgGbjFzD4IzAC+miO/syPARWY2FagEqiTN\niDlTq5uB1/tjRwO28AN3A9+kk6Ue42BmjVmbQ8iRbGb2lJk1R5t/A8rjzNPKzF43s5q4cwDTgW1m\ntsPMjgK/Az4RcyYAzOw5whoYOcXM9pjZy9HPbxMK2th4U0G0PO3BaHNQdIv9cyipHJgH3Ncf+xuQ\nhV/SfODfZrY57iztSfq+pF3A58idI/5sXwT+FHeIHDMW2JW1vZscKGKnCklnAucB6+NNEkRdKtVA\nClhnZrmQ66eEA9VMf+ysO4ut5yRJTwPv7eCh24FvAx/r30TBiXKZ2eNmdjtwe7R28UJgcS7kitrc\nTjhFX9EfmbqbKweog/tiP0o8FUg6HXgE+Hq7M97YRKsAVkbXsh6VdI6ZxXaNRNLlQMrMXpI0pz/2\necoWfjO7pKP7JZ0LTAQ2K6xuXw68LGm6me2NK1cHHgDW0E+Fv6tckq4DLgcutn4c49uD31ecdgPj\nsrbLgf/ElOWUIWkQoeivMLM/xJ2nPTNrkPRnwjWSOC+OzwTmS7oMKAbOkLTczK7uqx0OuK4eM3vV\nzEaa2ZlmdibhQ/uh/ij6XZF0VtbmfOAfcWXJJqkK+BYw38wOx50nB20AzpI0UVIhYU3pVTFnymkK\nR12/AV43s7viztNK0ojWUWuSBgOXEPPn0MxuM7PyqF4tAJ7py6IPA7Dw57gfSNoi6RVCV1RODHED\n7gWGAuuioaa/iDsQgKQrJO0GLgDWSFobR47owvdCYC3hIuWDZrY1jiztSVoJvAhUSNot6UtxZ4rM\nBK4BLoreU9XREW3cRgPPRp/BDYQ+/j4fPplr/Ju7zjmXZ/yI3znn8owXfuecyzNe+J1zLs944XfO\nuTzjhd855/KMF37nnMszXvidcy7PeOF3zrk88z9XnzbZF3JAvwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1c13d09d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from common.DensityRatio import kernelDensityRatio\n",
    "plt.plot(newdat,tw,lw=2)            # 真の密度比関数のプロット\n",
    "# データ点のプロット\n",
    "plt.scatter(nu.reshape(n,),np.repeat(0.3,n),marker='.',c='black',s=20)\n",
    "plt.scatter(de.reshape(m,),np.repeat(0.1,m),marker='x',c='gray',s=20)\n",
    "# 以下のカーネル幅で推定\n",
    "gammas = np.array([0.01, 0.1, 1])\n",
    "for g in gammas:\n",
    "    dr = kernelDensityRatio(gamma=g)\n",
    "    dr.fit(de,nu)                   # データへのフィッティング\n",
    "    drp = dr.predict(newdat)        # 密度比の予測値\n",
    "    plt.plot(newdat,drp, c='red',linestyle='dashed',lw=2)  # プロット\n",
    "plt.show()   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 密度比推定のための交差検証法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.spatial import distance        # distanceの計算\n",
    "cvk = 5                                                      # 交差検証法のK\n",
    "n, m = nu.shape[0], de.shape[0]          # データ数\n",
    "# カーネル幅パラメータの候補を生成\n",
    "idx = np.random.choice(n+m,round((n+m)/2))\n",
    "gammas = 1/np.percentile(distance.pdist(np.r_[nu,de][idx,:]),[1,99])**2\n",
    "gammas = np.logspace(np.log10(gammas.min()/100), np.log10(gammas.max()*100),10)\n",
    "# 正則化パラメータ lambda の候補を生成 \n",
    "lams = np.array([(min(n,m))**(-0.9)])\n",
    "# モデルパラメータの候補\n",
    "modelpars = np.array([(x,y) for x in gammas for y in lams])  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.058736149499858691, 0.015848931924611134)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from common.DensityRatio import kernelDensityRatio\n",
    "# それぞれのデータを5グループに分ける\n",
    "inu = np.repeat(np.arange(cvk),np.ceil(n/cvk))\n",
    "inu = inu[np.random.choice(n,n,replace=False)]\n",
    "ide = np.repeat(np.arange(cvk),np.ceil(m/cvk))\n",
    "ide = ide[np.random.choice(m,m,replace=False)]\n",
    "cvloss = []\n",
    "for gamma, lam in modelpars:\n",
    "    tcvloss = []\n",
    "    for k in np.arange(cvk):\n",
    "        # トレーニングデータ\n",
    "        trnu, trde = nu[inu!=k,:], de[ide!=k,:]\n",
    "        # テストデータ\n",
    "        tenu, tede = nu[inu==k,:], de[ide==k,:]\n",
    "        # 指定されたモデルパラメータで密度比を推定\n",
    "        kdr =kernelDensityRatio(gamma=gamma, lam=lam)\n",
    "        kdr.fit(trde, trnu)                             # 推定\n",
    "        wde = kdr.predict(tede)                         # de上で予測\n",
    "        wnu = kdr.predict(tenu)                         # nu上で予測\n",
    "        tcvloss.append(np.mean(wde**2)/2-np.mean(wnu))  # 二乗損失\n",
    "    cvloss.append(np.mean(tcvloss))\n",
    "# 最適なモデルパラメータ\n",
    "optgamma,optlam = modelpars[np.argmin(cvloss),:]\n",
    "optgamma,optlam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEOCAYAAABxdpuaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAFW5JREFUeJzt3X+QXWd93/H3B1mGrQmWXMu2JBAi\nwag4GKywJUMJLYOtikxapDgNkGGonIFqmtTtJJlqsErSmaTJoFb50YQk0yoOsWCmTQgVkpLQKLYC\noR0IwxqZKjajyiFQLGn8g6DwawOy+PaPe9beXd+V7pF077mrfb9mdu455z57znefvXs/e85zzzmp\nKiRJauNZXRcgSVp8DA9JUmuGhySpNcNDktSa4SFJas3wkCS1ZnhIklozPCRJrRkekqTWDA9JUmtX\ndF3ApXbttdfW+vXruy5DkhaV+++//4mqWjVo+07DI8kbgF8FlgF3V9Wuec8/G3gf8ErgS8Cbq+rz\n51rn+vXrmZqaGk7BknSZSvKFNu07O2yVZBnwG8D3AzcBP5LkpnnN3g58uapeDPwK8B9HW6UkqZ8u\nxzxeBTxcVZ+rqm8BvwtsmddmC7C3mf4gcGuSjLBGSVIfXYbHWuCLs+YfaZb1bVNVTwJ/A/zd+StK\nsj3JVJKpxx9/fEjlSpJmdBke/fYg5t9cZJA2VNWeqpqsqslVqwYe75EkXaAuw+MR4AWz5p8PnFyo\nTZIrgKuBvx5JdZKkBXUZHp8CbkzyoiRXAm8BDs5rcxDY1kz/M+BPy1sfSlLnOvuoblU9meRO4BC9\nj+q+t6oeTPJzwFRVHQR+G3h/kofp7XG8pat6JUlP6/Q8j6r6MPDhecv+/azpvwV+eNR1SZLOzcuT\nSJJau+wuTyJJS83+IyfYfegYJ09Ps2bFBDs2b2DrxvlnPlxahockLWL7j5xg576jTJ85C8CJ09Ps\n3HcUYKgB4mErSVrEdh869lRwzJg+c5bdh44NdbuGhyQtYidPT7dafql42EqSLlAXYw3zrVkxwYk+\nQbFmxcRQt+uehyRdgJmxhhOnpymeHmvYf+TESOvYsXkDE8uXzVk2sXwZOzZvGOp2DQ9JugBdjTXM\nt3XjWt59+82sXTFBgLUrJnj37Tf7aStJGkddjTX0s3Xj2pEfLnPPQ5IuwEJjCsMeaxgXhockXYCu\nxhrGhYetJOkCzBwm6vrTVl0xPCTpAnUx1jAuPGwlSWrNPQ9Ji844nJy31BkekhaVri4EqLk8bCVp\nURmXk/OWOsND0qIyTifnLWWGh6RFZamfnDcuDA9Ji8pSPzlvXDhgLmlRWeon540Lw0PSorOUT84b\nF50ctkpyTZJ7kxxvHlcu0O6Pk5xO8oejrlGStLCuxjzuAg5X1Y3A4Wa+n93A20ZWlSRpIF2FxxZg\nbzO9F9jar1FVHQa+OqqiJEmD6So8rq+qUwDN43Ud1SFJugBDGzBPch9wQ5+n3jWEbW0HtgOsW7fu\nUq9ekjTP0MKjqm5b6LkkjyZZXVWnkqwGHrvIbe0B9gBMTk7WxaxLknR+XR22Oghsa6a3AQc6qkOS\ndAG6Co9dwKYkx4FNzTxJJpPcPdMoyf8Cfh+4NckjSTZ3Uq0kaY5OThKsqi8Bt/ZZPgW8Y9b8a0dZ\nlyRpMF7bSpLUmuEhSWrN8JAktWZ4SJJaMzwkSa0ZHpKk1gwPSVJrhockqTXDQ5LUmuEhSWrN8JAk\ntWZ4SJJaMzwkSa0ZHpKk1gwPSVJrhockqTXDQ5LUWid3EpS0OO0/coLdh45x8vQ0a1ZMsGPzBrZu\nXNt1WeqA4SFpIPuPnGDnvqNMnzkLwInT0+zcdxTAAFmCPGwlaSC7Dx17KjhmTJ85y+5DxzqqSF0y\nPCQN5OTp6VbLdXkzPCQNZM2KiVbLdXnrJDySXJPk3iTHm8eVfdrckuQTSR5M8n+SvLmLWiX17Ni8\ngYnly+Ysm1i+jB2bN3RUkbrU1Z7HXcDhqroRONzMz/cN4J9X1XcDbwD+c5IVI6xR0ixbN67l3bff\nzNoVEwRYu2KCd99+s4PlS1RXn7baAryumd4LfBR45+wGVfV/Z02fTPIYsAo4PZoSJc23deNaw0JA\nd3se11fVKYDm8bpzNU7yKuBK4C9HUJsk6TyGtueR5D7ghj5PvavlelYD7we2VdW3F2izHdgOsG7d\nupaVSpLaGlp4VNVtCz2X5NEkq6vqVBMOjy3Q7nnAHwE/XVV/fo5t7QH2AExOTtbFVS5JOp+uDlsd\nBLY109uAA/MbJLkS+BDwvqr6/RHWJkk6j67CYxewKclxYFMzT5LJJHc3bd4E/EPgjiQPNF+3dFOu\nJGm2VF1eR3kmJydramqq6zIkaVFJcn9VTQ7a3jPMJUmtGR6SpNYMD0lSa4aHJKk1w0OS1JrhIUlq\nzfCQJLVmeEiSWjM8JEmtGR6SpNYMD0lSa4aHJKk1w0OS1JrhIUlqzfCQJLVmeEiSWjM8JEmtGR6S\npNYMD0lSa4aHJKm184ZHkquSPKuZfkmSNyZZPvzSJEnjapA9j48Bz0myFjgM/ChwzzCLkiSNt0HC\nI1X1DeB24D1V9YPATRez0STXJLk3yfHmcWWfNi9Mcn+SB5I8mORfXsw2JUmXzkDhkeTVwFuBP2qW\nXXGR270LOFxVN9Lbm7mrT5tTwD+oqluA7wXuSrLmIrcrSboEBgmPnwB2Ah+qqgeTfCfwkYvc7hZg\nbzO9F9g6v0FVfauqvtnMPnvAWiVJI3DePYiq+jPgzwCagfMnqurfXOR2r6+qU836TyW5rl+jJC+g\nt7fzYmBHVZ28yO1Kki6BQT5t9d+SPC/JVcBDwLEkOwb4vvuS/EWfry2DFldVX6yql9MLj21Jrl9g\nW9uTTCWZevzxxwddvSTpAg1yKOimqvoKvUNLHwbWAW873zdV1W1V9bI+XweAR5OsBmgeHzvPuk4C\nDwKvXeD5PVU1WVWTq1atGuBHkiRdjEHCY3lzXsdW4EBVnQHqIrd7ENjWTG8DDsxvkOT5SSaa6ZXA\na4BjF7ldSdIlMEh4/Ffg88BVwMeSvBD4ykVudxewKclxYFMzT5LJJHc3bV4KfDLJZ+iNufxiVR29\nyO1Kki6BVLXfiUhyRVU9OYR6Ltrk5GRNTU11XYYkLSpJ7q+qyUHbDzJgfnWSX54ZkE7yS/T2QiRJ\nS9Qgh63eC3wVeFPz9RXgd4ZZlCRpvA1ypvh3VdUPzZr/2SQPDKsgSdL4G2TPYzrJ983MJHkNMD28\nkiRJ426QPY8fA/YmuRoI8NfAHcMsSpI03ga5PMkDwCuSPK+Zv9iP6UqSFrkFwyPJTy2wHICq+uUh\n1SRJGnPn2vP4jpFVIUlaVBYMj6r62VEWIklaPLxHhiSpNcNDktTaIJcnWTaKQiRJi8cgex5/lWRP\nklsz81ErSdKSNkh4bADuA/4VvSD59dlnnEuSlp7zhkdVTVfVB6rqdmAj8Dyae5pLkpamgQbMk/yj\nJL8JfBp4Dr2r60qSlqjzXp4kyV8BDwAfAHZU1deHXpUkaawNcmHEV3g9K0nSbIMctnpPkhUzM0lW\nJnnvEGuSJI25QcLj5VV1emamqr5Mb+BckrREDRIez0qycmYmyTUMdrhLknSZGiQEfgn4eJIPAkXv\nk1a/MNSqJEljbZCbQb0vyRTwenp3Ery9qh4aemWSpLE10OGnJiwuWWA0h75+D1gPfB54UzOW0q/t\n84DPAh+qqjsvVQ2SpAvX1VV17wIOV9WNwOFmfiH/Ac9ol6Sx0lV4bAH2NtN7ga39GiV5JXA98Ccj\nqkuSNICuwuP6qjoF0DxeN79BkmfRG6zfcb6VJdmeZCrJ1OOPP37Ji5UkzTW0j9wmuQ+4oc9T7xpw\nFT8OfLiqvni+K8FX1R5gD8Dk5GS1qVOS1N7QwqOqblvouSSPJlldVaeSrAYe69Ps1cBrk/w48Fzg\nyiRfq6pzjY9Ikkagq5P9DgLbgF3N44H5DarqrTPTSe4AJg0OSRoPXY157AI2JTkObGrmSTKZ5O6O\napIkDShVl9cQweTkZE1NTXVdhiQtKknur6rJQdt3techSVrEDA9JUmteHVdaJPYfOcHuQ8c4eXqa\nNSsm2LF5A1s3ru26LC1Rhoe0COw/coKd+44yfeYsACdOT7Nz31EAA0Sd8LCVtAjsPnTsqeCYMX3m\nLLsPHeuoIi11hoe0CJw8Pd1quTRshoe0CKxZMdFquTRshoe0COzYvIGJ5cvmLJtYvowdmzd0VJGW\nOgfMpUVgZlDcT1tpXBge0iKxdeNaw0Jjw8NWkqTWDA9JUmuGhySpNcNDktSa4SFJas3wkCS1ZnhI\nklozPCRJrRkekqTWDA9JUmuGhySpNcNDktRaJ+GR5Jok9yY53jyuXKDd2SQPNF8HR12nJKm/rvY8\n7gIOV9WNwOFmvp/pqrql+Xrj6MqTJJ1LV+GxBdjbTO8FtnZUhyTpAnQVHtdX1SmA5vG6Bdo9J8lU\nkj9PYsBI0pgY2s2gktwH3NDnqXe1WM26qjqZ5DuBP01ytKr+ss+2tgPbAdatW3dB9UqSBje08Kiq\n2xZ6LsmjSVZX1akkq4HHFljHyebxc0k+CmwEnhEeVbUH2AMwOTlZl6B8SdI5dHXY6iCwrZneBhyY\n3yDJyiTPbqavBV4DPDSyCiVJC+oqPHYBm5IcBzY18ySZTHJ30+alwFSSzwAfAXZVleEhSWNgaIet\nzqWqvgTc2mf5FPCOZvrjwM0jLk2SNIBOwmMc7T9ygt2HjnHy9DRrVkywY/MGtm5c23VZkjSWDA96\nwbFz31Gmz5wF4MTpaXbuOwpggEhSH17bCth96NhTwTFj+sxZdh861lFFkjTeDA/g5OnpVsslaakz\nPIA1KyZaLZekpc7wAHZs3sDE8mVzlk0sX8aOzRs6qkiSxpsD5jw9KO6nrSRpMIZHY+vGtYaFJA3I\nw1aSpNYMD0lSa4aHJKk1w0OS1JrhIUlqzfCQJLVmeEiSWjM8JEmtGR6SpNYMD0lSa4aHJKk1w0OS\n1JrhIUlqzfCQJLXWSXgkuSbJvUmON48rF2i3LsmfJPlskoeSrB9tpZKkfrra87gLOFxVNwKHm/l+\n3gfsrqqXAq8CHhtRfZKkc+jqZlBbgNc103uBjwLvnN0gyU3AFVV1L0BVfW2E9UlP2X/khHeZlObp\nKjyur6pTAFV1Ksl1fdq8BDidZB/wIuA+4K6qOjvCOkfON6rxsv/ICXbuO8r0md7L7sTpaXbuOwrg\n70VL2tDCI8l9wA19nnrXgKu4AngtsBH4f8DvAXcAv91nW9uB7QDr1q27gGrHg29Uc41DkO4+dOyp\n38eM6TNn2X3o2JL8nUgzhjbmUVW3VdXL+nwdAB5Nshqgeew3lvEIcKSqPldVTwL7ge9ZYFt7qmqy\nqiZXrVo1rB9p6M71RrXUzATpidPTFE8H6f4jJ0Zax8nT062WS0tFVwPmB4FtzfQ24ECfNp8CViaZ\nSYPXAw+NoLbO+Eb1tHEJ0jUrJlotl5aKrsJjF7ApyXFgUzNPkskkdwM0Yxv/Fjic5CgQ4Lc6qnck\nfKN62rgE6Y7NG5hYvmzOsonly9ixecNI65DGTScD5lX1JeDWPsungHfMmr8XePkIS+vUjs0b5ox5\nwNJ9o1qzYoITfYJi1EE6M67R9diLNG66+rSV+vCN6mnjFKRbN65dkr8D6VwMjzHjG1WPQSqNN8ND\nY8sglcaX4aFnGIfzKySNN8NDc3iioqRBeEl2zTEu51dIGm+Gh+YYl/MrJI03w0NzeKKipEEYHprD\nM6olDcIBc83h+RWSBmF46Bk8v0LS+XjYSpLUmuEhSWrN8JAktWZ4SJJaMzwkSa0ZHpKk1gwPSVJr\nhockqTXDQ5LUWqqq6xouqSSPA1/ouo4huxZ4ousixox9Mpf9MZf9MVe//nhhVa0adAWXXXgsBUmm\nqmqy6zrGiX0yl/0xl/0x16XoDw9bSZJaMzwkSa0ZHovTnq4LGEP2yVz2x1z2x1wX3R+OeUiSWnPP\nQ5LUmuEhSWrN8JAktWZ4XGaSbE3yW0kOJPnHXdfThSRXJdnb9MNbu65nHPi6eKbmdXJ/kn/SdS1d\nS/KsJL+Q5D1Jtg3yPYbHGEny3iSPJfmLecvfkORYkoeT3HWudVTV/qr6F8AdwJuHWO5Iteyb24EP\nNv3wxpEXOyJt+uRyfV3MdgF/P+8EPjDaKkenZX9sAdYCZ4BHBlm/4TFe7gHeMHtBkmXAbwDfD9wE\n/EiSm5LcnOQP531dN+tbf7r5vsvFPQzYN8DzgS82zc6OsMZRu4fB+2TG5fa6mO0eBv/7uQ14CHh0\n1EWO0D0M/vrYAHyiqn4K+LFBVn7FJS1VF6WqPpZk/bzFrwIerqrPAST5XWBLVb0beMbudpIAu4D/\nWVWfHm7Fo9Omb+j95/R84AEu43+Q2vRJks9yGb4uZmv5GnkucBW9N9DpJB+uqm+PsNyha9kfXwS+\n1bQZ6B8uw2P8reXp/6Kh98b4vedo/6+B24Crk7y4qv7LMIvr2EJ982vAryf5AeAPuiisQwv1yVJ6\nXczWtz+q6k6AJHcAT1xuwXEOC70+fhV4T5LXAh8bZEWGx/hLn2ULntlZVb9G781zKejbN1X1deBH\nR13MmFioT5bS62K2c/79VNU9oytlLCz0+vgG8PY2K7psd+kvI48AL5g1/3zgZEe1jBv75pnsk7ns\nj7kuWX8YHuPvU8CNSV6U5ErgLcDBjmsaF/bNM9knc9kfc12y/jA8xkiS/w58AtiQ5JEkb6+qJ4E7\ngUPAZ4EPVNWDXdbZBfvmmeyTueyPuYbdH14YUZLUmnsekqTWDA9JUmuGhySpNcNDktSa4SFJas3w\nkCS1ZnhIklozPCRJrXlhRGkASX4GeCu9K5I+AdwP/A2wHbgSeBh4W1V9I8k9wDTw94AX0rtI4zbg\n1cAnq+qOZp1fo3dvhduALwP/DvhPwDrgJ6rqYHNJ7ffTu3w4wJ1V9fHh/rTS+bnnIZ1Hkkngh4CN\n9O5SONk8ta+q/n5VvYLepR5mX5V0JfB64CfpXRb+V4DvBm5OckvT5irgo1X1SuCrwM8Dm4AfBH6u\nafMYsKmqvofeHQCX4pVxNYbc85DO7/uAA1U1DZBk5h4hL0vy88AKejcXOjTre/6gqirJUeDRqjra\nfO+DwHp6N6r6FvDHTfujwDer6kzzPeub5cvp3ZvkFno36XnJcH5EqR3DQzq/fvdAgN5tPrdW1Wea\nmwq9btZz32wevz1remZ+5u/uTD19cbmn2lXVt5PMtPlJerdKfQW9IwV/e8E/hXQJedhKOr//DfzT\nJM9J8lzgB5rl3wGcSrKc3njIMFwNnGrudPc2YNmQtiO14p6HdB5V9akkB4HPAF8ApugNlv8M8Mlm\n2VF6YXKp/SbwP5L8MPAR4OtD2IbUmpdklwaQ5LlV9bUkf4fePZ63V9Wnu65L6op7HtJg9iS5CXgO\nsNfg0FLnnockqTUHzCVJrRkekqTWDA9JUmuGhySpNcNDktSa4SFJau3/A/uzG/RUxnMJAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c2e2438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xscale('log'); plt.xlabel('gamma'); plt.ylabel('cv loss')\n",
    "plt.scatter(gammas, cvloss)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.43549589e+00,   1.44049714e+00,   1.44547056e+00,\n",
       "         1.45041570e+00,   1.45533213e+00,   1.46021942e+00,\n",
       "         1.46507713e+00,   1.46990481e+00,   1.47470205e+00,\n",
       "         1.47946839e+00,   1.48420341e+00,   1.48890668e+00,\n",
       "         1.49357775e+00,   1.49821621e+00,   1.50282162e+00,\n",
       "         1.50739354e+00,   1.51193156e+00,   1.51643524e+00,\n",
       "         1.52090416e+00,   1.52533790e+00,   1.52973602e+00,\n",
       "         1.53409811e+00,   1.53842374e+00,   1.54271250e+00,\n",
       "         1.54696397e+00,   1.55117773e+00,   1.55535336e+00,\n",
       "         1.55949045e+00,   1.56358859e+00,   1.56764738e+00,\n",
       "         1.57166639e+00,   1.57564523e+00,   1.57958348e+00,\n",
       "         1.58348075e+00,   1.58733664e+00,   1.59115074e+00,\n",
       "         1.59492266e+00,   1.59865201e+00,   1.60233839e+00,\n",
       "         1.60598141e+00,   1.60958069e+00,   1.61313583e+00,\n",
       "         1.61664646e+00,   1.62011220e+00,   1.62353266e+00,\n",
       "         1.62690747e+00,   1.63023625e+00,   1.63351865e+00,\n",
       "         1.63675428e+00,   1.63994279e+00,   1.64308380e+00,\n",
       "         1.64617697e+00,   1.64922194e+00,   1.65221834e+00,\n",
       "         1.65516584e+00,   1.65806408e+00,   1.66091272e+00,\n",
       "         1.66371142e+00,   1.66645984e+00,   1.66915765e+00,\n",
       "         1.67180450e+00,   1.67440008e+00,   1.67694407e+00,\n",
       "         1.67943613e+00,   1.68187595e+00,   1.68426322e+00,\n",
       "         1.68659762e+00,   1.68887886e+00,   1.69110661e+00,\n",
       "         1.69328060e+00,   1.69540051e+00,   1.69746606e+00,\n",
       "         1.69947695e+00,   1.70143291e+00,   1.70333366e+00,\n",
       "         1.70517891e+00,   1.70696839e+00,   1.70870184e+00,\n",
       "         1.71037899e+00,   1.71199959e+00,   1.71356336e+00,\n",
       "         1.71507007e+00,   1.71651947e+00,   1.71791130e+00,\n",
       "         1.71924534e+00,   1.72052135e+00,   1.72173910e+00,\n",
       "         1.72289836e+00,   1.72399891e+00,   1.72504054e+00,\n",
       "         1.72602303e+00,   1.72694617e+00,   1.72780977e+00,\n",
       "         1.72861363e+00,   1.72935754e+00,   1.73004133e+00,\n",
       "         1.73066481e+00,   1.73122780e+00,   1.73173013e+00,\n",
       "         1.73217163e+00,   1.73255213e+00,   1.73287147e+00,\n",
       "         1.73312950e+00,   1.73332608e+00,   1.73346105e+00,\n",
       "         1.73353428e+00,   1.73354563e+00,   1.73349498e+00,\n",
       "         1.73338219e+00,   1.73320717e+00,   1.73296978e+00,\n",
       "         1.73266992e+00,   1.73230749e+00,   1.73188239e+00,\n",
       "         1.73139453e+00,   1.73084382e+00,   1.73023019e+00,\n",
       "         1.72955354e+00,   1.72881383e+00,   1.72801097e+00,\n",
       "         1.72714491e+00,   1.72621559e+00,   1.72522297e+00,\n",
       "         1.72416700e+00,   1.72304765e+00,   1.72186488e+00,\n",
       "         1.72061867e+00,   1.71930899e+00,   1.71793583e+00,\n",
       "         1.71649918e+00,   1.71499903e+00,   1.71343539e+00,\n",
       "         1.71180827e+00,   1.71011766e+00,   1.70836361e+00,\n",
       "         1.70654612e+00,   1.70466522e+00,   1.70272096e+00,\n",
       "         1.70071337e+00,   1.69864251e+00,   1.69650841e+00,\n",
       "         1.69431114e+00,   1.69205077e+00,   1.68972736e+00,\n",
       "         1.68734099e+00,   1.68489174e+00,   1.68237970e+00,\n",
       "         1.67980495e+00,   1.67716760e+00,   1.67446775e+00,\n",
       "         1.67170551e+00,   1.66888100e+00,   1.66599433e+00,\n",
       "         1.66304563e+00,   1.66003504e+00,   1.65696269e+00,\n",
       "         1.65382873e+00,   1.65063330e+00,   1.64737657e+00,\n",
       "         1.64405869e+00,   1.64067982e+00,   1.63724015e+00,\n",
       "         1.63373985e+00,   1.63017909e+00,   1.62655808e+00,\n",
       "         1.62287701e+00,   1.61913607e+00,   1.61533547e+00,\n",
       "         1.61147542e+00,   1.60755615e+00,   1.60357787e+00,\n",
       "         1.59954080e+00,   1.59544520e+00,   1.59129128e+00,\n",
       "         1.58707930e+00,   1.58280951e+00,   1.57848217e+00,\n",
       "         1.57409752e+00,   1.56965585e+00,   1.56515742e+00,\n",
       "         1.56060250e+00,   1.55599138e+00,   1.55132435e+00,\n",
       "         1.54660170e+00,   1.54182373e+00,   1.53699074e+00,\n",
       "         1.53210304e+00,   1.52716093e+00,   1.52216475e+00,\n",
       "         1.51711481e+00,   1.51201145e+00,   1.50685499e+00,\n",
       "         1.50164577e+00,   1.49638414e+00,   1.49107045e+00,\n",
       "         1.48570505e+00,   1.48028829e+00,   1.47482054e+00,\n",
       "         1.46930218e+00,   1.46373356e+00,   1.45811507e+00,\n",
       "         1.45244708e+00,   1.44673000e+00,   1.44096420e+00,\n",
       "         1.43515008e+00,   1.42928804e+00,   1.42337849e+00,\n",
       "         1.41742183e+00,   1.41141848e+00,   1.40536886e+00,\n",
       "         1.39927338e+00,   1.39313247e+00,   1.38694657e+00,\n",
       "         1.38071611e+00,   1.37444152e+00,   1.36812324e+00,\n",
       "         1.36176173e+00,   1.35535744e+00,   1.34891081e+00,\n",
       "         1.34242231e+00,   1.33589240e+00,   1.32932154e+00,\n",
       "         1.32271021e+00,   1.31605887e+00,   1.30936800e+00,\n",
       "         1.30263809e+00,   1.29586961e+00,   1.28906305e+00,\n",
       "         1.28221890e+00,   1.27533766e+00,   1.26841982e+00,\n",
       "         1.26146589e+00,   1.25447636e+00,   1.24745174e+00,\n",
       "         1.24039254e+00,   1.23329927e+00,   1.22617245e+00,\n",
       "         1.21901260e+00,   1.21182023e+00,   1.20459587e+00,\n",
       "         1.19734004e+00,   1.19005328e+00,   1.18273612e+00,\n",
       "         1.17538908e+00,   1.16801271e+00,   1.16060754e+00,\n",
       "         1.15317411e+00,   1.14571297e+00,   1.13822466e+00,\n",
       "         1.13070973e+00,   1.12316872e+00,   1.11560220e+00,\n",
       "         1.10801070e+00,   1.10039479e+00,   1.09275503e+00,\n",
       "         1.08509197e+00,   1.07740616e+00,   1.06969819e+00,\n",
       "         1.06196859e+00,   1.05421795e+00,   1.04644683e+00,\n",
       "         1.03865578e+00,   1.03084540e+00,   1.02301623e+00,\n",
       "         1.01516885e+00,   1.00730384e+00,   9.99421770e-01,\n",
       "         9.91523208e-01,   9.83608731e-01,   9.75678915e-01,\n",
       "         9.67734337e-01,   9.59775571e-01,   9.51803197e-01,\n",
       "         9.43817790e-01,   9.35819931e-01,   9.27810196e-01,\n",
       "         9.19789166e-01,   9.11757418e-01,   9.03715533e-01,\n",
       "         8.95664089e-01,   8.87603666e-01,   8.79534844e-01,\n",
       "         8.71458200e-01,   8.63374315e-01,   8.55283767e-01,\n",
       "         8.47187134e-01,   8.39084994e-01,   8.30977926e-01,\n",
       "         8.22866506e-01,   8.14751310e-01,   8.06632915e-01,\n",
       "         7.98511895e-01,   7.90388825e-01,   7.82264278e-01,\n",
       "         7.74138826e-01,   7.66013042e-01,   7.57887496e-01,\n",
       "         7.49762757e-01,   7.41639392e-01,   7.33517970e-01,\n",
       "         7.25399056e-01,   7.17283213e-01,   7.09171004e-01,\n",
       "         7.01062992e-01,   6.92959734e-01,   6.84861790e-01,\n",
       "         6.76769716e-01,   6.68684065e-01,   6.60605391e-01,\n",
       "         6.52534245e-01,   6.44471175e-01,   6.36416727e-01,\n",
       "         6.28371447e-01,   6.20335878e-01,   6.12310558e-01,\n",
       "         6.04296027e-01,   5.96292820e-01,   5.88301470e-01,\n",
       "         5.80322509e-01,   5.72356464e-01,   5.64403861e-01,\n",
       "         5.56465224e-01,   5.48541072e-01,   5.40631924e-01,\n",
       "         5.32738294e-01,   5.24860695e-01,   5.16999635e-01,\n",
       "         5.09155620e-01,   5.01329154e-01,   4.93520737e-01,\n",
       "         4.85730866e-01,   4.77960033e-01,   4.70208731e-01,\n",
       "         4.62477445e-01,   4.54766661e-01,   4.47076857e-01,\n",
       "         4.39408511e-01,   4.31762097e-01,   4.24138084e-01,\n",
       "         4.16536940e-01,   4.08959125e-01,   4.01405101e-01,\n",
       "         3.93875321e-01,   3.86370237e-01,   3.78890297e-01,\n",
       "         3.71435945e-01,   3.64007620e-01,   3.56605759e-01,\n",
       "         3.49230792e-01,   3.41883149e-01,   3.34563253e-01,\n",
       "         3.27271523e-01,   3.20008375e-01,   3.12774220e-01,\n",
       "         3.05569465e-01,   2.98394513e-01,   2.91249761e-01,\n",
       "         2.84135605e-01,   2.77052433e-01,   2.70000631e-01,\n",
       "         2.62980579e-01,   2.55992654e-01,   2.49037228e-01,\n",
       "         2.42114666e-01,   2.35225332e-01,   2.28369584e-01,\n",
       "         2.21547775e-01,   2.14760253e-01,   2.08007362e-01,\n",
       "         2.01289441e-01,   1.94606824e-01,   1.87959841e-01,\n",
       "         1.81348817e-01,   1.74774070e-01,   1.68235917e-01,\n",
       "         1.61734666e-01,   1.55270624e-01,   1.48844089e-01,\n",
       "         1.42455358e-01,   1.36104721e-01,   1.29792462e-01,\n",
       "         1.23518862e-01,   1.17284196e-01,   1.11088735e-01,\n",
       "         1.04932742e-01,   9.88164779e-02,   9.27401978e-02,\n",
       "         8.67041512e-02,   8.07085826e-02,   7.47537314e-02,\n",
       "         6.88398321e-02,   6.29671135e-02,   5.71357999e-02,\n",
       "         5.13461098e-02,   4.55982569e-02,   3.98924496e-02,\n",
       "         3.42288911e-02,   2.86077794e-02,   2.30293073e-02,\n",
       "         1.74936624e-02,   1.20010269e-02,   6.55157816e-03,\n",
       "         1.14548796e-03,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "         0.00000000e+00,   0.00000000e+00])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 密度比の推定\n",
    "kdr =kernelDensityRatio(gamma=optgamma, lam=optlam)\n",
    "kdr.fit(de,nu)\n",
    "kdr.predict(newdat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 共変量シフトの下での回帰分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 真の回帰関数 \n",
    "def f(x):            \n",
    "    return (x+2)*(x-3)*x\n",
    "# データ設定\n",
    "ntr, mtr, sdtr   = 100, -1.4, 0.7\n",
    "nte, mte, sdte = 100,  0.8, 0.8\n",
    "# トレーニングデータ生成\n",
    "xtr = np.random.normal(loc=mtr, scale=sdtr, size=ntr).reshape(ntr,1)\n",
    "ytr = f(xtr) + np.random.normal(scale=2, size=ntr).reshape(ntr,1)\n",
    "# テストデータ生成\n",
    "xte = np.random.normal(loc=mte, scale=sdte, size=nte).reshape(nte,1)\n",
    "yte = f(xte) + np.random.normal(scale=2, size=nte).reshape(nte,1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.91734913],\n",
       "       [-3.14963706]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# トレーニングデータ点上での密度比を推定\n",
    "kdr = kernelDensityRatio()\n",
    "kdr.fit(xtr,xte)\n",
    "pw = kdr.predict(xtr)\n",
    "# 重み付き最小二乗法で回帰パラメータを推定\n",
    "W = np.sqrt(np.diag(pw))\n",
    "X =  sm.add_constant(xtr)\n",
    "WX = np.dot(W,X); WY = np.dot(W, ytr)\n",
    "estTheta = np.linalg.solve(np.dot(WX.T,WX), np.dot(WX.T,WY))\n",
    "estTheta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 二標本検定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import scale\n",
    "from sklearn import datasets\n",
    "d = datasets.load_breast_cancer()        # データ読込み\n",
    "de = d.data[d.target==0]; de = scale(de) # スケーリング\n",
    "nu = d.data[d.target==1]; nu = scale(nu) # スケーリング"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.025499999999999998"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kdr =kernelDensityRatio() \n",
    "kdr.fit(de,nu)                               # 密度比の推定\n",
    "L1distEst = np.mean(abs(1-kdr.predict(de)))  # L1距離の推定値\n",
    "nperm = 10000                                # 並べ替え検定の繰り返し数\n",
    "nde, nnu = de.shape[0], nu.shape[0]\n",
    "dall = np.r_[de,nu]\n",
    "permL1dist = []\n",
    "for itr in np.arange(nperm):         \n",
    "    idx = np.random.choice(nde+nnu,nde,replace=False)   # データの並べ替え\n",
    "    perm_de = dall[idx,:]\n",
    "    perm_nu = np.delete(dall,idx,0)\n",
    "    pdr =kernelDensityRatio()\n",
    "    pdr.fit(perm_de, perm_nu)            # 並べ替えデータに対する密度比の推定\n",
    "    permL1dist.append(np.mean(abs(1-pdr.predict(perm_de)))) # L1距離の推定\n",
    "np.mean(L1distEst < np.array(permL1dist))          # 並べ替え検定による p値"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}