DPGMM experiment

dpgmm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "%matplotlib inline\n",
    "from dpgmm import DPGMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# パラメータ\n",
    "mu = [[5.0,-5.0], [5.0,5.0], [-5.0,-5.0]]\n",
    "var = 1.0\n",
    "D = len(mu[0])\n",
    "K = len(mu)\n",
    "N = 30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 真のクラスタ割り当て\n",
    "true_z = np.zeros(N, dtype=np.int)\n",
    "for i in range(N):\n",
    "    true_z[i] = i%K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# データ生成\n",
    "X = []\n",
    "for i in range(N):\n",
    "    k = true_z[i]\n",
    "    x = np.random.normal(loc=mu[k], scale=var, size=D)\n",
    "    X.append(x)\n",
    "X = np.array(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<dpgmm.DPGMM instance at 0x1071d27a0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# フィッティング\n",
    "alpha=1.\n",
    "mu0=0.\n",
    "rho0=1.\n",
    "a0=1.\n",
    "b0=1.\n",
    "max_iter=300\n",
    "model = DPGMM(alpha=alpha,mu0=mu0,rho0=rho0,a0=a0,b0=b0,max_iter=max_iter)\n",
    "model.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEDCAYAAADOc0QpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEaJJREFUeJzt3W2MXOV5xvHrWjskalJEpTSKi10BgtjNC8RLhBKlIYkI\n1XTXXtefQiJQlKhCSA0CpYYSLKKtZNSqSE0i9UVClAilUBIRrNh4swm0pOELCLGLw4tfcFFVTIwC\nTaNGqVpkzd0Ps8OOd2e8c+acPS/P/H/SyN6zZ2buY2uv8+xz7nmOI0IAgLRMVF0AAKB4hDsAJIhw\nB4AEEe4AkCDCHQASRLgDQIII9yHZnq26hiJwHPWSynFI6RxLMsdBn/twbEdEuOo68uI46iWV45DS\nOZZUjoOROwAkiHAHgARVNi1jm/kgABjBMNNGG8soZBDm+wEgG3u4ywFMywBAggh3AEgQ4Q4ACSLc\nASBBhDsAJKjSbhkAiWm3pfl5aWFBmpyUWi1pgjFkFQh3AMVot6Xdu6UDB5a3zcxI+/cT8BXgXxxA\nMebnzwx2qfP1/Hw19Yw5wh1AMRYW+m9fXCy3Dkgi3AEUZXKy//bt28utA5IIdwBFabU6c+y9ZmY6\n21G6ShcOY20ZIDHdbpnFxc6InW6ZwtkeauEwwh1ANrQ7VmrYcKcVEsDwqmx35KSSCSN3AMObm5Om\np1dvP3RImppav/elh/4tw47cx+tfBUA+VbU70kOfGeEOYHhVtTvSQ59ZIeFu+zzbD9k+avuI7Y8V\n8boARtRud6ZQ9u3r/NluF/O6VbU70kOfWSFz7rbvk/RERNxj+xxJvxERv1zjOcy5A+thveenq2h3\nZM79LaW1Qto+V9JhSRdlSWvCHVgnVV30XG/00EsqtxXyIkmvS/qW7cskPSPppoj4dQGvDSCrs81P\nNzncJyY69Tf5GEpUxGlvo6RJSX8fEdsl/VrSbSt3sj1rO7qPAt4XQD/MTyevN0ttz/bdp4BpmfdK\nejIiLlj6+hOSbouIPr8Xnlkc0zLAOmB+OmmlTctExGu2X7G9NSKOSbpK0ot5XxfAiCYmOkHO/PRY\nK6pb5sOS7pF0jqSXJX0xIv5rjecwcgeAjFg4DAASxPIDADDGCHcASBDhDgAJItwBIEGEOwAkiHAH\ngARxmz0AHU24jV0TaqwJwh1AM5YsaEKNNcK/CIBm3MauCTXWCOEOoBm3sWtCjTXCtAyAei8T3J1n\nP368//frUGMNsbYMgPrOZ/erq1cdaiwZC4cByKaOt7EbdMvA666TrrmmHjWWrMzb7AFIQR1vYzdo\nnn3r1nrVWUPjdcoD0Cx1vhZQc4Q7gPpqtTrz6r1mZjrbcVbMuQOotzpeC6gQF1QBIEFcUAWAPBq+\njg3hDgAr1bXvP4NmVAkAZUpgHRvCHQBWSmAdG8IdAFZKoL++sHC3vcH2ou1HinpNAKhEAv31hbVC\n2v6KpI9IOjcidgyxP62QAOqrpv31pfa5294s6T5Jd0r6CuEOAOtj2HAv6jT0DUm3SmoX9HoAgBxy\nh7vtHZJ+HhHPrLHfrO3oPvK+LwCMq94stT3bd5+8UyO2/0LSdZJOS3qHpHMlPRwR165VHNMyAJr+\nSdCyVbK2jO1PSdrDnDuAoSTwSdCylT3nDgDZJfBJ0LoqNNwj4sfDjNoBQFISnwStK0buAKqTwCdB\n64pwB1CdBD4JWlfcrANAtWr6SdC64k5MAJAg7sQEYLzQL38Gwh1A89Evv8p4HjWAtNAvvwrhDqD5\n6JdfhXAH0Hz0y69CuANoPvrlV6EVEkAaxqRfnj53AEgQq0ICwBgj3AEgQYQ7ACSIcAeABBHuAJAg\nwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkKHe4295i+3HbR2y/YPumIgoDAIwu98JhtjdJ2hQR\nC7Z/U9Izkv4oIl5c43ksHAYAGZW2cFhEnIqIhaW//0rSEUnn531dAMDoCp1zt32BpO2SnurzvVnb\n0X0U+b4AME56s9T2bN99ipoasf0uSf8q6c6IeHiY4piWAYBsSl3P3fbbJH1P0v3DBDsAYH0VcUHV\nku6T9IuIuDnD8xi5A0BGpd1mz/bvS3pC0nOS2kubb4+IuTWeR7gDQEbcQxUAEsQ9VAFgjBHuAJAg\nwh0AEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLc\nASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAkqJNxtt2wfs33C9m1FvCYAYHSOiHwv\nYG+QdFzS1ZJOSnpa0uci4sU1nhd53xsAxo1tRYTX2q+IkfsVkk5ExMsR8aakByXtKuB1AQAjKiLc\nz5f0Ss/XJ5e2AQAqUkS49/v1YNV8i+1Z29F9FPC+ADCWerPU9my/fTYW8D4nJW3p+XqzpJ+t3Cki\nZiW9VQQBDwCjKWvO/WlJl9i+0PY5kq6RdKCA1wUAjCj3yD0iTtv+sqQfStog6d6IeCF3ZQCAkeVu\nhRz5jWmFBIDMymyFBADUDOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCCHcASBDhDgAJ\nItwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDC\nHQASlCvcbd9l+6jtn9reb/u8ogoDAIwu78j9UUkfjIhLJR2X9NX8JQEA8soV7hHxo4g4vfTlk5I2\n5y8JAJBXkXPuX5L0gwJfDwAwojXD3fZjtp/v89jVs89eSacl3X+W15m1Hd1HMeUDwPjpzVLbs333\niciXs7a/IOkGSVdFxP9kKS7vewPAuLGtiPBa+23M+SYtSX8m6ZNZgh0AsL5yjdxtn5D0dkn/ubTp\nyYi4YcjnMnIHgIyGHbnnnpYZFeEOANkNG+58QhUAEkS4A0CCCHcASBDhDgAJItwBIEGEOwAkiHAH\ngAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASlOseqsiu3Zbm\n56WFBWlyUmq1pAlOsQAKRriXqN2Wdu+WDhxY3jYzI+3fT8ADa2lHW/Mn5rVwakGTmybVurilCfOD\nMwj3UC3R3Jw0Pb16+6FD0tRU+fUATdGOtnZ/Z7cOHFseGc1sndH+z+4fu4DnHqo1tLDQf/viYrl1\nAE0zf2L+jGCXpAPHDmj+xHxFFdUf4V6iycn+27dvL7cOoGkWTvUfGS2eYmQ0COFekHa7M+2yb1/n\nz3Z79T6tVmeOvdfMTGc7gMEmN/UfGW3fxMhokELm3G3vkXSXpN+OiDeGfE4yc+5ZLpR2u2UWFzsj\ndrplgLUx575s2Dn33OFue4ukeyRtk3R5k8N91DZFLpQidXXoVOnWsHhqUds3bR/bbplhw72IVsiv\nS7pV0vcLeK3K5GlTPNuFUsIdTVfkqDnPSWLCE5q6ZEpTl/BDNYxcpz3bM5JejYjDQ+w7azu6jzzv\nKw03x53F/PyZwS51vp4f4mI8F0qRsqI6VboniekHpnXH43do+oFp7f7ObrUj5w/vGOrNUtuz/fZZ\nM9xtP2b7+T6PXZL2SvraMMVExGxEuPvIdCQrdEfZ09PSHXd0/ty9O1/AZ21T7D25tNvSzp1nfp8L\npUhFUZ0qtDMWpzdLI2K23z5rTstExGf6bbf9IUkXSjpsW5I2S1qwfUVEvDZ62Ws72yh71GmQLKPv\nflM4O3dKBw9Khw9zoRRpKapT5WwnCaZaijdy/ETEcxHxnoi4ICIukHRS0uR6B7s0eJT94IOjT9Nk\naVPsd3I5eLAT5nv3dk4wBDtS0bq4pZmtZ/5wzGydUevibL+aFt3O2I625l6a076f7NPcS3NM76zQ\nyLVlBo2yv/3t5b9nXbNlYqKz/zBtilxAxTiZ8IT2f3Z/7k6V7kli5YXZrCcJidbIYTRubZnuXPct\nt0hHj5593/VqRRzU+njwoLRjR/HvB6SiqHbGuZfmNP3A6h/CQ58/lPwUT5Jry3TnunfuXA72bduk\na6/tv/96rdnSaq2+gCpJd9+dv2sHSFm3nXHvlXs1dcnUyKNsliNYW6PCvd9c99Gj0sUX999/vVoR\nJyak669fvf3gweFaJwHkw3IEa2tUuA+a656YKH/Nlmef7b+dFR6B9VfURd6UNeqC6qALqZdf3ulS\nKXPNFj64BFSnqIu8ZSp7CYdGXVCt052M6lQLgHorsruntIXDRpWnW6YuqyrWqRagVx0W+sKyR44/\nop3/tLoLY5TunmTDHcDZ0QNeL+1o6wN/9wEdfWN17/a+T+/T3iv3Znq9JFshAayNNVzqZf7EfN9g\nl9a3u4dwBxJDD3i9DPr/2Pbubeva3UO4A4mhB7xeBv1/3HX1Xes6TUa4A4mhB7xeBv1/rPcyCUlc\nUB319nhAqrglXb0U+f8xNt0yg9ZWv/76zqdICXsAKRmbcB+0QmMvPlwEjI/Ue/zLvEF2pQatN9Mr\n712aADQDPf7LGn+0g9Z4WYkFvYD00eO/rPHh3u/2eP2woBeQPnr8lzV+Wmbl7fEuu6xz04yDB5f3\nWe/lfwHUAz3+yxp/QbUfFvQC6qeMC53jMOc+Nt0yAOqvzNBNvcefcAdQG+N8Q+uisSokgNrgQmf5\ncoe77RttH7P9gu2/KqIoAMNrR1tzL81p30/2ae6lObWjXXVJq3Chs3y5umVsf1rSLkmXRsT/2X5P\nMWUBGEZTLiB2F89aWSeLma2fXHPutr8r6e6IeGyE5zLnDuTUpLns1C90lqWs5QfeJ+kTtu+U9L+S\n9kTE0zlfE8CQzjaXXbdwn/CEpi6Zql1dqVrztGn7MdvP93nsUufk8FuSPirpFknftd33jGJ71nZ0\nH4UeBTCmmMseT71Zanu27z45p2XmJf1lRPx46et/k/TRiHh9mOKYlgHyacqcO4pTSp+77Rsk/U5E\nfM32+yT9s6TfHSa1CXegGMxlj5eywv0cSfdK+rCkN9WZc/+XIZ9LuANARnxCFQASxCdUAWCMEe4A\nkCDCHQASRLgDQIIIdwBIUKW32RvwYVYAQE6VtUI2zVLrZuPPRhxHvaRyHFI6x5LKcTAtAwAJItwB\nIEGE+/D+vOoCCsJx1EsqxyGlcyxJHAdz7gCQIEbuAJAgwh0AEkS4Z2R7z9LdT95ddS2jsn2X7aO2\nf2p7v+3zqq4pC9st28dsn7B9W9X1jML2FtuP2z5i+wXbN1VdUx62N9hetP1I1bWMyvZ5th9a+tk4\nYvtjVdeUB+Gege0tkq6W9B9V15LTo5I+GBGXSjou6asV1zM02xsk/a2kP5T0fkmfs/3+aqsayWlJ\nfxoRv6fObSr/pKHH0XWTpCNVF5HTNyXNR8Q2SZep4cdDuGfzdUm3Smr0VeiI+FFEnF768klJm6us\nJ6MrJJ2IiJcj4k1JD0raVXFNmUXEqYhYWPr7r9QJkvOrrWo0tjdLmpZ0T9W1jMr2uZKulPQPkhQR\nb0bEL6utKh/CfUi2ZyS9GhGHq66lYF+S9IOqi8jgfEmv9Hx9Ug0NxS7bF0jaLumpaisZ2TfUGfS0\nqy4kh4skvS7pW0vTS/fYfmfVReVR6doydWP7MUnv7fOtvZJul/QH5VY0urMdS0R8f2mfvepMD9xf\nZm059ftYeGN/k7L9Lknfk3RzRPx31fVkZXuHpJ9HxDO2P1V1PTlslDQp6caIeMr2NyXdJumOassa\nHeHeIyI+02+77Q9JulDS4aXFzjZLWrB9RUS8VmKJQxt0LF22vyBph6SrGna/w5OStvR8vVnSzyqq\nJRfbb1Mn2O+PiIerrmdEH5c0Y3tK0jsknWv7HyPi2orryuqkpJMR0f3t6SF1wr2x+BDTCGz/u6SP\nRMQbVdcyCtstSX8t6ZMR8XrV9WRhe6M6F4GvkvSqpKclfT4iXqi0sIzcGSXcJ+kXEXFz1fUUYWnk\nvicidlRdyyhsPyHpjyPimO1ZSe+MiFsqLmtkjNzH099IerukR5d+E3kyIm6otqThRMRp21+W9ENJ\nGyTd27RgX/JxSddJes72s0vbbo+IuQprGnc3Srrf9jmSXpb0xYrryYWROwAkiG4ZAEgQ4Q4ACSLc\nASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIL+H69Qd9hOS8gdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1073a9ad0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# プロット\n",
    "colors = ['red', 'green', 'blue', 'yellow', 'black']\n",
    "ci = 0\n",
    "for k in model.n:\n",
    "    if model.n[k]==0: continue\n",
    "    assign = (model.z==k)\n",
    "    plt.scatter(X[assign,0], X[assign,1], c=colors[ci])\n",
    "    ci += 1\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}