Bayesian_Probabilistic_Matrix_Factorization.ipynb

bpmf_python.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.linalg import inv\n",
    "from numpy.random import multivariate_normal\n",
    "from scipy.stats import wishart\n",
    "from scipy.sparse import lil_matrix, coo_matrix\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bpmf_gibbs_sampling(R, U, V, N, M, D, n_sample):\n",
    "    # 初期値(BPMFの論文と同じ)\n",
    "    beta0 = 2\n",
    "    mu0 = 0\n",
    "    nu0 = D\n",
    "    W0 = np.identity(D)\n",
    "    alpha = 2\n",
    "\n",
    "    for t_ in range(n_sample - 1):\n",
    "        # sample lam_u\n",
    "        S_bar = np.sum([np.outer(U[t_, i, :], U[t_, i, :].T) for i in range(N)], axis=0) / N\n",
    "        U_bar = np.sum(U[t_], axis=0) / N\n",
    "        W0_ast = inv(inv(W0) + N * S_bar +\n",
    "                     (beta0 * N / (beta0 + N)) * np.outer(mu0 - U_bar, (mu0 - U_bar).T))\n",
    "        lam_u = wishart.rvs(df=nu0 + N, scale=W0_ast)\n",
    "\n",
    "        # sample mu_u\n",
    "        mu0_ast =(beta0 * mu0 + N * U_bar) / (beta0 + N)\n",
    "        mu_u = multivariate_normal(mu0_ast, inv((beta0 + N) * lam_u))\n",
    "\n",
    "        # sample lam_v\n",
    "        S_bar = np.sum([np.outer(V[t_, :, j], V[t_, :, j].T) for j in range(M)], axis=0) / M\n",
    "        V_bar = np.sum(V[t_], axis=1) / M\n",
    "        W0_ast = inv(inv(W0) + M * S_bar +\n",
    "                     (beta0 * M / (beta0 + M)) * np.outer(mu0 - V_bar, (mu0 - V_bar).T))\n",
    "        lam_v = wishart.rvs(df=nu0 + M, scale=W0_ast)\n",
    "\n",
    "        # sample mu_v\n",
    "        mu0_ast = (beta0 * mu0 + M * V_bar) / (beta0 + M)\n",
    "        mu_v = multivariate_normal(mu0_ast, inv((beta0 + M) * lam_v))\n",
    "\n",
    "        # sample U\n",
    "        for i in range(N):\n",
    "            V_VT_I = np.sum([np.outer(V[t_, :, j], V[t_, :, j].T)\n",
    "                             for j in R.getrow(i).nonzero()[1]], axis=0)\n",
    "            lam_ast_inv = inv(lam_u + alpha * V_VT_I)\n",
    "\n",
    "            V_R_I = np.sum([V[t_, :, j] * r for j, r in zip(R.getrow(i).nonzero()[1],\n",
    "                                                            R.getrow(i).data[0])], axis=0)\n",
    "            mu_ast = lam_ast_inv.dot((alpha * V_R_I + lam_u.dot(mu_u.T)).T)\n",
    "\n",
    "            U[t_ + 1, i, :] = multivariate_normal(mu_ast, lam_ast_inv)\n",
    "\n",
    "        # sample V\n",
    "        for j in range(M):\n",
    "            U_UT_I = np.sum([np.outer(U[t_ + 1, i, :], U[t_ + 1, i, :].T)\n",
    "                             for i in R.getcol(j).nonzero()[0]], axis=0)\n",
    "            lam_ast_inv = inv(lam_v + alpha * U_UT_I)\n",
    "\n",
    "            U_R_I = np.sum([U[t_ + 1, i, :] * r for i, r in zip(R.getcol(j).nonzero()[0],\n",
    "                                                                R.getcol(j).data)], axis=0)\n",
    "            mu_ast = lam_ast_inv.dot((alpha * U_R_I + lam_v.dot(mu_v.T)).T)\n",
    "\n",
    "            V[t_ + 1, :, j] = multivariate_normal(mu_ast, lam_ast_inv)\n",
    "\n",
    "    return U, V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.array([\n",
    "    [0, 0, 7],\n",
    "    [0, 1, 6],\n",
    "    [0, 2, 7],\n",
    "    [0, 3, 4],\n",
    "    [0, 4, 5],\n",
    "    [0, 5, 4],\n",
    "    [1, 0, 6],\n",
    "    [1, 1, 7],\n",
    "    [1, 3, 4],\n",
    "    [1, 4, 3],\n",
    "    [1, 5, 4],\n",
    "    [2, 1, 3],\n",
    "    [2, 2, 3],\n",
    "    [2, 3, 1],\n",
    "    [2, 4, 1],\n",
    "    [3, 0, 1],\n",
    "    [3, 1, 2],\n",
    "    [3, 2, 2],\n",
    "    [3, 3, 3],\n",
    "    [3, 4, 3],\n",
    "    [3, 5, 4],\n",
    "    [4, 0, 1],\n",
    "    [4, 2, 1],\n",
    "    [4, 3, 2],\n",
    "    [4, 4, 3],\n",
    "    [4, 5, 3]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = coo_matrix((r[:, 2], (r[:, 0], r[:, 1])))\n",
    "r = lil_matrix(r)\n",
    "\n",
    "np.random.seed(1)\n",
    "N = 5\n",
    "M = 6\n",
    "D = 3\n",
    "n_sample = 300\n",
    "\n",
    "U = np.zeros((n_sample, N, D))\n",
    "V = np.zeros((n_sample, D, M))\n",
    "U[0, :, :] = np.random.rand(N, D)\n",
    "V[0, :, :] = np.random.rand(D, M)\n",
    "U, V = bpmf_gibbs_sampling(r, U, V, N, M, D, n_sample)\n",
    "\n",
    "burn_in = 100\n",
    "p = np.empty((N, M))\n",
    "for u in range(N):\n",
    "    for i in range(M):\n",
    "        p[u, i] = np.mean(np.sum(U[burn_in:, u, :] * V[burn_in:, :, i], axis=1))\n",
    "\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# movie lens"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"ml-100k/u.data\", header = None, sep = \"\\t\")\n",
    "\n",
    "#ID取得\n",
    "N = len(np.unique(df.loc[:, 0]))\n",
    "M = len(np.unique(df.loc[:, 1]))\n",
    "\n",
    "#idを0から始まるようにする\n",
    "df.iloc[:, 0] = df.iloc[:, 0] - 1\n",
    "df.iloc[:, 1] = df.iloc[:, 1] - 1\n",
    "\n",
    "#df = df.loc[1:1000, :]\n",
    "r = np.array(df) #今回は簡易的な実験のため学習と精度検証用のデータはわけない\n",
    "\n",
    "r = coo_matrix((r[:, 2], (r[:, 0], r[:, 1])))\n",
    "r = lil_matrix(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4599.483886003494"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(1)\n",
    "\n",
    "import time\n",
    "start = time.time()\n",
    "\n",
    "\n",
    "D = 10\n",
    "n_sample = 50\n",
    "\n",
    "U = np.zeros((n_sample, N, D))\n",
    "V = np.zeros((n_sample, D, M))\n",
    "U[0, :, :] = np.random.rand(N, D)\n",
    "V[0, :, :] = np.random.rand(D, M)\n",
    "U, V = bpmf_gibbs_sampling(r, U, V, N, M, D, n_sample)\n",
    "\n",
    "burn_in = 10\n",
    "p = np.empty((N, M))\n",
    "for u in range(N):\n",
    "    for i in range(M):\n",
    "        p[u, i] = np.mean(np.sum(U[burn_in:, u, :] * V[burn_in:, :, i], axis=1))\n",
    "\n",
    "#print(p)\n",
    "\n",
    "time.time() - start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.71530658114895218"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#精度確認\n",
    "new_data = np.array(df)\n",
    "ret = np.zeros(new_data.shape[0])\n",
    "\n",
    "for i in range(len(ret)):\n",
    "    ret[i] = p[new_data[i, :][0], new_data[i, :][1]]\n",
    "    \n",
    "np.sqrt(np.mean(pow((new_data[:, 2] - ret), 2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}