chore: import upstream snapshot with attribution

2026-07-13 13:30:25 +08:00
commit f19b2512d7
562 changed files with 38082 additions and 0 deletions
@@ -0,0 +1,93 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+env/
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*,cover
+.hypothesis/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# IPython Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# dotenv
+.env
+
+# virtualenv
+venv/
+ENV/
+
+# Spyder project settings
+.spyderproject
+
+# Rope project settings
+.ropeproject
+
+
+# VSCode project settings
+.vscode/
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2018 ctgk
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
@@ -0,0 +1,23 @@
+# PRML
+Python codes implementing algorithms described in Bishop's book "Pattern Recognition and Machine Learning"
+
+## Required Packages
+- python 3
+- numpy
+- scipy
+- jupyter (optional: to run jupyter notebooks)
+- matplotlib (optional: to plot results in the notebooks)
+- sklearn (optional: to fetch data)
+
+## Notebooks
+- [ch1. Introduction](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch01_Introduction.ipynb)
+- [ch2. Probability Distributions](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch02_Probability_Distributions.ipynb)
+- [ch3. Linear Models for Regression](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch03_Linear_Models_for_Regression.ipynb)
+- [ch4. Linear Models for Classification](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch04_Linear_Models_for_Classfication.ipynb)
+- [ch5. Neural Networks](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch05_Neural_Networks.ipynb)
+- [ch6. Kernel Methods](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch06_Kernel_Methods.ipynb)
+- [ch7. Sparse Kernel Machines](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch07_Sparse_Kernel_Machines.ipynb)
+- [ch9. Mixture Models and EM](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch09_Mixture_Models_and_EM.ipynb)
+- [ch10. Approximate Inference](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch10_Approximate_Inference.ipynb)
+- [ch11. Sampling Methods](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch11_Sampling_Methods.ipynb)
+- [ch12. Continuous Latent Variables](https://nbviewer.jupyter.org/github/ctgk/PRML/blob/master/notebooks/ch12_Continuous_Latent_Variables.ipynb)
@@ -0,0 +1,286 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 8. Graphical Models"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import itertools\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "from sklearn.datasets import fetch_mldata\n",
+    "from prml import bayesnet as bn\n",
+    "\n",
+    "\n",
+    "np.random.seed(1234)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "b = bn.discrete([0.1, 0.9])\n",
+    "f = bn.discrete([0.1, 0.9])\n",
+    "\n",
+    "g = bn.discrete([[[0.9, 0.8], [0.8, 0.2]], [[0.1, 0.2], [0.2, 0.8]]], b, f)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "b: DiscreteVariable(proba=[0.1 0.9])\n",
+      "f: DiscreteVariable(proba=[0.1 0.9])\n",
+      "g: DiscreteVariable(proba=[0.315 0.685])\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(\"b:\", b)\n",
+    "print(\"f:\", f)\n",
+    "print(\"g:\", g)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "g.observe(0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "b: DiscreteVariable(proba=[0.25714286 0.74285714])\n",
+      "f: DiscreteVariable(proba=[0.25714286 0.74285714])\n",
+      "g: DiscreteVariable(observed=[1. 0.])\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(\"b:\", b)\n",
+    "print(\"f:\", f)\n",
+    "print(\"g:\", g)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "b.observe(0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "b: DiscreteVariable(observed=[1. 0.])\n",
+      "f: DiscreteVariable(proba=[0.11111111 0.88888889])\n",
+      "g: DiscreteVariable(observed=[1. 0.])\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(\"b:\", b)\n",
+    "print(\"f:\", f)\n",
+    "print(\"g:\", g)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 8.3.3 Illustration: Image de-noising"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.image.AxesImage at 0x2ada5ce4c18>"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAC0tJREFUeJzt3UGoZGeZxvH/M1E3MYsOIU0Tk4kjYTYu4tC4UYaehZJx03GRwaxaZtEuJqA7g5sERJBBndkJGWzsgTESiJomDBODOBNXIZ0gpmNPTJCe2KbpJvTCZCWadxb3tFw7996qW1WnTt1+/z8oqurcunXerttPfd853znnS1UhqZ+/mLoASdMw/FJThl9qyvBLTRl+qSnDLzVl+KWmDL/UlOGXmnrfOleWxMMJpZFVVeZ53VItf5L7krya5PUkDy/zXpLWK4se25/kJuBXwKeAi8ALwINV9cs9fseWXxrZOlr+jwOvV9Wvq+r3wPeB40u8n6Q1Wib8dwC/2fb84rDszyQ5meRskrNLrEvSii2zw2+nrsV7uvVV9RjwGNjtlzbJMi3/ReDObc8/BLy5XDmS1mWZ8L8A3JPkw0k+AHwOOLOasiSNbeFuf1X9IclDwDPATcCpqnplZZVJGtXCQ30Lrcxtfml0aznIR9LBZfilpgy/1JThl5oy/FJThl9qyvBLTRl+qSnDLzVl+KWmDL/UlOGXmjL8UlOGX2rK8EtNGX6pKcMvNWX4paYMv9SU4ZeaMvxSU2udolv9jHl16GSui9RqF7b8UlOGX2rK8EtNGX6pKcMvNWX4paYMv9TUUuP8SS4AbwN/BP5QVUdXUZQOjnXO8rzfdXscwN5WcZDP31XVWyt4H0lrZLdfamrZ8Bfw4yQvJjm5ioIkrcey3f5PVNWbSW4Hnk3yv1X13PYXDF8KfjFIGyar2mGT5FHgnar6xh6vmW7vkEYx5Q6/Wbru8Kuquf7hC3f7k9yc5JZrj4FPA+cWfT9J67VMt/8w8MPh2/V9wPeq6r9WUpWk0a2s2z/Xyuz2Hzib3K2fxW7/3hzqk5oy/FJThl9qyvBLTRl+qSnDLzXlpbub2+ShvFlDdbNq3+vnXYcBt7Pll5oy/FJThl9qyvBLTRl+qSnDLzVl+KWmHOe/AWzyWP2m8rLftvxSW4ZfasrwS00Zfqkpwy81Zfilpgy/1JTj/FrKsufcH9R13whs+aWmDL/UlOGXmjL8UlOGX2rK8EtNGX6pqZnhT3IqyZUk57YtuzXJs0leG+4PjVtmb1W1521MSfa8Lfv7y7y3ljNPy/9d4L7rlj0M/KSq7gF+MjyXdIDMDH9VPQdcvW7xceD08Pg0cP+K65I0skW3+Q9X1SWA4f721ZUkaR1GP7Y/yUng5NjrkbQ/i7b8l5McARjur+z2wqp6rKqOVtXRBdclaQSLhv8McGJ4fAJ4ajXlSFqXzHEJ48eBY8BtwGXgEeBHwBPAXcAbwANVdf1OwZ3ey3MsFzDlqakHechtmc/tgP+75yp+ZvhXyfDvzHCPY8prCUxp3vB7hJ/UlOGXmjL8UlOGX2rK8EtNGX6pKS/dvQbrOO1W2i9bfqkpwy81Zfilpgy/1JThl5oy/FJThl9qynH+A8BxfI3Bll9qyvBLTRl+qSnDLzVl+KWmDL/UlOGXmnKcfwWmvPS2tChbfqkpwy81Zfilpgy/1JThl5oy/FJThl9qamb4k5xKciXJuW3LHk3y2yQ/H26fGbfMG1uSPW/SGOZp+b8L3LfD8n+pqnuH23+utixJY5sZ/qp6Dri6hlokrdEy2/wPJfnFsFlwaGUVSVqLRcP/beAjwL3AJeCbu70wyckkZ5OcXXBdkkaQeU5KSXI38HRVfXQ/P9vhtTfkGTDLntjjTr1xjHnC1Sb/zapqruIWavmTHNn29LPAud1eK2kzzTylN8njwDHgtiQXgUeAY0nuBQq4AHxhxBoljWCubv/KVma3f0eb3IXcZF279bOM2u2XdPAZfqkpwy81Zfilpgy/1JThl5ry0t3aWA7ljcuWX2rK8EtNGX6pKcMvNWX4paYMv9SU4Zeacpxfkxn7dHLH8vdmyy81Zfilpgy/1JThl5oy/FJThl9qyvBLTTnOr1F5Tv7msuWXmjL8UlOGX2rK8EtNGX6pKcMvNWX4paZmhj/JnUl+muR8kleSfHFYfmuSZ5O8NtwfGr/cG1NV7Xmbct3L3paRZM+blpNZf6AkR4AjVfVSkluAF4H7gc8DV6vq60keBg5V1ZdnvNe4/5MncpAvSjF27csw4Iupqrk+uJktf1VdqqqXhsdvA+eBO4DjwOnhZafZ+kKQdEDsa5s/yd3Ax4DngcNVdQm2viCA21ddnKTxzH1sf5IPAk8CX6qq383bJUtyEji5WHmSxjJzmx8gyfuBp4Fnqupbw7JXgWNVdWnYL/DfVfXXM95nczcwl+A2/zjc5l/Myrb5s/UX+A5w/lrwB2eAE8PjE8BT+y1S0nTm2dv/SeBnwMvAu8Pir7C13f8EcBfwBvBAVV2d8V6b28wsYZNbz01myz6OeVv+ubr9q2L4tZ3hH8fKuv2SbkyGX2rK8EtNGX6pKcMvNWX4paa8dLeW4nDdwWXLLzVl+KWmDL/UlOGXmjL8UlOGX2rK8EtNOc6/ArPGujf5lF/H6fuy5ZeaMvxSU4ZfasrwS00Zfqkpwy81ZfilphznXwPH0rWJbPmlpgy/1JThl5oy/FJThl9qyvBLTRl+qamZ4U9yZ5KfJjmf5JUkXxyWP5rkt0l+Ptw+M365klYlsy40keQIcKSqXkpyC/AicD/wD8A7VfWNuVeWbO5VLaQbRFXNdVTZzCP8quoScGl4/HaS88Ady5UnaWr72uZPcjfwMeD5YdFDSX6R5FSSQ7v8zskkZ5OcXapSSSs1s9v/pxcmHwT+B/haVf0gyWHgLaCAr7K1afCPM97Dbr80snm7/XOFP8n7gaeBZ6rqWzv8/G7g6ar66Iz3MfzSyOYN/zx7+wN8Bzi/PfjDjsBrPguc22+RkqYzz97+TwI/A14G3h0WfwV4ELiXrW7/BeALw87Bvd7Lll8a2Uq7/ati+KXxrazbL+nGZPilpgy/1JThl5oy/FJThl9qyvBLTRl+qSnDLzVl+KWmDL/UlOGXmjL8UlOGX2pq3VN0vwX837bntw3LNtGm1rapdYG1LWqVtf3lvC9c6/n871l5craqjk5WwB42tbZNrQusbVFT1Wa3X2rK8EtNTR3+xyZe/142tbZNrQusbVGT1DbpNr+k6Uzd8kuayCThT3JfkleTvJ7k4Slq2E2SC0leHmYennSKsWEatCtJzm1bdmuSZ5O8NtzvOE3aRLVtxMzNe8wsPelnt2kzXq+925/kJuBXwKeAi8ALwINV9cu1FrKLJBeAo1U1+Zhwkr8F3gH+/dpsSEn+GbhaVV8fvjgPVdWXN6S2R9nnzM0j1bbbzNKfZ8LPbpUzXq/CFC3/x4HXq+rXVfV74PvA8Qnq2HhV9Rxw9brFx4HTw+PTbP3nWbtdatsIVXWpql4aHr8NXJtZetLPbo+6JjFF+O8AfrPt+UU2a8rvAn6c5MUkJ6cuZgeHr82MNNzfPnE915s5c/M6XTez9MZ8dovMeL1qU4R/p9lENmnI4RNV9TfA3wP/NHRvNZ9vAx9haxq3S8A3pyxmmFn6SeBLVfW7KWvZboe6Jvncpgj/ReDObc8/BLw5QR07qqo3h/srwA/Z2kzZJJevTZI63F+ZuJ4/qarLVfXHqnoX+Dcm/OyGmaWfBP6jqn4wLJ78s9uprqk+tynC/wJwT5IPJ/kA8DngzAR1vEeSm4cdMSS5Gfg0mzf78BngxPD4BPDUhLX8mU2ZuXm3maWZ+LPbtBmvJznIZxjK+FfgJuBUVX1t7UXsIMlfsdXaw9YZj9+bsrYkjwPH2Drr6zLwCPAj4AngLuAN4IGqWvuOt11qO8Y+Z24eqbbdZpZ+ngk/u1XOeL2SejzCT+rJI/ykpgy/1JThl5oy/FJThl9qyvBLTRl+qSnDLzX1//kg+AXbkzT0AAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "mnist = fetch_mldata(\"MNIST original\")\n",
+    "x = mnist.data[0]\n",
+    "binarized_img = (x > 127).astype(np.int).reshape(28, 28)\n",
+    "plt.imshow(binarized_img, cmap=\"gray\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.image.AxesImage at 0x2ada5d84898>"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADMBJREFUeJzt3U+IpHedx/H3d6NeYg4JMbNDTHZcCV5yiNuNF2UZD0oUYeIhwZxGhG0PBtabIZcEFiEs6q4nYXYdHEGjgagZwrJRRI2nkJ6wmOj4J8hsnE0zo4xgchLNdw/9jLST7qrqep6nnuep7/sFTVfX1NTzrafqU8+vnu/z1C8yE0n1/M3QBUgahuGXijL8UlGGXyrK8EtFGX6pKMMvFWX4paIMv1TUm1a5sIhYy8MJNzY2Zv77uXPnBlt+38ueqqGfsz5lZixyu2hzeG9E3A18EbgO+M/MfHTO7dcy/PPWYcRCz0Uvy+972VM19HPWp97DHxHXAb8EPgBcBJ4D7s/Mn834P4Z/xcuf8ou4T0M/Z31aNPxtPvO/B3gpM3+dmX8EvgGcaHF/klaoTfhvBX6z5++LzXV/JSK2ImI7IrZbLEtSx9rs8NtvaPGGsVRmngJOwfoO+6UparPlvwjctufvtwOvtCtH0qq0Cf9zwB0R8Y6IeAvwMeBsN2VJ6tvSw/7M/FNEPAA8zW6r73Rm/rSzyiZk6D3Ds5a/znu126j6uPdq1ec/9ML8zL9yhr+eVbT6JE2Y4ZeKMvxSUYZfKsrwS0UZfqmolZ7PP6S2Lc0+W2J9tuNs5a2fWa+Xzc3Nhe/HLb9UlOGXijL8UlGGXyrK8EtFGX6pqDKtvnktr1We3Xgt23E6jK5eL275paIMv1SU4ZeKMvxSUYZfKsrwS0UZfqmoMn3+eda11z7lU5nb8ivNZ3PLLxVl+KWiDL9UlOGXijL8UlGGXyrK8EtFterzR8QF4FXgz8CfMnPx7w3uWN9921n33/a7AvrsKffdr27z2IZcbxX6+PN0cZDP+zPzdx3cj6QVctgvFdU2/Al8NyLORcRWFwVJWo22w/73ZuYrEXEL8L2I+HlmPrP3Bs2bgm8M0shEVydXRMQjwGuZ+bkZt+ntTA53+A2jz8e2zuutT5m50IpZetgfEddHxA1XLwMfBF5c9v4krVabYf8R4NvNu++bgK9n5n93UpWk3nU27F9oYXOG/W2GeWMeIo65tnmGnM+grTGv1z71PuyXNG2GXyrK8EtFGX6pKMMvFWX4paJG1eqrasiv1x5zK6/PadXXuQ1oq0/STIZfKsrwS0UZfqkowy8VZfilogy/VNTaTNE95dNm2/azx9yrH6spv1664pZfKsrwS0UZfqkowy8VZfilogy/VJThl4pamz7/OuvzvPa2hqxtzOtlCtzyS0UZfqkowy8VZfilogy/VJThl4oy/FJRc/v8EXEa+AhwOTPvbK67CfgmcAy4ANyXmb/vr8xds/q2Yz7/uu2542Pu44/1vsduDK/lRbb8XwHuvua6B4HvZ+YdwPebvyVNyNzwZ+YzwJVrrj4BnGkunwHu6bguST1b9jP/kczcAWh+39JdSZJWofdj+yNiC9jqezmSDmfZLf+liDgK0Py+fNANM/NUZm5m5uaSy5LUg2XDfxY42Vw+CTzZTTmSVmXuFN0R8RhwHLgZuAQ8DHwHeBy4HXgZuDczr90puN99tepZjaE9sozKrb4hn7MxT+Hd53pZdIruueHvUtvwr6u2z8GY3/j61OZNdZ3X+aLh9wg/qSjDLxVl+KWiDL9UlOGXijL8UlGT+uruqfb5K+vzORvz8RFT4JZfKsrwS0UZfqkowy8VZfilogy/VJThl4qaVJ9/qtqezz+ktrW3eWxTXm99mrVeNjcX/8Ist/xSUYZfKsrwS0UZfqkowy8VZfilogy/VNRK+/wbGxtsb28f+O9tesZD9oTbnjduP3t/Uz5fv8/ntKvXg1t+qSjDLxVl+KWiDL9UlOGXijL8UlGGXypq7hTdEXEa+AhwOTPvbK57BPgn4LfNzR7KzP+au7CiU3Tbxx9Gn8cBjPk563KK7q8Ad+9z/b9l5l3Nz9zgSxqXueHPzGeAKyuoRdIKtfnM/0BE/CQiTkfEjZ1VJGkllg3/l4B3AncBO8DnD7phRGxFxHZEHHxQv6SVm7vDDyAijgFPXd3ht+i/7XNbd/jtY8w7j6bMHX6zLbXlj4ije/78KPDiMvcjaThzT+mNiMeA48DNEXEReBg4HhF3AQlcAD7ZY42SerDQsL+zha3psL/tOhzzEHLK2jwvU35Oeh32S5o+wy8VZfilogy/VJThl4oy/FJRo5qi2yPhujfldVr1CL1VccsvFWX4paIMv1SU4ZeKMvxSUYZfKsrwS0WNqs8/1d7rmKeSHvM67Xu9TPWxr6put/xSUYZfKsrwS0UZfqkowy8VZfilogy/VNRK+/wbGxtsbx88a1ebfnnbXvuYe/V96vt8/3U9J7/tehvDMQhu+aWiDL9UlOGXijL8UlGGXyrK8EtFGX6pqLlTdEfEbcBXgb8FXgdOZeYXI+Im4JvAMeACcF9m/n7Ofa1ls3zK56WP+fiFKffxh7ToFN2LhP8ocDQzn4+IG4BzwD3Ax4ErmfloRDwI3JiZn5lzX+N9pbVg+Pth+JezaPjnDvszcyczn28uvwqcB24FTgBnmpudYfcNQdJEHOozf0QcA94NPAscycwd2H2DAG7pujhJ/Vn42P6IeCvwBPDpzPzDosOeiNgCtpYrT1Jf5n7mB4iINwNPAU9n5hea634BHM/MnWa/wA8z811z7me8HzBb8DN/P/zMv5zOPvPH7qP8MnD+avAbZ4GTzeWTwJOHLVLScBbZ2/8+4MfAC+y2+gAeYvdz/+PA7cDLwL2ZeWXOfbXazAz5dcdj3kJO1Zi3nlPWWauvS4Zfexn+fnQ27Je0ngy/VJThl4oy/FJRhl8qyvBLRY1qiu55bA2Nj8/JdLnll4oy/FJRhl8qyvBLRRl+qSjDLxVl+KWiJtXnH1KbfvaYTweu2qef8jf1dMUtv1SU4ZeKMvxSUYZfKsrwS0UZfqkowy8VZZ9/BcY880xVY+7jr+oYBLf8UlGGXyrK8EtFGX6pKMMvFWX4paIMv1TU3PBHxG0R8YOIOB8RP42If26ufyQi/i8i/qf5+XD/5fYnM2f+TFVEzPzR+KzqOYsFDig4ChzNzOcj4gbgHHAPcB/wWmZ+buGFRYw2RX65g9ZFZi70Yp17hF9m7gA7zeVXI+I8cGu78iQN7VCf+SPiGPBu4Nnmqgci4icRcToibjzg/2xFxHZEbLeqVFKn5g77/3LDiLcCPwI+m5nfiogjwO+ABP6F3Y8Gn5hzHw77pZ4tOuxfKPwR8WbgKeDpzPzCPv9+DHgqM++ccz+GX+rZouFfZG9/AF8Gzu8NfrMj8KqPAi8etkhJw1lkb//7gB8DLwCvN1c/BNwP3MXusP8C8Mlm5+Cs+5rsln+WtqMCRx3jM+XnpNNhf1cM/3LLHvMLbV1N+TnpbNgvaT0Zfqkowy8VZfilogy/VJThl4pa6Vd3b2xssL3dzyH+bVsvQ7Zu+mwV2oZczro+rr3c8ktFGX6pKMMvFWX4paIMv1SU4ZeKMvxSUas+pfe3wP/uuepmdr8KbIzGWttY6wJrW1aXtf1dZr5tkRuuNPxvWHjEdmZuDlbADGOtbax1gbUta6jaHPZLRRl+qaihw39q4OXPMtbaxloXWNuyBqlt0M/8koYz9JZf0kAGCX9E3B0Rv4iIlyLiwSFqOEhEXIiIF5qZhwedYqyZBu1yRLy457qbIuJ7EfGr5ve+06QNVNsoZm6eMbP0oOtubDNer3zYHxHXAb8EPgBcBJ4D7s/Mn620kANExAVgMzMH7wlHxD8CrwFfvTobUkT8K3AlMx9t3jhvzMzPjKS2RzjkzM091XbQzNIfZ8B11+WM110YYsv/HuClzPx1Zv4R+AZwYoA6Ri8znwGuXHP1CeBMc/kMuy+elTugtlHIzJ3MfL65/CpwdWbpQdfdjLoGMUT4bwV+s+fvi4xryu8EvhsR5yJia+hi9nHk6sxIze9bBq7nWnNnbl6la2aWHs26W2bG664NEf79vh9pTC2H92bmPwAfAj7VDG+1mC8B72R3Grcd4PNDFtPMLP0E8OnM/MOQtey1T12DrLchwn8RuG3P328HXhmgjn1l5ivN78vAt9n9mDIml65Oktr8vjxwPX+RmZcy88+Z+TrwHwy47pqZpZ8AvpaZ32quHnzd7VfXUOttiPA/B9wREe+IiLcAHwPODlDHG0TE9c2OGCLieuCDjG/24bPAyebySeDJAWv5K2OZufmgmaUZeN2NbcbrQQ7yaVoZ/w5cB5zOzM+uvIh9RMTfs7u1h91vNv76kLVFxGPAcXbP+roEPAx8B3gcuB14Gbg3M1e+4+2A2o5zyJmbe6rtoJmln2XAddfljNed1OMRflJNHuEnFWX4paIMv1SU4ZeKMvxSUYZfKsrwS0UZfqmo/wfVHcvvUKWKAwAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "indices = np.random.choice(binarized_img.size, size=int(binarized_img.size * 0.1), replace=False)\n",
+    "noisy_img = np.copy(binarized_img)\n",
+    "noisy_img.ravel()[indices] = 1 - noisy_img.ravel()[indices]\n",
+    "plt.imshow(noisy_img, cmap=\"gray\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "markov_random_field = np.array([\n",
+    "        [[bn.discrete([0.5, 0.5], name=f\"p(z_({i},{j}))\") for j in range(28)] for i in range(28)], \n",
+    "        [[bn.DiscreteVariable(2) for _ in range(28)] for _ in range(28)]])\n",
+    "a = 0.9\n",
+    "b = 0.9\n",
+    "pa = [[a, 1 - a], [1 - a, a]]\n",
+    "pb = [[b, 1 - b], [1 - b, b]]\n",
+    "for i, j in itertools.product(range(28), range(28)):\n",
+    "    bn.discrete(pb, markov_random_field[0, i, j], out=markov_random_field[1, i, j], name=f\"p(x_({i},{j})|z_({i},{j}))\")\n",
+    "    if i != 27:\n",
+    "        bn.discrete(pa, out=[markov_random_field[0, i, j], markov_random_field[0, i + 1, j]], name=f\"p(z_({i},{j}), z_({i+1},{j}))\")\n",
+    "    if j != 27:\n",
+    "        bn.discrete(pa, out=[markov_random_field[0, i, j], markov_random_field[0, i, j + 1]], name=f\"p(z_({i},{j}), z_({i},{j+1}))\")\n",
+    "    markov_random_field[1, i, j].observe(noisy_img[i, j], proprange=0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.image.AxesImage at 0x2ada61c9f28>"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAC4FJREFUeJzt3UHMZWV9x/Hvr6AbZAGhTCcIHWtINy6wEDeaZrrQUDfgAiOrMV2Mi5LUncQNJI2JadS2KxMaiWNSsSSoTEhTJKYtrggDMTJKUWJGHJnM1EwTYWWUfxfvGfI6vO9777z3nnvOO//vJ7m59z1z55z/e977u89zznPufVJVSOrnj6YuQNI0DL/UlOGXmjL8UlOGX2rK8EtNGX6pKcMvNWX4paau3eTGkng5oTSyqsoyz1up5U9yd5JXkrya5MFV1iVps7Lfa/uTXAP8FPgocBZ4Hri/qn6yx/+x5ZdGtomW/0PAq1X186r6LfAt4J4V1idpg1YJ/y3AL7f9fHZY9geSHE9yKsmpFbYlac1WOeG3U9fiHd36qnoEeATs9ktzskrLfxa4ddvP7wVeX60cSZuySvifB25P8r4k7wY+BZxcT1mSxrbvbn9V/S7JA8DTwDXAo1X147VVJmlU+x7q29fGPOaXRreRi3wkHVyGX2rK8EtNGX6pKcMvNWX4paY2+nn+OVs05JksNXoiHRi2/FJThl9qyvBLTRl+qSnDLzVl+KWmHOobOJSnbmz5paYMv9SU4ZeaMvxSU4ZfasrwS00Zfqkpx/k1W37Mely2/FJThl9qyvBLTRl+qSnDLzVl+KWmDL/U1Erj/EnOAG8Avwd+V1V3raMo9bDqDNFeB7CadVzk81dV9es1rEfSBtntl5paNfwFfC/JC0mOr6MgSZuxarf/w1X1epKbgWeS/E9VPbv9CcObgm8M0sxk1ZMub68oeRh4s6q+tMdz1rMxXRXW9drbTdcTflW11C++725/kuuSXH/pMfAx4PR+1ydps1bp9h8CvjO8u14LfLOq/mMtVUka3dq6/UttzG7/gbPJ18e62e3fm0N9UlOGX2rK8EtNGX6pKcMvNWX4pab86u4DYMyPrh7kobxF9vrdug4DbmfLLzVl+KWmDL/UlOGXmjL8UlOGX2rK8EtNOc5/Fbhax+oXjcWv8nv7td+2/FJbhl9qyvBLTRl+qSnDLzVl+KWmDL/UlOP8OrDGvA6gA1t+qSnDLzVl+KWmDL/UlOGXmjL8UlOGX2pqYfiTPJrkQpLT25bdmOSZJD8b7m8Yt8yrW1XteZtSksluGtcyLf/XgbsvW/Yg8P2quh34/vCzpANkYfir6lng4mWL7wFODI9PAPeuuS5JI9vvMf+hqjoHMNzfvL6SJG3C6Nf2JzkOHB97O5KuzH5b/vNJDgMM9xd2e2JVPVJVd1XVXfvclqQR7Df8J4Fjw+NjwJPrKUfSpmSJrzB+DDgK3AScBx4Cvgs8DtwGvAbcV1WXnxTcaV1+xnIHUw/n7eUgD7mtsl8P+O+9VPELw79Ohn9nq/4NDvILdUxjvrbnvM+XDb9X+ElNGX6pKcMvNWX4paYMv9SU4ZeaMvxSU4ZfasrwS00Zfqkpwy81Zfilpgy/1JThl5pyiu4DYM4fH53SnL8H4SCw5ZeaMvxSU4ZfasrwS00Zfqkpwy81Zfilphzn34Al5kbYUCVXl0X7zesA9mbLLzVl+KWmDL/UlOGXmjL8UlOGX2rK8EtNLQx/kkeTXEhyetuyh5P8KskPh9vHxy1z3qpqz5s0R8u0/F8H7t5h+T9W1R3D7d/XW5aksS0Mf1U9C1zcQC2SNmiVY/4HkvxoOCy4YW0VSdqI/Yb/q8D7gTuAc8CXd3tikuNJTiU5tc9tSRpBljkhleQI8FRVfeBK/m2H516VZ79WPannB3vGMebJ1jn/zapqqeL21fInObztx08Ap3d7rqR5WviR3iSPAUeBm5KcBR4Cjia5AyjgDPCZEWuUNIKluv1r25jd/h3NuQt5kK3ydznIf5NRu/2SDj7DLzVl+KWmDL/UlOGXmjL8UlN+dbdmq+sVeptiyy81Zfilpgy/1JThl5oy/FJThl9qyvBLTTnOr8mM/XFyx/L3ZssvNWX4paYMv9SU4ZeaMvxSU4ZfasrwS005zq9R+Zn8+bLll5oy/FJThl9qyvBLTRl+qSnDLzVl+KWmFo7zJ7kV+AbwJ8BbwCNV9c9JbgT+DTgCnAE+WVX/N16p0xpzvHrRusccz97kFO1XynH8cWWJF95h4HBVvZjkeuAF4F7g08DFqvpikgeBG6rqcwvWNd9X2gJThsTw60pU1VI7bmG3v6rOVdWLw+M3gJeBW4B7gBPD006w9YYg6YC4omP+JEeADwLPAYeq6hxsvUEAN6+7OEnjWfra/iTvAZ4APltVv1m2S5bkOHB8f+VJGsvCY36AJO8CngKerqqvDMteAY5W1bnhvMB/VdWfL1jPfA8wF/CYf/M85t+ftR3zZ+sv8DXg5UvBH5wEjg2PjwFPXmmRkqazzNn+jwA/AF5ia6gP4PNsHfc/DtwGvAbcV1UXF6xrvs3MAnNuIQ8qW/ZxLNvyL9XtXxfDr+0M/zjW1u2XdHUy/FJThl9qyvBLTRl+qSnDLzXlV3drJQ7XHVy2/FJThl9qyvBLTRl+qSnDLzVl+KWmDL/UlOP8S1plPHvOHwd2nL4vW36pKcMvNWX4paYMv9SU4ZeaMvxSU4Zfaspx/g1wLF1zZMsvNWX4paYMv9SU4ZeaMvxSU4ZfasrwS00tDH+SW5P8Z5KXk/w4yd8Nyx9O8qskPxxuHx+/XEnrkkVfNJHkMHC4ql5Mcj3wAnAv8Engzar60tIbS+b7rRbSVaKqlrqqbOEVflV1Djg3PH4jycvALauVJ2lqV3TMn+QI8EHguWHRA0l+lOTRJDfs8n+OJzmV5NRKlUpaq4Xd/refmLwH+G/gC1X17SSHgF8DBfw9W4cGf7NgHXb7pZEt2+1fKvxJ3gU8BTxdVV/Z4d+PAE9V1QcWrMfwSyNbNvzLnO0P8DXg5e3BH04EXvIJ4PSVFilpOsuc7f8I8APgJeCtYfHngfuBO9jq9p8BPjOcHNxrXbNt+ZfYDxuqRFrNWrv962L4pfGtrdsv6epk+KWmDL/UlOGXmjL8UlOGX2pqo+G/8847qapRbqtKsudNutrY8ktNGX6pKcMvNWX4paYMv9SU4ZeaMvxSU5v+SO//Ar/Ytugmtr4KbI7mWttc6wJr26911vanVfXHyzxxo+F/x8aTU1V112QF7GGutc21LrC2/ZqqNrv9UlOGX2pq6vA/MvH29zLX2uZaF1jbfk1S26TH/JKmM3XLL2kik4Q/yd1JXknyapIHp6hhN0nOJHlpmHl40inGhmnQLiQ5vW3ZjUmeSfKz4X7HadImqm0WMzfvMbP0pPtubjNeb7zbn+Qa4KfAR4GzwPPA/VX1k40WsoskZ4C7qmryMeEkfwm8CXzj0mxISf4BuFhVXxzeOG+oqs/NpLaHucKZm0eqbbeZpT/NhPtunTNer8MULf+HgFer6udV9VvgW8A9E9Qxe1X1LHDxssX3ACeGxyfYevFs3C61zUJVnauqF4fHbwCXZpaedN/tUdckpgj/LcAvt/18lnlN+V3A95K8kOT41MXs4NClmZGG+5snrudyC2du3qTLZpaezb7bz4zX6zZF+Hf6Tqw5DTl8uKr+Avhr4G+H7q2W81Xg/WxN43YO+PKUxQwzSz8BfLaqfjNlLdvtUNck+22K8J8Fbt3283uB1yeoY0dV9fpwfwH4DluHKXNy/tIkqcP9hYnreVtVna+q31fVW8C/MOG+G2aWfgL416r69rB48n23U11T7bcpwv88cHuS9yV5N/Ap4OQEdbxDkuuGEzEkuQ74GPObffgkcGx4fAx4csJa/sBcZm7ebWZpJt53c5vxepKLfIahjH8CrgEeraovbLyIHST5M7Zae4BrgW9OWVuSx4CjbH3q6zzwEPBd4HHgNuA14L6q2viJt11qO8oVztw8Um27zSz9HBPuu3XOeL2WerzCT+rJK/ykpgy/1JThl5oy/FJThl9qyvBLTRl+qSnDLzX1/2wNFDzYGLDhAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "for _ in range(10000):\n",
+    "    i, j = np.random.choice(28, 2)\n",
+    "    markov_random_field[1, i, j].send_message(proprange=3)\n",
+    "restored_img = np.zeros_like(noisy_img)\n",
+    "for i, j in itertools.product(range(28), range(28)):\n",
+    "    restored_img[i, j] = np.argmax(markov_random_field[0, i, j].proba)\n",
+    "plt.imshow(restored_img, cmap=\"gray\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
@@ -0,0 +1,24 @@
+from prml import (
+    bayesnet,
+    clustering,
+    dimreduction,
+    kernel,
+    linear,
+    markov,
+    nn,
+    rv,
+    sampling
+)
+
+
+__all__ = [
+    "bayesnet",
+    "clustering",
+    "dimreduction",
+    "kernel",
+    "linear",
+    "markov",
+    "nn",
+    "rv",
+    "sampling"
+]
@@ -0,0 +1,7 @@
+from prml.bayesnet.discrete import discrete, DiscreteVariable
+
+
+__all__ = [
+    "DiscreteVariable",
+    "discrete"
+]
@@ -0,0 +1,242 @@
+import numpy as np
+from prml.bayesnet.probability_function import ProbabilityFunction
+from prml.bayesnet.random_variable import RandomVariable
+
+
+class DiscreteVariable(RandomVariable):
+    """
+    Discrete random variable
+    """
+
+    def __init__(self, n_class:int):
+        """
+        intialize a discrete random variable
+
+        parameters
+        ----------
+        n_class : int
+            number of classes
+
+        Attributes
+        ----------
+        parent : DiscreteProbability, optional
+            parent node this variable came out from
+        message_from : dict
+            dictionary of message from neighbor node and itself
+        child : list of DiscreteProbability
+            probability function this variable is conditioning
+        proba : np.ndarray
+            current estimate
+        """
+        self.n_class = n_class
+        self.parent = []
+        self.message_from = {self: np.ones(n_class)}
+        self.child = []
+        self.is_observed = False
+
+    def __repr__(self):
+        string = f"DiscreteVariable("
+        if self.is_observed:
+            string += f"observed={self.proba})"
+        else:
+            string += f"proba={self.proba})"
+        return string
+
+    def add_parent(self, parent):
+        self.parent.append(parent)
+
+    def add_child(self, child):
+        self.child.append(child)
+        self.message_from[child] = np.ones(self.n_class)
+
+    @property
+    def proba(self):
+        return self.posterior
+
+    def receive_message(self, message, giver, proprange):
+        self.message_from[giver] = message
+        self.summarize_message()
+        self.send_message(proprange, exclude=giver)
+
+    def summarize_message(self):
+        if self.is_observed:
+            self.prior = self.message_from[self]
+            self.likelihood = self.prior
+            self.posterior = self.prior
+            return
+
+        self.prior = np.ones(self.n_class)
+        for func in self.parent:
+            self.prior *= self.message_from[func]
+        self.prior /= np.sum(self.prior, keepdims=True)
+
+        self.likelihood = np.copy(self.message_from[self])
+        for func in self.child:
+            self.likelihood *= self.message_from[func]
+
+        self.posterior = self.prior * self.likelihood
+        self.posterior /= self.posterior.sum()
+
+    def send_message(self, proprange=-1, exclude=None):
+        for func in self.parent:
+            if func is not exclude:
+                func.receive_message(self.likelihood, self, proprange)
+        for func in self.child:
+            if func is not exclude:
+                func.receive_message(self.prior, self, proprange)
+
+    def observe(self, data:int, proprange=-1):
+        """
+        set observed data of this variable
+
+        Parameters
+        ----------
+        data : int
+            observed data of this variable
+            This must be smaller than n_class and must be non-negative
+        propagate : int, optional
+            Range to propagate the observation effect to the other random variable using belief propagation alg.
+            If proprange=1, the effect only propagate to the neighboring random variables.
+            Default is -1, which is infinite range.
+        """
+        assert(0 <= data < self.n_class)
+        self.is_observed = True
+        self.receive_message(np.eye(self.n_class)[data], self, proprange=proprange)
+
+
+class DiscreteProbability(ProbabilityFunction):
+    """
+    Discrete probability function
+    """
+
+    def __init__(self, table, *condition, out=None, name=None):
+        """
+        initialize discrete probability function
+
+        Parameters
+        ----------
+        table : (K, ...) np.ndarray or array-like
+            probability table
+            If a discrete variable A is conditioned with B and C,
+            table[a,b,c] give probability of A=a when B=b and C=c.
+            Thus, the sum along the first axis should equal to 1.
+            If a table is 1 dimensional, the variable is not conditioned.
+        condition : tuple of DiscreteVariable, optional
+            parent node, discrete variable this function is conidtioned by
+            len(condition) should equal to (table.ndim - 1)
+            (Default is (), which means no condition)
+        out : DiscreteVariable or list of DiscreteVariable, optional
+            output of this discrete probability function
+            Default is None which construct a new output instance
+        name : str
+            name of this discrete probability function
+        """
+        self.table = np.asarray(table)
+        self.condition = condition
+        if condition:
+            for var in condition:
+                var.add_child(self)
+        self.message_from = {var: var.prior for var in condition}
+
+        if out is None:
+            self.out = [DiscreteVariable(len(table))]
+        elif isinstance(out, DiscreteVariable):
+            self.out = [out]
+        else:
+            self.out = out
+
+        for i, random_variable in enumerate(self.out):
+            random_variable.add_parent(self)
+            self.message_from[random_variable] = np.ones(np.size(self.table, i))
+
+        for random_variable in self.out:
+            self.send_message_to(random_variable, proprange=0)
+
+        self.name = name
+
+    def __repr__(self):
+        if self.name is not None:
+            return self.name
+        else:
+            return super().__repr__()
+
+    def receive_message(self, message, giver, proprange):
+        self.message_from[giver] = message
+        if proprange:
+            self.send_message(proprange, exclude=giver)
+
+    @staticmethod
+    def expand_dims(x, ndim, axis):
+        shape = [-1 if i == axis else 1 for i in range(ndim)]
+        return x.reshape(*shape)
+
+    def compute_message_to(self, destination):
+        proba = np.copy(self.table)
+        for i, random_variable in enumerate(self.out):
+            if random_variable is destination:
+                index = i
+                continue
+            message = self.message_from[random_variable]
+            proba *= self.expand_dims(message, proba.ndim, i)
+        for i, random_variable in enumerate(self.condition, len(self.out)):
+            if random_variable is destination:
+                index = i
+                continue
+            message = self.message_from[random_variable]
+            proba *= self.expand_dims(message, proba.ndim, i)
+        axis = list(range(proba.ndim))
+        axis.remove(index)
+        message = np.sum(proba, axis=tuple(axis))
+        message /= np.sum(message, keepdims=True)
+        return message
+
+    def send_message_to(self, destination, proprange=-1):
+        message = self.compute_message_to(destination)
+        destination.receive_message(message, self, proprange)
+
+    def send_message(self, proprange, exclude=None):
+        proprange = proprange - 1
+
+        for random_variable in self.out:
+            if random_variable is not exclude:
+                self.send_message_to(random_variable, proprange)
+
+        if proprange == 0: return
+
+        for random_variable in self.condition:
+            if random_variable is not exclude:
+                self.send_message_to(random_variable, proprange - 1)
+
+
+def discrete(table, *condition, out=None, name=None):
+    """
+    discrete probability function
+
+    Parameters
+    ----------
+    table : (K, ...) np.ndarray or array-like
+        probability table
+        If a discrete variable A is conditioned with B and C,
+        table[a,b,c] give probability of A=a when B=b and C=c.
+        Thus, the sum along the first axis should equal to 1.
+        If a table is 1 dimensional, the variable is not conditioned.
+    condition : tuple of DiscreteVariable, optional
+        parent node, discrete variable this function is conidtioned by
+        len(condition) should equal to (table.ndim - 1)
+        (Default is (), which means no condition)
+    out : DiscreteVariable, optional
+        output of this discrete probability function
+        Default is None which construct a new output instance
+    name : str
+        name of the discrete probability function
+
+    Returns
+    -------
+    DiscreteVariable
+        output discrete random variable of discrete probability function
+    """
+    function = DiscreteProbability(table, *condition, out=out, name=name)
+    if len(function.out) == 1:
+        return function.out[0]
+    else:
+        return function.out
@@ -0,0 +1,2 @@
+class ProbabilityFunction(object):
+    pass
@@ -0,0 +1,4 @@
+class RandomVariable(object):
+    """
+    Base class for random variable
+    """
@@ -0,0 +1,6 @@
+from .k_means import KMeans
+
+
+__all__ = [
+    "KMeans"
+]
@@ -0,0 +1,53 @@
+import numpy as np
+from scipy.spatial.distance import cdist
+
+
+class KMeans(object):
+
+    def __init__(self, n_clusters):
+        self.n_clusters = n_clusters
+
+    def fit(self, X, iter_max=100):
+        """
+        perform k-means algorithm
+
+        Parameters
+        ----------
+        X : (sample_size, n_features) ndarray
+            input data
+        iter_max : int
+            maximum number of iterations
+
+        Returns
+        -------
+        centers : (n_clusters, n_features) ndarray
+            center of each cluster
+        """
+        I = np.eye(self.n_clusters)
+        centers = X[np.random.choice(len(X), self.n_clusters, replace=False)]
+        for _ in range(iter_max):
+            prev_centers = np.copy(centers)
+            D = cdist(X, centers)
+            cluster_index = np.argmin(D, axis=1)
+            cluster_index = I[cluster_index]
+            centers = np.sum(X[:, None, :] * cluster_index[:, :, None], axis=0) / np.sum(cluster_index, axis=0)[:, None]
+            if np.allclose(prev_centers, centers):
+                break
+        self.centers = centers
+
+    def predict(self, X):
+        """
+        calculate closest cluster center index
+
+        Parameters
+        ----------
+        X : (sample_size, n_features) ndarray
+            input data
+
+        Returns
+        -------
+        index : (sample_size,) ndarray
+            indicates which cluster they belong
+        """
+        D = cdist(X, self.centers)
+        return np.argmin(D, axis=1)
@@ -0,0 +1,10 @@
+from prml.dimreduction.autoencoder import Autoencoder
+from prml.dimreduction.bayesian_pca import BayesianPCA
+from prml.dimreduction.pca import PCA
+
+
+__all__ = [
+    "Autoencoder",
+    "BayesianPCA",
+    "PCA",
+]
@@ -0,0 +1,38 @@
+import numpy as np
+from prml import nn
+
+
+class Autoencoder(nn.Network):
+
+    def __init__(self, *args):
+        self.n_unit = len(args)
+        super().__init__()
+        for i in range(self.n_unit - 1):
+            self.parameter[f"w_encode{i}"] = nn.Parameter(np.random.randn(args[i], args[i + 1]))
+            self.parameter[f"b_encode{i}"] = nn.Parameter(np.zeros(args[i + 1]))
+            self.parameter[f"w_decode{i}"] = nn.Parameter(np.random.randn(args[i + 1], args[i]))
+            self.parameter[f"b_decode{i}"] = nn.Parameter(np.zeros(args[i]))
+
+    def transform(self, x):
+        h = x
+        for i in range(self.n_unit - 1):
+            h = nn.tanh(h @ self.parameter[f"w_encode{i}"] + self.parameter[f"b_encode{i}"])
+        return h.value
+
+    def forward(self, x):
+        h = x
+        for i in range(self.n_unit - 1):
+            h = nn.tanh(h @ self.parameter[f"w_encode{i}"] + self.parameter[f"b_encode{i}"])
+        for i in range(self.n_unit - 2, 0, -1):
+            h = nn.tanh(h @ self.parameter[f"w_decode{i}"] + self.parameter[f"b_decode{i}"])
+        x_ = h @ self.parameter["w_decode0"] + self.parameter["b_decode0"]
+        self.px = nn.random.Gaussian(x_, 1., data=x)
+
+    def fit(self, x, n_iter=100, learning_rate=1e-3):
+        optimizer = nn.optimizer.Adam(self.parameter, learning_rate)
+        for _ in range(n_iter):
+            self.clear()
+            self.forward(x)
+            log_likelihood = self.log_pdf()
+            log_likelihood.backward()
+            optimizer.update()
@@ -0,0 +1,59 @@
+import numpy as np
+from prml.dimreduction.pca import PCA
+
+
+class BayesianPCA(PCA):
+
+    def fit(self, X, iter_max=100, initial="random"):
+        """
+        empirical bayes estimation of pca parameters
+
+        Parameters
+        ----------
+        X : (sample_size, n_features) ndarray
+            input data
+        iter_max : int
+            maximum number of em steps
+
+        Returns
+        -------
+        mean : (n_features,) ndarray
+            sample mean fo the input data
+        W : (n_features, n_components) ndarray
+            projection matrix
+        var : float
+            variance of observation noise
+        """
+        initial_list = ["random", "eigen"]
+        self.mean = np.mean(X, axis=0)
+        self.I = np.eye(self.n_components)
+        if initial not in initial_list:
+            print("availabel initializations are {}".format(initial_list))
+        if initial == "random":
+            self.W = np.eye(np.size(X, 1), self.n_components)
+            self.var = 1.
+        elif initial == "eigen":
+            self.eigen(X)
+        self.alpha = len(self.mean) / np.sum(self.W ** 2, axis=0).clip(min=1e-10)
+        for i in range(iter_max):
+            W = np.copy(self.W)
+            stats = self._expectation(X - self.mean)
+            self._maximization(X - self.mean, *stats)
+            self.alpha = len(self.mean) / np.sum(self.W ** 2, axis=0).clip(min=1e-10)
+            if np.allclose(W, self.W):
+                break
+        self.n_iter = i + 1
+
+    def _maximization(self, X, Ez, Ezz):
+        self.W = X.T @ Ez @ np.linalg.inv(np.sum(Ezz, axis=0) + self.var * np.diag(self.alpha))
+        self.var = np.mean(
+            np.mean(X ** 2, axis=-1)
+            - 2 * np.mean(Ez @ self.W.T * X, axis=-1)
+            + np.trace((Ezz @ self.W.T @ self.W).T) / len(self.mean))
+
+    def maximize(self, D, Ez, Ezz):
+        self.W = D.T.dot(Ez).dot(np.linalg.inv(np.sum(Ezz, axis=0) + self.var * np.diag(self.alpha)))
+        self.var = np.mean(
+            np.mean(D ** 2, axis=-1)
+            - 2 * np.mean(Ez.dot(self.W.T) * D, axis=-1)
+            + np.trace(Ezz.dot(self.W.T).dot(self.W).T) / self.ndim)
@@ -0,0 +1,156 @@
+import numpy as np
+
+
+class PCA(object):
+
+    def __init__(self, n_components):
+        """
+        construct principal component analysis
+
+        Parameters
+        ----------
+        n_components : int
+            number of components
+        """
+        assert isinstance(n_components, int)
+        self.n_components = n_components
+
+    def fit(self, X, method="eigen", iter_max=100):
+        """
+        maximum likelihood estimate of pca parameters
+        x ~ \int_z N(x|Wz+mu,sigma^2)N(z|0,I)dz
+
+        Parameters
+        ----------
+        X : (sample_size, n_features) ndarray
+            input data
+        method : str
+            method to estimate the parameters
+            ["eigen", "em"]
+        iter_max : int
+            maximum number of iterations for em algorithm
+
+        Attributes
+        ----------
+        mean : (n_features,) ndarray
+            sample mean of the data
+        W : (n_features, n_components) ndarray
+            projection matrix
+        var : float
+            variance of observation noise
+        C : (n_features, n_features) ndarray
+            variance of the marginal dist N(x|mean,C)
+        Cinv : (n_features, n_features) ndarray
+            precision of the marginal dist N(x|mean, C)
+        """
+        method_list = ["eigen", "em"]
+        if method not in method_list:
+            print("availabel methods are {}".format(method_list))
+        self.mean = np.mean(X, axis=0)
+        getattr(self, method)(X - self.mean, iter_max)
+
+    def eigen(self, X, *arg):
+        sample_size, n_features = X.shape
+        if sample_size >= n_features:
+            cov = np.cov(X, rowvar=False)
+            values, vectors = np.linalg.eigh(cov)
+            index = n_features - self.n_components
+        else:
+            cov = np.cov(X)
+            values, vectors = np.linalg.eigh(cov)
+            vectors = (X.T @ vectors) / np.sqrt(sample_size * values)
+            index = sample_size - self.n_components
+        self.I = np.eye(self.n_components)
+        if index == 0:
+            self.var = 0
+        else:
+            self.var = np.mean(values[:index])
+
+        self.W = vectors[:, index:].dot(np.sqrt(np.diag(values[index:]) - self.var * self.I))
+        self.__M = self.W.T @ self.W + self.var * self.I
+        self.C = self.W @ self.W.T + self.var * np.eye(n_features)
+        if index == 0:
+            self.Cinv = np.linalg.inv(self.C)
+        else:
+            self.Cinv = np.eye(n_features) / np.sqrt(self.var) - self.W @ np.linalg.inv(self.__M) @ self.W.T / self.var
+
+    def em(self, X, iter_max):
+        self.I = np.eye(self.n_components)
+        self.W = np.eye(np.size(X, 1), self.n_components)
+        self.var = 1.
+        for i in range(iter_max):
+            W = np.copy(self.W)
+            stats = self._expectation(X)
+            self._maximization(X, *stats)
+            if np.allclose(W, self.W):
+                break
+        self.C = self.W @ self.W.T + self.var * np.eye(np.size(X, 1))
+        self.Cinv = np.linalg.inv(self.C)
+
+    def _expectation(self, X):
+        self.__M = self.W.T @ self.W + self.var * self.I
+        Minv = np.linalg.inv(self.__M)
+        Ez = X @ self.W @ Minv
+        Ezz = self.var * Minv + Ez[:, :, None] * Ez[:, None, :]
+        return Ez, Ezz
+
+    def _maximization(self, X, Ez, Ezz):
+        self.W = X.T @ Ez @ np.linalg.inv(np.sum(Ezz, axis=0))
+        self.var = np.mean(
+            np.mean(X ** 2, axis=1)
+            - 2 * np.mean(Ez @ self.W.T * X, axis=1)
+            + np.trace((Ezz @ self.W.T @ self.W).T) / np.size(X, 1))
+
+    def transform(self, X):
+        """
+        project input data into latent space
+        p(Z|X) = N(Z|(X-mu)WMinv, sigma^-2M)
+
+        Parameters
+        ----------
+        X : (sample_size, n_features) ndarray
+            input data
+
+        Returns
+        -------
+        Z : (sample_size, n_components) ndarray
+            projected input data
+        """
+        return np.linalg.solve(self.__M, ((X - self.mean) @ self.W).T).T
+
+    def fit_transform(self, X, method="eigen"):
+        """
+        perform pca and whiten the input data
+
+        Parameters
+        ----------
+        X : (sample_size, n_features) ndarray
+            input data
+
+        Returns
+        -------
+        Z : (sample_size, n_components) ndarray
+            projected input data
+        """
+        self.fit(X, method)
+        return self.transform(X)
+
+    def proba(self, X):
+        """
+        the marginal distribution of the observed variable
+
+        Parameters
+        ----------
+        X : (sample_size, n_features) ndarray
+            input data
+
+        Returns
+        -------
+        p : (sample_size,) ndarray
+            value of the marginal distribution
+        """
+        d = X - self.mean
+        return (
+            np.exp(-0.5 * np.sum(d @ self.Cinv * d, axis=-1))
+            / np.sqrt(np.linalg.det(self.C))
+            / np.power(2 * np.pi, 0.5 * np.size(X, 1)))
@@ -0,0 +1,19 @@
+from prml.kernel.polynomial import PolynomialKernel
+from prml.kernel.rbf import RBF
+
+from prml.kernel.gaussian_process_classifier import GaussianProcessClassifier
+from prml.kernel.gaussian_process_regressor import GaussianProcessRegressor
+from prml.kernel.relevance_vector_classifier import RelevanceVectorClassifier
+from prml.kernel.relevance_vector_regressor import RelevanceVectorRegressor
+from prml.kernel.support_vector_classifier import SupportVectorClassifier
+
+
+__all__ = [
+    "PolynomialKernel",
+    "RBF",
+    "GaussianProcessClassifier",
+    "GaussianProcessRegressor",
+    "RelevanceVectorClassifier",
+    "RelevanceVectorRegressor",
+    "SupportVectorClassifier"
+]
@@ -0,0 +1,37 @@
+import numpy as np
+
+
+class GaussianProcessClassifier(object):
+
+    def __init__(self, kernel, noise_level=1e-4):
+        """
+        construct gaussian process classifier
+
+        Parameters
+        ----------
+        kernel
+            kernel function to be used to compute Gram matrix
+        noise_level : float
+            parameter to ensure the matrix to be positive
+        """
+        self.kernel = kernel
+        self.noise_level = noise_level
+
+    def _sigmoid(self, a):
+        return np.tanh(a * 0.5) * 0.5 + 0.5
+
+    def fit(self, X, t):
+        if X.ndim == 1:
+            X = X[:, None]
+        self.X = X
+        self.t = t
+        Gram = self.kernel(X, X)
+        self.covariance = Gram + np.eye(len(Gram)) * self.noise_level
+        self.precision = np.linalg.inv(self.covariance)
+
+    def predict(self, X):
+        if X.ndim == 1:
+            X = X[:, None]
+        K = self.kernel(X, self.X)
+        a_mean = K @ self.precision @ self.t
+        return self._sigmoid(a_mean)
@@ -0,0 +1,105 @@
+import numpy as np
+
+
+class GaussianProcessRegressor(object):
+
+    def __init__(self, kernel, beta=1.):
+        """
+        construct gaussian process regressor
+
+        Parameters
+        ----------
+        kernel
+            kernel function
+        beta : float
+            precision parameter of observation noise
+        """
+        self.kernel = kernel
+        self.beta = beta
+
+    def fit(self, X, t, iter_max=0, learning_rate=0.1):
+        """
+        maximum likelihood estimation of parameters in kernel function
+
+        Parameters
+        ----------
+        X : ndarray (sample_size, n_features)
+            input
+        t : ndarray (sample_size,)
+            corresponding target
+        iter_max : int
+            maximum number of iterations updating hyperparameters
+        learning_rate : float
+            updation coefficient
+
+        Attributes
+        ----------
+        covariance : ndarray (sample_size, sample_size)
+            variance covariance matrix of gaussian process
+        precision : ndarray (sample_size, sample_size)
+            precision matrix of gaussian process
+
+        Returns
+        -------
+        log_likelihood_list : list
+            list of log likelihood value at each iteration
+        """
+        if X.ndim == 1:
+            X = X[:, None]
+        log_likelihood_list = [-np.Inf]
+        self.X = X
+        self.t = t
+        I = np.eye(len(X))
+        Gram = self.kernel(X, X)
+        self.covariance = Gram + I / self.beta
+        self.precision = np.linalg.inv(self.covariance)
+        for i in range(iter_max):
+            gradients = self.kernel.derivatives(X, X)
+            updates = np.array(
+                [-np.trace(self.precision.dot(grad)) + t.dot(self.precision.dot(grad).dot(self.precision).dot(t)) for grad in gradients])
+            for j in range(iter_max):
+                self.kernel.update_parameters(learning_rate * updates)
+                Gram = self.kernel(X, X)
+                self.covariance = Gram + I / self.beta
+                self.precision = np.linalg.inv(self.covariance)
+                log_like = self.log_likelihood()
+                if log_like > log_likelihood_list[-1]:
+                    log_likelihood_list.append(log_like)
+                    break
+                else:
+                    self.kernel.update_parameters(-learning_rate * updates)
+                    learning_rate *= 0.9
+        log_likelihood_list.pop(0)
+        return log_likelihood_list
+
+    def log_likelihood(self):
+        return -0.5 * (
+            np.linalg.slogdet(self.covariance)[1]
+            + self.t @ self.precision @ self.t
+            + len(self.t) * np.log(2 * np.pi))
+
+    def predict(self, X, with_error=False):
+        """
+        mean of the gaussian process
+
+        Parameters
+        ----------
+        X : ndarray (sample_size, n_features)
+            input
+
+        Returns
+        -------
+        mean : ndarray (sample_size,)
+            predictions of corresponding inputs
+        """
+        if X.ndim == 1:
+            X = X[:, None]
+        K = self.kernel(X, self.X)
+        mean = K @ self.precision @ self.t
+        if with_error:
+            var = (
+                self.kernel(X, X, False)
+                + 1 / self.beta
+                - np.sum(K @ self.precision * K, axis=1))
+            return mean.ravel(), np.sqrt(var.ravel())
+        return mean
@@ -0,0 +1,28 @@
+import numpy as np
+
+
+class Kernel(object):
+    """
+    Base class for kernel function
+    """
+
+    def _pairwise(self, x, y):
+        """
+        all pairs of x and y
+
+        Parameters
+        ----------
+        x : (sample_size, n_features)
+            input
+        y : (sample_size, n_features)
+            another input
+
+        Returns
+        -------
+        output : tuple
+            two array with shape (sample_size, sample_size, n_features)
+        """
+        return (
+            np.tile(x, (len(y), 1, 1)).transpose(1, 0, 2),
+            np.tile(y, (len(x), 1, 1))
+        )
@@ -0,0 +1,43 @@
+import numpy as np
+from prml.kernel.kernel import Kernel
+
+
+class PolynomialKernel(Kernel):
+    """
+    Polynomial kernel
+    k(x,y) = (x @ y + c)^M
+    """
+
+    def __init__(self, degree=2, const=0.):
+        """
+        construct Polynomial kernel
+
+        Parameters
+        ----------
+        const : float
+            a constant to be added
+        degree : int
+            degree of polynomial order
+        """
+        self.const = const
+        self.degree = degree
+
+    def __call__(self, x, y, pairwise=True):
+        """
+        calculate pairwise polynomial kernel
+
+        Parameters
+        ----------
+        x : (..., ndim) ndarray
+            input
+        y : (..., ndim) ndarray
+            another input with the same shape
+
+        Returns
+        -------
+        output : ndarray
+            polynomial kernel
+        """
+        if pairwise:
+            x, y = self._pairwise(x, y)
+        return (np.sum(x * y, axis=-1) + self.const) ** self.degree
@@ -0,0 +1,58 @@
+import numpy as np
+from prml.kernel.kernel import Kernel
+
+
+class RBF(Kernel):
+
+    def __init__(self, params):
+        """
+        construct Radial basis kernel function
+
+        Parameters
+        ----------
+        params : (ndim + 1,) ndarray
+            parameters of radial basis function
+
+        Attributes
+        ----------
+        ndim : int
+            dimension of expected input data
+        """
+        assert params.ndim == 1
+        self.params = params
+        self.ndim = len(params) - 1
+
+    def __call__(self, x, y, pairwise=True):
+        """
+        calculate radial basis function
+        k(x, y) = c0 * exp(-0.5 * c1 * (x1 - y1) ** 2 ...)
+
+        Parameters
+        ----------
+        x : ndarray [..., ndim]
+            input of this kernel function
+        y : ndarray [..., ndim]
+            another input
+
+        Returns
+        -------
+        output : ndarray
+            output of this radial basis function
+        """
+        assert x.shape[-1] == self.ndim
+        assert y.shape[-1] == self.ndim
+        if pairwise:
+            x, y = self._pairwise(x, y)
+        d = self.params[1:] * (x - y) ** 2
+        return self.params[0] * np.exp(-0.5 * np.sum(d, axis=-1))
+
+    def derivatives(self, x, y, pairwise=True):
+        if pairwise:
+            x, y = self._pairwise(x, y)
+        d = self.params[1:] * (x - y) ** 2
+        delta = np.exp(-0.5 * np.sum(d, axis=-1))
+        deltas = -0.5 * (x - y) ** 2 * (delta * self.params[0])[:, :, None]
+        return np.concatenate((np.expand_dims(delta, 0), deltas.T))
+
+    def update_parameters(self, updates):
+        self.params += updates
@@ -0,0 +1,122 @@
+import numpy as np
+
+
+class RelevanceVectorClassifier(object):
+
+    def __init__(self, kernel, alpha=1.):
+        """
+        construct relevance vector classifier
+
+        Parameters
+        ----------
+        kernel : Kernel
+            kernel function to compute components of feature vectors
+        alpha : float
+            initial precision of prior weight distribution
+        """
+        self.kernel = kernel
+        self.alpha = alpha
+
+    def _sigmoid(self, a):
+        return np.tanh(a * 0.5) * 0.5 + 0.5
+
+    def _map_estimate(self, X, t, w, n_iter=10):
+        for _ in range(n_iter):
+            y = self._sigmoid(X @ w)
+            g = X.T @ (y - t) + self.alpha * w
+            H = (X.T * y * (1 - y)) @ X + np.diag(self.alpha)
+            w -= np.linalg.solve(H, g)
+        return w, np.linalg.inv(H)
+
+    def fit(self, X, t, iter_max=100):
+        """
+        maximize evidence with respect ot hyperparameter
+
+        Parameters
+        ----------
+        X : (sample_size, n_features) ndarray
+            input
+        t : (sample_size,) ndarray
+            corresponding target
+        iter_max : int
+            maximum number of iterations
+
+        Attributes
+        ----------
+        X : (N, n_features) ndarray
+            relevance vector
+        t : (N,) ndarray
+            corresponding target
+        alpha : (N,) ndarray
+            hyperparameter for each weight or training sample
+        cov : (N, N) ndarray
+            covariance matrix of weight
+        mean : (N,) ndarray
+            mean of each weight
+        """
+        if X.ndim == 1:
+            X = X[:, None]
+        assert X.ndim == 2
+        assert t.ndim == 1
+        Phi = self.kernel(X, X)
+        N = len(t)
+        self.alpha = np.zeros(N) + self.alpha
+        mean = np.zeros(N)
+        for _ in range(iter_max):
+            param = np.copy(self.alpha)
+            mean, cov = self._map_estimate(Phi, t, mean, 10)
+            gamma = 1 - self.alpha * np.diag(cov)
+            self.alpha = gamma / np.square(mean)
+            np.clip(self.alpha, 0, 1e10, out=self.alpha)
+            if np.allclose(param, self.alpha):
+                break
+        mask = self.alpha < 1e8
+        self.X = X[mask]
+        self.t = t[mask]
+        self.alpha = self.alpha[mask]
+        Phi = self.kernel(self.X, self.X)
+        mean = mean[mask]
+        self.mean, self.covariance = self._map_estimate(Phi, self.t, mean, 100)
+
+    def predict(self, X):
+        """
+        predict class label
+
+        Parameters
+        ----------
+        X : (sample_size, n_features)
+            input
+
+        Returns
+        -------
+        label : (sample_size,) ndarray
+            predicted label
+        """
+        if X.ndim == 1:
+            X = X[:, None]
+        assert X.ndim == 2
+        phi = self.kernel(X, self.X)
+        label = (phi @ self.mean > 0).astype(np.int)
+        return label
+
+    def predict_proba(self, X):
+        """
+        probability of input belonging class one
+
+        Parameters
+        ----------
+        X : (sample_size, n_features) ndarray
+            input
+
+        Returns
+        -------
+        proba : (sample_size,) ndarray
+            probability of predictive distribution p(C1|x)
+        """
+        if X.ndim == 1:
+            X = X[:, None]
+        assert X.ndim == 2
+        phi = self.kernel(X, self.X)
+        mu_a = phi @ self.mean
+        var_a = np.sum(phi @ self.covariance * phi, axis=1)
+        return self._sigmoid(mu_a / np.sqrt(1 + np.pi * var_a / 8))
@@ -0,0 +1,102 @@
+import numpy as np
+
+
+class RelevanceVectorRegressor(object):
+
+    def __init__(self, kernel, alpha=1., beta=1.):
+        """
+        construct relevance vector regressor
+
+        Parameters
+        ----------
+        kernel : Kernel
+            kernel function to compute components of feature vectors
+        alpha : float
+            initial precision of prior weight distribution
+        beta : float
+            precision of observation
+        """
+        self.kernel = kernel
+        self.alpha = alpha
+        self.beta = beta
+
+    def fit(self, X, t, iter_max=1000):
+        """
+        maximize evidence with respect to hyperparameter
+
+        Parameters
+        ----------
+        X : (sample_size, n_features) ndarray
+            input
+        t : (sample_size,) ndarray
+            corresponding target
+        iter_max : int
+            maximum number of iterations
+
+        Attributes
+        -------
+        X : (N, n_features) ndarray
+            relevance vector
+        t : (N,) ndarray
+            corresponding target
+        alpha : (N,) ndarray
+            hyperparameter for each weight or training sample
+        cov : (N, N) ndarray
+            covariance matrix of weight
+        mean : (N,) ndarray
+            mean of each weight
+        """
+        if X.ndim == 1:
+            X = X[:, None]
+        assert X.ndim == 2
+        assert t.ndim == 1
+        N = len(t)
+        Phi = self.kernel(X, X)
+        self.alpha = np.zeros(N) + self.alpha
+        for _ in range(iter_max):
+            params = np.hstack([self.alpha, self.beta])
+            precision = np.diag(self.alpha) + self.beta * Phi.T @ Phi
+            covariance = np.linalg.inv(precision)
+            mean = self.beta * covariance @ Phi.T @ t
+            gamma = 1 - self.alpha * np.diag(covariance)
+            self.alpha = gamma / np.square(mean)
+            np.clip(self.alpha, 0, 1e10, out=self.alpha)
+            self.beta = (N - np.sum(gamma)) / np.sum((t - Phi.dot(mean)) ** 2)
+            if np.allclose(params, np.hstack([self.alpha, self.beta])):
+                break
+        mask = self.alpha < 1e9
+        self.X = X[mask]
+        self.t = t[mask]
+        self.alpha = self.alpha[mask]
+        Phi = self.kernel(self.X, self.X)
+        precision = np.diag(self.alpha) + self.beta * Phi.T @ Phi
+        self.covariance = np.linalg.inv(precision)
+        self.mean = self.beta * self.covariance @ Phi.T @ self.t
+
+    def predict(self, X, with_error=True):
+        """
+        predict output with this model
+
+        Parameters
+        ----------
+        X : (sample_size, n_features)
+            input
+        with_error : bool
+            if True, predict with standard deviation of the outputs
+
+        Returns
+        -------
+        mean : (sample_size,) ndarray
+            mean of predictive distribution
+        std : (sample_size,) ndarray
+            standard deviation of predictive distribution
+        """
+        if X.ndim == 1:
+            X = X[:, None]
+        assert X.ndim == 2
+        phi = self.kernel(X, self.X)
+        mean = phi @ self.mean
+        if with_error:
+            var = 1 / self.beta + np.sum(phi @ self.covariance * phi, axis=1)
+            return mean, np.sqrt(var)
+        return mean
@@ -0,0 +1,107 @@
+import numpy as np
+
+
+class SupportVectorClassifier(object):
+
+    def __init__(self, kernel, C=np.Inf):
+        """
+        construct support vector classifier
+
+        Parameters
+        ----------
+        kernel : Kernel
+            kernel function to compute inner products
+        C : float
+            penalty of misclassification
+        """
+        self.kernel = kernel
+        self.C = C
+
+    def fit(self, X:np.ndarray, t:np.ndarray, tol:float=1e-8):
+        """
+        estimate support vectors and their parameters
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training independent variable
+        t : (N,) np.ndarray
+            training dependent variable
+            binary -1 or 1
+        tol : float, optional
+            numerical tolerance (the default is 1e-8)
+        """
+
+        N = len(t)
+        coef = np.zeros(N)
+        grad = np.ones(N)
+        Gram = self.kernel(X, X)
+
+        while True:
+            tg = t * grad
+            mask_up = (t == 1) & (coef < self.C - tol)
+            mask_up |= (t == -1) & (coef > tol)
+            mask_down = (t == -1) & (coef < self.C - tol)
+            mask_down |= (t == 1) & (coef > tol)
+            i = np.where(mask_up)[0][np.argmax(tg[mask_up])]
+            j = np.where(mask_down)[0][np.argmin(tg[mask_down])]
+            if tg[i] < tg[j] + tol:
+                self.b = 0.5 * (tg[i] + tg[j])
+                break
+            else:
+                A = self.C - coef[i] if t[i] == 1 else coef[i]
+                B = coef[j] if t[j] == 1 else self.C - coef[j]
+                direction = (tg[i] - tg[j]) / (Gram[i, i] - 2 * Gram[i, j] + Gram[j, j])
+                direction = min(A, B, direction)
+                coef[i] += direction * t[i]
+                coef[j] -= direction * t[j]
+                grad -= direction * t * (Gram[i] - Gram[j])
+        support_mask = coef > tol
+        self.a = coef[support_mask]
+        self.X = X[support_mask]
+        self.t = t[support_mask]
+
+    def lagrangian_function(self):
+        return (
+            np.sum(self.a)
+            - self.a
+            @ (self.t * self.t[:, None] * self.kernel(self.X, self.X))
+            @ self.a)
+
+    def predict(self, x):
+        """
+        predict labels of the input
+
+        Parameters
+        ----------
+        x : (sample_size, n_features) ndarray
+            input
+
+        Returns
+        -------
+        label : (sample_size,) ndarray
+            predicted labels
+        """
+        y = self.distance(x)
+        label = np.sign(y)
+        return label
+
+    def distance(self, x):
+        """
+        calculate distance from the decision boundary
+
+        Parameters
+        ----------
+        x : (sample_size, n_features) ndarray
+            input
+
+        Returns
+        -------
+        distance : (sample_size,) ndarray
+            distance from the boundary
+        """
+        distance = np.sum(
+            self.a * self.t
+            * self.kernel(x, self.X),
+            axis=-1) + self.b
+        return distance
@@ -0,0 +1,28 @@
+from prml.linear.bayesian_logistic_regression import BayesianLogisticRegression
+from prml.linear.bayesian_regression import BayesianRegression
+from prml.linear.emprical_bayes_regression import EmpiricalBayesRegression
+from prml.linear.least_squares_classifier import LeastSquaresClassifier
+from prml.linear.linear_regression import LinearRegression
+from prml.linear.fishers_linear_discriminant import FishersLinearDiscriminant
+from prml.linear.logistic_regression import LogisticRegression
+from prml.linear.perceptron import Perceptron
+from prml.linear.ridge_regression import RidgeRegression
+from prml.linear.softmax_regression import SoftmaxRegression
+from prml.linear.variational_linear_regression import VariationalLinearRegression
+from prml.linear.variational_logistic_regression import VariationalLogisticRegression
+
+
+__all__ = [
+    "BayesianLogisticRegression",
+    "BayesianRegression",
+    "EmpiricalBayesRegression",
+    "LeastSquaresClassifier",
+    "LinearRegression",
+    "FishersLinearDiscriminant",
+    "LogisticRegression",
+    "Perceptron",
+    "RidgeRegression",
+    "SoftmaxRegression",
+    "VariationalLinearRegression",
+    "VariationalLogisticRegression"
+]
@@ -0,0 +1,66 @@
+import numpy as np
+from prml.linear.logistic_regression import LogisticRegression
+
+
+class BayesianLogisticRegression(LogisticRegression):
+    """
+    Logistic regression model
+
+    w ~ Gaussian(0, alpha^(-1)I)
+    y = sigmoid(X @ w)
+    t ~ Bernoulli(t|y)
+    """
+
+    def __init__(self, alpha:float=1.):
+        self.alpha = alpha
+
+    def fit(self, X:np.ndarray, t:np.ndarray, max_iter:int=100):
+        """
+        bayesian estimation of logistic regression model
+        using Laplace approximation
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training data independent variable
+        t : (N,) np.ndarray
+            training data dependent variable
+            binary 0 or 1
+        max_iter : int, optional
+            maximum number of paramter update iteration (the default is 100)
+        """
+        w = np.zeros(np.size(X, 1))
+        eye = np.eye(np.size(X, 1))
+        self.w_mean = np.copy(w)
+        self.w_precision = self.alpha * eye
+        for _ in range(max_iter):
+            w_prev = np.copy(w)
+            y = self._sigmoid(X @ w)
+            grad = X.T @ (y - t) + self.w_precision @ (w - self.w_mean)
+            hessian = (X.T * y * (1 - y)) @ X + self.w_precision
+            try:
+                w -= np.linalg.solve(hessian, grad)
+            except np.linalg.LinAlgError:
+                break
+            if np.allclose(w, w_prev):
+                break
+        self.w_mean = w
+        self.w_precision = hessian
+
+    def proba(self, X:np.ndarray):
+        """
+        compute probability of input belonging class 1
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training data independent variable
+
+        Returns
+        -------
+        (N,) np.ndarray
+            probability of positive
+        """
+        mu_a = X @ self.w_mean
+        var_a = np.sum(np.linalg.solve(self.w_precision, X.T).T * X, axis=1)
+        return self._sigmoid(mu_a / np.sqrt(1 + np.pi * var_a / 8))
@@ -0,0 +1,87 @@
+import numpy as np
+from prml.linear.regression import Regression
+
+
+class BayesianRegression(Regression):
+    """
+    Bayesian regression model
+
+    w ~ N(w|0, alpha^(-1)I)
+    y = X @ w
+    t ~ N(t|X @ w, beta^(-1))
+    """
+
+    def __init__(self, alpha:float=1., beta:float=1.):
+        self.alpha = alpha
+        self.beta = beta
+        self.w_mean = None
+        self.w_precision = None
+
+    def _is_prior_defined(self) -> bool:
+        return self.w_mean is not None and self.w_precision is not None
+
+    def _get_prior(self, ndim:int) -> tuple:
+        if self._is_prior_defined():
+            return self.w_mean, self.w_precision
+        else:
+            return np.zeros(ndim), self.alpha * np.eye(ndim)
+
+    def fit(self, X:np.ndarray, t:np.ndarray):
+        """
+        bayesian update of parameters given training dataset
+
+        Parameters
+        ----------
+        X : (N, n_features) np.ndarray
+            training data independent variable
+        t : (N,) np.ndarray
+            training data dependent variable
+        """
+
+        mean_prev, precision_prev = self._get_prior(np.size(X, 1))
+
+        w_precision = precision_prev + self.beta * X.T @ X
+        w_mean = np.linalg.solve(
+            w_precision,
+            precision_prev @ mean_prev + self.beta * X.T @ t
+        )
+        self.w_mean = w_mean
+        self.w_precision = w_precision
+        self.w_cov = np.linalg.inv(self.w_precision)
+
+    def predict(self, X:np.ndarray, return_std:bool=False, sample_size:int=None):
+        """
+        return mean (and standard deviation) of predictive distribution
+
+        Parameters
+        ----------
+        X : (N, n_features) np.ndarray
+            independent variable
+        return_std : bool, optional
+            flag to return standard deviation (the default is False)
+        sample_size : int, optional
+            number of samples to draw from the predictive distribution
+            (the default is None, no sampling from the distribution)
+
+        Returns
+        -------
+        y : (N,) np.ndarray
+            mean of the predictive distribution
+        y_std : (N,) np.ndarray
+            standard deviation of the predictive distribution
+        y_sample : (N, sample_size) np.ndarray
+            samples from the predictive distribution
+        """
+
+        if sample_size is not None:
+            w_sample = np.random.multivariate_normal(
+                self.w_mean, self.w_cov, size=sample_size
+            )
+            y_sample = X @ w_sample.T
+            return y_sample
+        y = X @ self.w_mean
+        if return_std:
+            y_var = 1 / self.beta + np.sum(X @ self.w_cov * X, axis=1)
+            y_std = np.sqrt(y_var)
+            return y, y_std
+        return y
@@ -0,0 +1,5 @@
+class Classifier(object):
+    """
+    Base class for classifiers
+    """
+    pass
@@ -0,0 +1,86 @@
+import numpy as np
+from prml.linear.bayesian_regression import BayesianRegression
+
+
+class EmpiricalBayesRegression(BayesianRegression):
+    """
+    Empirical Bayes Regression model
+    a.k.a.
+    type 2 maximum likelihood,
+    generalized maximum likelihood,
+    evidence approximation
+
+    w ~ N(w|0, alpha^(-1)I)
+    y = X @ w
+    t ~ N(t|X @ w, beta^(-1))
+    evidence function p(t|X,alpha,beta) = S p(t|w;X,beta)p(w|0;alpha) dw
+    """
+
+    def __init__(self, alpha:float=1., beta:float=1.):
+        super().__init__(alpha, beta)
+
+    def fit(self, X:np.ndarray, t:np.ndarray, max_iter:int=100):
+        """
+        maximization of evidence function with respect to
+        the hyperparameters alpha and beta given training dataset
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training independent variable
+        t : (N,) np.ndarray
+            training dependent variable
+        max_iter : int
+            maximum number of iteration
+        """
+        M = X.T @ X
+        eigenvalues = np.linalg.eigvalsh(M)
+        eye = np.eye(np.size(X, 1))
+        N = len(t)
+        for _ in range(max_iter):
+            params = [self.alpha, self.beta]
+
+            w_precision = self.alpha * eye + self.beta * X.T @ X
+            w_mean = self.beta * np.linalg.solve(w_precision, X.T @ t)
+
+            gamma = np.sum(eigenvalues / (self.alpha + eigenvalues))
+            self.alpha = float(gamma / np.sum(w_mean ** 2).clip(min=1e-10))
+            self.beta = float(
+                (N - gamma) / np.sum(np.square(t - X @ w_mean))
+            )
+            if np.allclose(params, [self.alpha, self.beta]):
+                break
+        self.w_mean = w_mean
+        self.w_precision = w_precision
+        self.w_cov = np.linalg.inv(w_precision)
+
+    def _log_prior(self, w):
+        return -0.5 * self.alpha * np.sum(w ** 2)
+
+    def _log_likelihood(self, X, t, w):
+        return -0.5 * self.beta * np.square(t - X @ w).sum()
+
+    def _log_posterior(self, X, t, w):
+        return self._log_likelihood(X, t, w) + self._log_prior(w)
+
+    def log_evidence(self, X:np.ndarray, t:np.ndarray):
+        """
+        logarithm or the evidence function
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            indenpendent variable
+        t : (N,) np.ndarray
+            dependent variable
+        Returns
+        -------
+        float
+            log evidence
+        """
+        N = len(t)
+        D = np.size(X, 1)
+        return 0.5 * (
+            D * np.log(self.alpha) + N * np.log(self.beta)
+            - np.linalg.slogdet(self.w_precision)[1] - D * np.log(2 * np.pi)
+        ) + self._log_posterior(X, t, self.w_mean)
@@ -0,0 +1,80 @@
+import numpy as np
+from prml.linear.classifier import Classifier
+from prml.rv.gaussian import Gaussian
+
+
+class FishersLinearDiscriminant(Classifier):
+    """
+    Fisher's Linear discriminant model
+    """
+
+    def __init__(self, w:np.ndarray=None, threshold:float=None):
+        self.w = w
+        self.threshold = threshold
+
+    def fit(self, X:np.ndarray, t:np.ndarray):
+        """
+        estimate parameter given training dataset
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training dataset independent variable
+        t : (N,) np.ndarray
+            training dataset dependent variable
+            binary 0 or 1
+        """
+        X0 = X[t == 0]
+        X1 = X[t == 1]
+        m0 = np.mean(X0, axis=0)
+        m1 = np.mean(X1, axis=0)
+        cov_inclass = np.cov(X0, rowvar=False) + np.cov(X1, rowvar=False)
+        self.w = np.linalg.solve(cov_inclass, m1 - m0)
+        self.w /= np.linalg.norm(self.w).clip(min=1e-10)
+
+        g0 = Gaussian()
+        g0.fit((X0 @ self.w))
+        g1 = Gaussian()
+        g1.fit((X1 @ self.w))
+        root = np.roots([
+            g1.var - g0.var,
+            2 * (g0.var * g1.mu - g1.var * g0.mu),
+            g1.var * g0.mu ** 2 - g0.var * g1.mu ** 2
+            - g1.var * g0.var * np.log(g1.var / g0.var)
+        ])
+        if g0.mu < root[0] < g1.mu or g1.mu < root[0] < g0.mu:
+            self.threshold = root[0]
+        else:
+            self.threshold = root[1]
+
+    def transform(self, X:np.ndarray):
+        """
+        project data
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            independent variable
+
+        Returns
+        -------
+        y : (N,) np.ndarray
+            projected data
+        """
+        return X @ self.w
+
+    def classify(self, X:np.ndarray):
+        """
+        classify input data
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            independent variable to be classified
+
+        Returns
+        -------
+        (N,) np.ndarray
+            binary class for each input
+        """
+        return (X @ self.w > self.threshold).astype(np.int)
@@ -0,0 +1,48 @@
+import numpy as np
+from prml.linear.classifier import Classifier
+from prml.preprocess.label_transformer import LabelTransformer
+
+
+class LeastSquaresClassifier(Classifier):
+    """
+    Least squares classifier model
+
+    X : (N, D)
+    W : (D, K)
+    y = argmax_k X @ W
+    """
+
+    def __init__(self, W:np.ndarray=None):
+        self.W = W
+
+    def fit(self, X:np.ndarray, t:np.ndarray):
+        """
+        least squares fitting for classification
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training independent variable
+        t : (N,) or (N, K) np.ndarray
+            training dependent variable
+            in class index (N,) or one-of-k coding (N,K)
+        """
+        if t.ndim == 1:
+            t = LabelTransformer().encode(t)
+        self.W = np.linalg.pinv(X) @ t
+
+    def classify(self, X:np.ndarray):
+        """
+        classify input data
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            independent variable to be classified
+
+        Returns
+        -------
+        (N,) np.ndarray
+            class index for each input
+        """
+        return np.argmax(X @ self.W, axis=-1)
@@ -0,0 +1,48 @@
+import numpy as np
+from prml.linear.regression import Regression
+
+
+class LinearRegression(Regression):
+    """
+    Linear regression model
+    y = X @ w
+    t ~ N(t|X @ w, var)
+    """
+
+    def fit(self, X:np.ndarray, t:np.ndarray):
+        """
+        perform least squares fitting
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training independent variable
+        t : (N,) np.ndarray
+            training dependent variable
+        """
+        self.w = np.linalg.pinv(X) @ t
+        self.var = np.mean(np.square(X @ self.w - t))
+
+    def predict(self, X:np.ndarray, return_std:bool=False):
+        """
+        make prediction given input
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            samples to predict their output
+        return_std : bool, optional
+            returns standard deviation of each predition if True
+
+        Returns
+        -------
+        y : (N,) np.ndarray
+            prediction of each sample
+        y_std : (N,) np.ndarray
+            standard deviation of each predition
+        """
+        y = X @ self.w
+        if return_std:
+            y_std = np.sqrt(self.var) + np.zeros_like(y)
+            return y, y_std
+        return y
@@ -0,0 +1,77 @@
+import numpy as np
+from prml.linear.classifier import Classifier
+
+
+class LogisticRegression(Classifier):
+    """
+    Logistic regression model
+
+    y = sigmoid(X @ w)
+    t ~ Bernoulli(t|y)
+    """
+
+    @staticmethod
+    def _sigmoid(a):
+        return np.tanh(a * 0.5) * 0.5 + 0.5
+
+    def fit(self, X:np.ndarray, t:np.ndarray, max_iter:int=100):
+        """
+        maximum likelihood estimation of logistic regression model
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training data independent variable
+        t : (N,) np.ndarray
+            training data dependent variable
+            binary 0 or 1
+        max_iter : int, optional
+            maximum number of paramter update iteration (the default is 100)
+        """
+        w = np.zeros(np.size(X, 1))
+        for _ in range(max_iter):
+            w_prev = np.copy(w)
+            y = self._sigmoid(X @ w)
+            grad = X.T @ (y - t)
+            hessian = (X.T * y * (1 - y)) @ X
+            try:
+                w -= np.linalg.solve(hessian, grad)
+            except np.linalg.LinAlgError:
+                break
+            if np.allclose(w, w_prev):
+                break
+        self.w = w
+
+    def proba(self, X:np.ndarray):
+        """
+        compute probability of input belonging class 1
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training data independent variable
+
+        Returns
+        -------
+        (N,) np.ndarray
+            probability of positive
+        """
+        return self._sigmoid(X @ self.w)
+
+    def classify(self, X:np.ndarray, threshold:float=0.5):
+        """
+        classify input data
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            independent variable to be classified
+        threshold : float, optional
+            threshold of binary classification (default is 0.5)
+
+        Returns
+        -------
+        (N,) np.ndarray
+            binary class for each input
+        """
+        return (self.proba(X) > threshold).astype(np.int)
@@ -0,0 +1,52 @@
+import numpy as np
+from prml.linear.classifier import Classifier
+
+
+class Perceptron(Classifier):
+    """
+    Perceptron model
+    """
+
+    def fit(self, X, t, max_epoch=100):
+        """
+        fit perceptron model on given input pair
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training independent variable
+        t : (N,)
+            training dependent variable
+            binary -1 or 1
+        max_epoch : int, optional
+            maximum number of epoch (the default is 100)
+        """
+        self.w = np.zeros(np.size(X, 1))
+        for _ in range(max_epoch):
+            N = len(t)
+            index = np.random.permutation(N)
+            X = X[index]
+            t = t[index]
+            for x, label in zip(X, t):
+                self.w += x * label
+                if (X @ self.w * t > 0).all():
+                    break
+            else:
+                continue
+            break
+
+    def classify(self, X):
+        """
+        classify input data
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            independent variable to be classified
+
+        Returns
+        -------
+        (N,) np.ndarray
+            binary class (-1 or 1) for each input
+        """
+        return np.sign(X @ self.w).astype(np.int)
@@ -0,0 +1,5 @@
+class Regression(object):
+    """
+    Base class for regressors
+    """
+    pass
@@ -0,0 +1,44 @@
+import numpy as np
+from prml.linear.regression import Regression
+
+
+class RidgeRegression(Regression):
+    """
+    Ridge regression model
+
+    w* = argmin |t - X @ w| + alpha * |w|_2^2
+    """
+
+    def __init__(self, alpha:float=1.):
+        self.alpha = alpha
+
+    def fit(self, X:np.ndarray, t:np.ndarray):
+        """
+        maximum a posteriori estimation of parameter
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training data independent variable
+        t : (N,) np.ndarray
+            training data dependent variable
+        """
+
+        eye = np.eye(np.size(X, 1))
+        self.w = np.linalg.solve(self.alpha * eye + X.T @ X, X.T @ t)
+
+    def predict(self, X:np.ndarray):
+        """
+        make prediction given input
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            samples to predict their output
+
+        Returns
+        -------
+        (N,) np.ndarray
+            prediction of each input
+        """
+        return X @ self.w
@@ -0,0 +1,83 @@
+import numpy as np
+from prml.linear.classifier import Classifier
+from prml.preprocess.label_transformer import LabelTransformer
+
+
+class SoftmaxRegression(Classifier):
+    """
+    Softmax regression model
+    aka
+    multinomial logistic regression,
+    multiclass logistic regression,
+    maximum entropy classifier.
+
+    y = softmax(X @ W)
+    t ~ Categorical(t|y)
+    """
+
+    @staticmethod
+    def _softmax(a):
+        a_max = np.max(a, axis=-1, keepdims=True)
+        exp_a = np.exp(a - a_max)
+        return exp_a / np.sum(exp_a, axis=-1, keepdims=True)
+
+    def fit(self, X:np.ndarray, t:np.ndarray, max_iter:int=100, learning_rate:float=0.1):
+        """
+        maximum likelihood estimation of the parameter
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training independent variable
+        t : (N,) or (N, K) np.ndarray
+            training dependent variable
+            in class index or one-of-k encoding
+        max_iter : int, optional
+            maximum number of iteration (the default is 100)
+        learning_rate : float, optional
+            learning rate of gradient descent (the default is 0.1)
+        """
+        if t.ndim == 1:
+            t = LabelTransformer().encode(t)
+        self.n_classes = np.size(t, 1)
+        W = np.zeros((np.size(X, 1), self.n_classes))
+        for _ in range(max_iter):
+            W_prev = np.copy(W)
+            y = self._softmax(X @ W)
+            grad = X.T @ (y - t)
+            W -= learning_rate * grad
+            if np.allclose(W, W_prev):
+                break
+        self.W = W
+
+    def proba(self, X:np.ndarray):
+        """
+        compute probability of input belonging each class
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            independent variable
+
+        Returns
+        -------
+        (N, K) np.ndarray
+            probability of each class
+        """
+        return self._softmax(X @ self.W)
+
+    def classify(self, X:np.ndarray):
+        """
+        classify input data
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            independent variable to be classified
+
+        Returns
+        -------
+        (N,) np.ndarray
+            class index for each input
+        """
+        return np.argmax(self.proba(X), axis=-1)
@@ -0,0 +1,93 @@
+import numpy as np
+from prml.linear.regression import Regression
+
+
+class VariationalLinearRegression(Regression):
+    """
+    variational bayesian estimation of linear regression model
+    p(w,alpha|X,t)
+    ~ q(w)q(alpha)
+    = N(w|w_mean, w_var)Gamma(alpha|a,b)
+
+    Attributes
+    ----------
+    a : float
+        a parameter of variational posterior gamma distribution
+    b : float
+        another parameter of variational posterior gamma distribution
+    w_mean : (n_features,) ndarray
+        mean of variational posterior gaussian distribution
+    w_var : (n_features, n_feautures) ndarray
+        variance of variational posterior gaussian distribution
+    n_iter : int
+        number of iterations performed
+    """
+
+    def __init__(self, beta:float=1., a0:float=1., b0:float=1.):
+        """
+        construct variational linear regressor
+        Parameters
+        ----------
+        beta : float
+            precision of observation noise
+        a0 : float
+            a parameter of prior gamma distribution
+            Gamma(alpha|a0,b0)
+        b0 : float
+            another parameter of prior gamma distribution
+            Gamma(alpha|a0,b0)
+        """
+        self.beta = beta
+        self.a0 = a0
+        self.b0 = b0
+
+    def fit(self, X:np.ndarray, t:np.ndarray, iter_max:int=100):
+        """
+        variational bayesian estimation of parameter
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training independent variable
+        t : (N,) np.ndarray
+            training dependent variable
+        iter_max : int, optional
+            maximum number of iteration (the default is 100)
+        """
+        D = np.size(X, 1)
+        self.a = self.a0 + 0.5 * D
+        self.b = self.b0
+        I = np.eye(D)
+        for _ in range(iter_max):
+            param = self.b
+            self.w_var = np.linalg.inv(self.a * I / self.b + self.beta * X.T @ X)
+            self.w_mean = self.beta * self.w_var @ X.T @ t
+            self.b = self.b0 + 0.5 * (np.sum(self.w_mean ** 2) + np.trace(self.w_var))
+            if np.allclose(self.b, param):
+                break
+
+    def predict(self, X:np.ndarray, return_std:bool=False):
+        """
+        make prediction of input
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            independent variable
+        return_std : bool, optional
+            return standard deviation of predictive distribution if True
+            (the default is False)
+
+        Returns
+        -------
+        y : (N,) np.ndarray
+            mean of predictive distribution
+        y_std : (N,) np.ndarray
+            standard deviation of predictive distribution
+        """
+        y = X @ self.w_mean
+        if return_std:
+            y_var = 1 / self.beta + np.sum(X @ self.w_var * X, axis=1)
+            y_std = np.sqrt(y_var)
+            return y, y_std
+        return y
@@ -0,0 +1,88 @@
+import numpy as np
+from prml.linear.logistic_regression import LogisticRegression
+
+
+class VariationalLogisticRegression(LogisticRegression):
+
+    def __init__(self, alpha:float=None, a0:float=1., b0:float=1.):
+        """
+        construct variational logistic regressor
+
+        Parameters
+        ----------
+        alpha : float
+            precision parameter of the prior
+            if None, this is also the subject to estimate
+        a0 : float
+            a parameter of hyper prior Gamma dist.
+            Gamma(alpha|a0,b0)
+            if alpha is not None, this argument will be ignored
+        b0 : float
+            another parameter of hyper prior Gamma dist.
+            Gamma(alpha|a0,b0)
+            if alpha is not None, this argument will be ignored
+        """
+        if alpha is not None:
+            self.__alpha = alpha
+        else:
+            self.a0 = a0
+            self.b0 = b0
+
+    def fit(self, X:np.ndarray, t:np.ndarray, iter_max:int=1000):
+        """
+        variational bayesian estimation of the parameter
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training independent variable
+        t : (N,) np.ndarray
+            training dependent variable
+        iter_max : int, optional
+            maximum number of iteration (the default is 1000)
+        """
+        N, D = X.shape
+        if hasattr(self, "a0"):
+            self.a = self.a0 + 0.5 * D
+        xi = np.random.uniform(-1, 1, size=N)
+        I = np.eye(D)
+        param = np.copy(xi)
+        for _ in range(iter_max):
+            lambda_ = np.tanh(xi) * 0.25 / xi
+            self.w_var = np.linalg.inv(I / self.alpha + 2 * (lambda_ * X.T) @ X)
+            self.w_mean = self.w_var @ np.sum(X.T * (t - 0.5), axis=1)
+            xi = np.sqrt(np.sum(X @ (self.w_var + self.w_mean * self.w_mean[:, None]) * X, axis=-1))
+            if np.allclose(xi, param):
+                break
+            else:
+                param = np.copy(xi)
+
+    @property
+    def alpha(self):
+        if hasattr(self, "__alpha"):
+            return self.__alpha
+        else:
+            try:
+                self.b = self.b0 + 0.5 * (np.sum(self.w_mean ** 2) + np.trace(self.w_var))
+            except AttributeError:
+                self.b = self.b0
+            return self.a / self.b
+
+    def proba(self, X:np.ndarray):
+        """
+        compute probability of input belonging class 1
+
+        Parameters
+        ----------
+        X : (N, D) np.ndarray
+            training data independent variable
+
+        Returns
+        -------
+        (N,) np.ndarray
+            probability of positive
+        """
+        mu_a = X @ self.w_mean
+        var_a = np.sum(X @ self.w_var * X, axis=1)
+        y = self._sigmoid(mu_a / np.sqrt(1 + np.pi * var_a / 8))
+        return y
@@ -0,0 +1,14 @@
+from .categorical_hmm import CategoricalHMM
+from .gaussian_hmm import GaussianHMM
+from prml.markov.kalman import Kalman, kalman_filter, kalman_smoother
+from .particle import Particle
+
+
+__all__ = [
+    "GaussianHMM",
+    "CategoricalHMM",
+    "Kalman",
+    "kalman_filter",
+    "kalman_smoother",
+    "Particle"
+]
@@ -0,0 +1,65 @@
+import numpy as np
+from .hmm import HiddenMarkovModel
+
+
+class CategoricalHMM(HiddenMarkovModel):
+    """
+    Hidden Markov Model with categorical emission model
+    """
+
+    def __init__(self, initial_proba, transition_proba, means):
+        """
+        construct hidden markov model with categorical emission model
+
+        Parameters
+        ----------
+        initial_proba : (n_hidden,) np.ndarray
+            probability of initial latent state
+        transition_proba : (n_hidden, n_hidden) np.ndarray
+            transition probability matrix
+            (i, j) component denotes the transition probability from i-th to j-th hidden state
+        means : (n_hidden, ndim) np.ndarray
+            mean parameters of categorical distribution
+
+        Returns
+        -------
+        ndim : int
+            number of observation categories
+        n_hidden : int
+            number of hidden states
+        """
+        assert initial_proba.size == transition_proba.shape[0] == transition_proba.shape[1] == means.shape[0]
+        assert np.allclose(means.sum(axis=1), 1)
+        super().__init__(initial_proba, transition_proba)
+        self.ndim = means.shape[1]
+        self.means = means
+
+    def draw(self, n=100):
+        """
+        draw random sequence from this model
+
+        Parameters
+        ----------
+        n : int
+            length of the random sequence
+
+        Returns
+        -------
+        seq : (n,) np.ndarray
+            generated random sequence
+        """
+        hidden_state = np.random.choice(self.n_hidden, p=self.initial_proba)
+        seq = []
+        while len(seq) < n:
+            seq.append(np.random.choice(self.ndim, p=self.means[hidden_state]))
+            hidden_state = np.random.choice(self.n_hidden, p=self.transition_proba[hidden_state])
+        return np.asarray(seq)
+
+    def likelihood(self, X):
+        return self.means[X]
+
+    def maximize(self, seq, p_hidden, p_transition):
+        self.initial_proba = p_hidden[0] / np.sum(p_hidden[0])
+        self.transition_proba = np.sum(p_transition, axis=0) / np.sum(p_transition, axis=(0, 2))
+        x = p_hidden[:, None, :] * (np.eye(self.ndim)[seq])[:, :, None]
+        self.means = np.sum(x, axis=0) / np.sum(p_hidden, axis=0)
@@ -0,0 +1,76 @@
+import numpy as np
+from prml.rv import MultivariateGaussian
+from .hmm import HiddenMarkovModel
+
+
+class GaussianHMM(HiddenMarkovModel):
+    """
+    Hidden Markov Model with Gaussian emission model
+    """
+
+    def __init__(self, initial_proba, transition_proba, means, covs):
+        """
+        construct hidden markov model with Gaussian emission model
+
+        Parameters
+        ----------
+        initial_proba : (n_hidden,) np.ndarray or None
+            probability of initial states
+        transition_proba : (n_hidden, n_hidden) np.ndarray or None
+            transition probability matrix
+            (i, j) component denotes the transition probability from i-th to j-th hidden state
+        means : (n_hidden, ndim) np.ndarray
+            mean of each gaussian component
+        covs : (n_hidden, ndim, ndim) np.ndarray
+            covariance matrix of each gaussian component
+
+        Attributes
+        ----------
+        ndim : int
+            dimensionality of observation space
+        n_hidden : int
+            number of hidden states
+        """
+        assert initial_proba.size == transition_proba.shape[0] == transition_proba.shape[1] == means.shape[0] == covs.shape[0]
+        assert means.shape[1] == covs.shape[1] == covs.shape[2]
+        super().__init__(initial_proba, transition_proba)
+        self.ndim = means.shape[1]
+        self.means = means
+        self.covs = covs
+        self.precisions = np.linalg.inv(self.covs)
+        self.gaussians = [MultivariateGaussian(m, cov) for m, cov in zip(means, covs)]
+
+    def draw(self, n=100):
+        """
+        draw random sequence from this model
+
+        Parameters
+        ----------
+        n : int
+            length of the random sequence
+
+        Returns
+        -------
+        seq : (n, ndim) np.ndarray
+            generated random sequence
+        """
+        hidden_state = np.random.choice(self.n_hidden, p=self.initial_proba)
+        seq = []
+        while len(seq) < n:
+            seq.extend(self.gaussians[hidden_state].draw())
+            hidden_state = np.random.choice(self.n_hidden, p=self.transition_proba[hidden_state])
+        return np.asarray(seq)
+
+    def likelihood(self, X):
+        diff = X[:, None, :] - self.means
+        exponents = np.sum(
+            np.einsum('nki,kij->nkj', diff, self.precisions) * diff, axis=-1)
+        return np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(self.covs) * (2 * np.pi) ** self.ndim)
+
+    def maximize(self, seq, p_hidden, p_transition):
+        self.initial_proba = p_hidden[0] / np.sum(p_hidden[0])
+        self.transition_proba = np.sum(p_transition, axis=0) / np.sum(p_transition, axis=(0, 2))
+        Nk = np.sum(p_hidden, axis=0)
+        self.means = (seq.T @ p_hidden / Nk).T
+        diffs = seq[:, None, :] - self.means
+        self.covs = np.einsum('nki,nkj->kij', diffs, diffs * p_hidden[:, :, None]) / Nk[:, None, None]
@@ -0,0 +1,178 @@
+import numpy as np
+
+
+class HiddenMarkovModel(object):
+    """
+    Base class of Hidden Markov models
+    """
+
+    def __init__(self, initial_proba, transition_proba):
+        """
+        construct hidden markov model
+
+        Parameters
+        ----------
+        initial_proba : (n_hidden,) np.ndarray
+            initial probability of each hidden state
+        transition_proba : (n_hidden, n_hidden) np.ndarray
+            transition probability matrix
+            (i, j) component denotes the transition probability from i-th to j-th hidden state
+
+        Attribute
+        ---------
+        n_hidden : int
+            number of hidden state
+        """
+        self.n_hidden = initial_proba.size
+        self.initial_proba = initial_proba
+        self.transition_proba = transition_proba
+
+    def fit(self, seq, iter_max=100):
+        """
+        perform EM algorithm to estimate parameter of emission model and hidden variables
+
+        Parameters
+        ----------
+        seq : (N, ndim) np.ndarray
+            observed sequence
+        iter_max : int
+            maximum number of EM steps
+
+        Returns
+        -------
+        posterior : (N, n_hidden) np.ndarray
+            posterior distribution of each latent variable
+        """
+        params = np.hstack(
+            (self.initial_proba.ravel(), self.transition_proba.ravel()))
+        for i in range(iter_max):
+            p_hidden, p_transition = self.expect(seq)
+            self.maximize(seq, p_hidden, p_transition)
+            params_new = np.hstack(
+                (self.initial_proba.ravel(), self.transition_proba.ravel()))
+            if np.allclose(params, params_new):
+                break
+            else:
+                params = params_new
+        return self.forward_backward(seq)
+
+    def expect(self, seq):
+        """
+        estimate posterior distributions of hidden states and
+        transition probability between adjacent latent variables
+
+        Parameters
+        ----------
+        seq : (N, ndim) np.ndarray
+            observed sequence
+
+        Returns
+        -------
+        p_hidden : (N, n_hidden) np.ndarray
+            posterior distribution of each hidden variable
+        p_transition : (N - 1, n_hidden, n_hidden) np.ndarray
+            posterior transition probability between adjacent latent variables
+        """
+        likelihood = self.likelihood(seq)
+
+        f = self.initial_proba * likelihood[0]
+        constant = [f.sum()]
+        forward = [f / f.sum()]
+        for like in likelihood[1:]:
+            f = forward[-1] @ self.transition_proba * like
+            constant.append(f.sum())
+            forward.append(f / f.sum())
+        forward = np.asarray(forward)
+        constant = np.asarray(constant)
+
+        backward = [np.ones(self.n_hidden)]
+        for like, c in zip(likelihood[-1:0:-1], constant[-1:0:-1]):
+            backward.insert(0, self.transition_proba @ (like * backward[0]) / c)
+        backward = np.asarray(backward)
+
+        p_hidden = forward * backward
+        p_transition = self.transition_proba * likelihood[1:, None, :] * backward[1:, None, :] * forward[:-1, :, None]
+        return p_hidden, p_transition
+
+    def forward_backward(self, seq):
+        """
+        estimate posterior distributions of hidden states
+
+        Parameters
+        ----------
+        seq : (N, ndim) np.ndarray
+            observed sequence
+
+        Returns
+        -------
+        posterior : (N, n_hidden) np.ndarray
+            posterior distribution of hidden states
+        """
+        likelihood = self.likelihood(seq)
+
+        f = self.initial_proba * likelihood[0]
+        constant = [f.sum()]
+        forward = [f / f.sum()]
+        for like in likelihood[1:]:
+            f = forward[-1] @ self.transition_proba * like
+            constant.append(f.sum())
+            forward.append(f / f.sum())
+
+        backward = [np.ones(self.n_hidden)]
+        for like, c in zip(likelihood[-1:0:-1], constant[-1:0:-1]):
+            backward.insert(0, self.transition_proba @ (like * backward[0]) / c)
+
+        forward = np.asarray(forward)
+        backward = np.asarray(backward)
+        posterior = forward * backward
+        return posterior
+
+    def filtering(self, seq):
+        """
+        bayesian filtering
+
+        Parameters
+        ----------
+        seq : (N, ndim) np.ndarray
+            observed sequence
+
+        Returns
+        -------
+        posterior : (N, n_hidden) np.ndarray
+            posterior distributions of each latent variables
+        """
+        likelihood = self.likelihood(seq)
+        p = self.initial_proba * likelihood[0]
+        posterior = [p / np.sum(p)]
+        for like in likelihood[1:]:
+            p = posterior[-1] @ self.transition_proba * like
+            posterior.append(p / np.sum(p))
+        posterior = np.asarray(posterior)
+        return posterior
+
+    def viterbi(self, seq):
+        """
+        viterbi algorithm (a.k.a. max-sum algorithm)
+
+        Parameters
+        ----------
+        seq : (N, ndim) np.ndarray
+            observed sequence
+
+        Returns
+        -------
+        seq_hid : (N,) np.ndarray
+            the most probable sequence of hidden variables
+        """
+        nll = -np.log(self.likelihood(seq))
+        cost_total = nll[0]
+        from_list = []
+        for i in range(1, len(seq)):
+            cost_temp = cost_total[:, None] - np.log(self.transition_proba) + nll[i]
+            cost_total = np.min(cost_temp, axis=0)
+            index = np.argmin(cost_temp, axis=0)
+            from_list.append(index)
+        seq_hid = [np.argmin(cost_total)]
+        for source in from_list[::-1]:
+            seq_hid.insert(0, source[seq_hid[0]])
+        return seq_hid
@@ -0,0 +1,268 @@
+import numpy as np
+from prml.rv.multivariate_gaussian import MultivariateGaussian as Gaussian
+from prml.markov.state_space_model import StateSpaceModel
+
+
+class Kalman(StateSpaceModel):
+    """
+    A class to perform kalman filtering or smoothing
+    z : internal state
+    x : observation
+
+    z_1 ~ N(z_1|mu_0, P_0)\n
+    z_n ~ N(z_n|A z_n-1, P)\n
+    x_n ~ N(x_n|C z_n, S)
+
+    Parameters
+    ----------
+    system : (Dz, Dz) np.ndarray
+        system matrix aka transition matrix (A)
+    cov_system : (Dz, Dz) np.ndarray
+        covariance matrix of process noise
+    measure : (Dx, Dz) np.ndarray
+        measurement matrix aka observation matrix (C)
+    cov_measure : (Dx, Dx) np.ndarray
+        covariance matrix of measurement noise
+    mu0 : (Dz,) np.ndarray
+        mean parameter of initial hidden variable
+    P0 : (Dz, Dz) np.ndarray
+        covariance parameter of initial hidden variable
+
+    Attributes
+    ----------
+    Dz : int
+        dimensionality of hidden variable
+    Dx : int
+        dimensionality of observed variable
+    """
+
+
+    def __init__(self, system, cov_system, measure, cov_measure, mu0, P0):
+        """
+        construct Kalman model
+
+        z_1 ~ N(z_1|mu_0, P_0)\n
+        z_n ~ N(z_n|A z_n-1, P)\n
+        x_n ~ N(x_n|C z_n, S)
+
+        Parameters
+        ----------
+        system : (Dz, Dz) np.ndarray
+            system matrix aka transition matrix (A)
+        cov_system : (Dz, Dz) np.ndarray
+            covariance matrix of process noise
+        measure : (Dx, Dz) np.ndarray
+            measurement matrix aka observation matrix (C)
+        cov_measure : (Dx, Dx) np.ndarray
+            covariance matrix of measurement noise
+        mu0 : (Dz,) np.ndarray
+            mean parameter of initial hidden variable
+        P0 : (Dz, Dz) np.ndarray
+            covariance parameter of initial hidden variable
+
+        Attributes
+        ----------
+        hidden_mean : list of (Dz,) np.ndarray
+            list of mean of hidden state starting from the given hidden state
+        hidden_cov : list of (Dz, Dz) np.ndarray
+            list of covariance of hidden state starting from the given hidden state
+        """
+        self.system = system
+        self.cov_system = cov_system
+        self.measure = measure
+        self.cov_measure = cov_measure
+
+        self.hidden_mean = [mu0]
+        self.hidden_cov = [P0]
+        self.hidden_cov_predicted = [None]
+
+        self.smoothed_until = -1
+        self.smoothing_gain = [None]
+
+    def predict(self):
+        """
+        predict hidden state at current step given estimate at previous step
+
+        Returns
+        -------
+        tuple ((Dz,) np.ndarray, (Dz, Dz) np.ndarray)
+            tuple of mean and covariance of the estimate at current step
+        """
+        mu_prev, cov_prev = self.hidden_mean[-1], self.hidden_cov[-1]
+        mu = self.system @ mu_prev
+        cov = self.system @ cov_prev @ self.system.T + self.cov_system
+        self.hidden_mean.append(mu)
+        self.hidden_cov.append(cov)
+        self.hidden_cov_predicted.append(np.copy(cov))
+        return mu, cov
+
+    def filter(self, observed):
+        """
+        bayesian update of current estimate given current observation
+
+        Parameters
+        ----------
+        observed : (Dx,) np.ndarray
+            current observation
+
+        Returns
+        -------
+        tuple ((Dz,) np.ndarray, (Dz, Dz) np.ndarray)
+            tuple of mean and covariance of the updated estimate
+        """
+        mu, cov = self.hidden_mean[-1], self.hidden_cov[-1]
+        innovation = observed - self.measure @ mu
+        cov_innovation = self.cov_measure + self.measure @ cov @ self.measure.T
+        kalman_gain = np.linalg.solve(cov_innovation, self.measure @ cov).T
+        mu += kalman_gain @ innovation
+        cov -= kalman_gain @ self.measure @ cov
+        return mu, cov
+
+    def filtering(self, observed_sequence):
+        """
+        perform kalman filtering given observed sequence
+
+        Parameters
+        ----------
+        observed_sequence : (T, Dx) np.ndarray
+            sequence of observations
+
+        Returns
+        -------
+        tuple ((T, Dz) np.ndarray, (T, Dz, Dz) np.ndarray)
+            seuquence of mean and covariance of hidden variable at each time step
+        """
+        for obs in observed_sequence:
+            self.predict()
+            self.filter(obs)
+        mean_sequence = np.asarray(self.hidden_mean[1:])
+        cov_sequence = np.asarray(self.hidden_cov[1:])
+        return mean_sequence, cov_sequence
+
+    def smooth(self):
+        """
+        bayesian update of current estimate with future observations
+        """
+        mean_smoothed_next = self.hidden_mean[self.smoothed_until]
+        cov_smoothed_next = self.hidden_cov[self.smoothed_until]
+        cov_pred_next = self.hidden_cov_predicted[self.smoothed_until]
+
+        self.smoothed_until -= 1
+        mean = self.hidden_mean[self.smoothed_until]
+        cov = self.hidden_cov[self.smoothed_until]
+        gain = np.linalg.solve(cov_pred_next, self.system @ cov).T
+        mean += gain @ (mean_smoothed_next - self.system @ mean)
+        cov += gain @ (cov_smoothed_next - cov_pred_next) @ gain.T
+        self.smoothing_gain.insert(0, gain)
+
+    def smoothing(self, observed_sequence:np.ndarray=None):
+        """
+        perform Kalman smoothing (given observed sequence)
+
+        Parameters
+        ----------
+        observed_sequence : (T, Dx) np.ndarray, optional
+            sequence of observation
+            run Kalman filter if given (the default is None)
+
+        Returns
+        -------
+        tuple ((T, Dz) np.ndarray, (T, Dz, Dz) np.ndarray)
+            sequence of mean and covariance of hidden variable at each time step
+        """
+        if observed_sequence is not None:
+            self.filtering(observed_sequence)
+        while self.smoothed_until != -len(self.hidden_mean):
+            self.smooth()
+        mean_sequence = np.asarray(self.hidden_mean[1:])
+        cov_sequence = np.asarray(self.hidden_cov[1:])
+        return mean_sequence, cov_sequence
+
+    def update_parameter(self, observation_sequence):
+        """
+        maximization step of EM algorithm
+        """
+        mu0 = self.hidden_mean[1]
+        P0 = self.hidden_cov[1]
+
+        Ezn = np.asarray(self.hidden_mean)
+        Eznzn = np.asarray(self.hidden_cov) + Ezn[..., None] * Ezn[:, None, :]
+        Eznzn_1 = np.einsum("nij,nkj->nik", self.hidden_cov[2:], self.smoothing_gain[1:-1]) + Ezn[2:, :, None] * Ezn[1:-1, None, :]
+        self.system = np.linalg.solve(np.sum(Eznzn[2:], axis=0), np.sum(Eznzn_1, axis=0).T).T
+        self.cov_system = np.mean(
+            Eznzn[2:]
+            - np.einsum("ij,nkj->nik", self.system, Eznzn_1)
+            - np.einsum("nij,kj->nik", Eznzn_1, self.system)
+            + np.einsum("ij,njk,lk->nil", self.system, Eznzn[1:-1], self.system),
+            axis=0
+        )
+        self.measure = np.linalg.solve(
+            np.sum(Eznzn[1:], axis=0),
+            np.sum(np.einsum("ni,nj->nij", Ezn[1:], observation_sequence), axis=0)
+        ).T
+        self.cov_measure = np.mean(
+            np.einsum("ni,nj->nij", observation_sequence, observation_sequence)
+            - np.einsum("ij,nj,nk->nik", self.measure, Ezn[1:], observation_sequence)
+            - np.einsum("ni,nj,kj->nik", observation_sequence, Ezn[1:], self.measure)
+            + np.einsum("ij,njk,lk->nil", self.measure, Eznzn[1:], self.measure),
+            axis=0
+        )
+        return self.system, self.cov_system, self.measure, self.cov_measure, mu0, P0
+
+    def fit(self, sequence, max_iter=10):
+        for _ in range(max_iter):
+            kalman_smoother(self, sequence)
+            param = self.update_parameter(sequence)
+            self.__init__(*param)
+        return kalman_smoother(self, sequence)
+
+
+def kalman_filter(kalman:Kalman, observed_sequence:np.ndarray)->tuple:
+    """
+    perform kalman filtering given Kalman model and observed sequence
+
+    Parameters
+    ----------
+    kalman : Kalman
+        Kalman model
+    observed_sequence : (T, Dx) np.ndarray
+        sequence of observations
+
+    Returns
+    -------
+    tuple ((T, Dz) np.ndarray, (T, Dz, Dz) np.ndarray)
+        seuquence of mean and covariance of hidden variable at each time step
+    """
+    for obs in observed_sequence:
+        kalman.predict()
+        kalman.filter(obs)
+    mean_sequence = np.asarray(kalman.hidden_mean[1:])
+    cov_sequence = np.asarray(kalman.hidden_cov[1:])
+    return mean_sequence, cov_sequence
+
+
+def kalman_smoother(kalman:Kalman, observed_sequence:np.ndarray=None):
+    """
+    perform Kalman smoothing given Kalman model (and observed sequence)
+
+    Parameters
+    ----------
+    kalman : Kalman
+        Kalman model
+    observed_sequence : (T, Dx) np.ndarray, optional
+        sequence of observation
+        run Kalman filter if given (the default is None)
+
+    Returns
+    -------
+    tuple ((T, Dz) np.ndarray, (T, Dz, Dz) np.ndarray)
+        seuqnce of mean and covariance of hidden variable at each time step
+    """
+
+    if observed_sequence is not None:
+        kalman_filter(kalman, observed_sequence)
+    while kalman.smoothed_until != -len(kalman.hidden_mean):
+        kalman.smooth()
+    mean_sequence = np.asarray(kalman.hidden_mean[1:])
+    cov_sequence = np.asarray(kalman.hidden_cov[1:])
+    return mean_sequence, cov_sequence
@@ -0,0 +1,124 @@
+import numpy as np
+from scipy.misc import logsumexp
+from scipy.spatial.distance import cdist
+from .state_space_model import StateSpaceModel
+
+
+class Particle(StateSpaceModel):
+    """
+    A class to perform particle filtering, smoothing
+
+    z_1 ~ p(z_1)\n
+    z_n ~ p(z_n|z_n-1)\n
+    x_n ~ p(x_n|z_n)
+
+    Parameters
+    ----------
+    init_particle : (n_particle, ndim_hidden)
+        initial hidden state
+    sampler : callable (particles)
+        function to sample particles at current step given previous state
+    nll : callable (observation, particles)
+        function to compute negative log likelihood for each particle
+
+    Attribute
+    ---------
+    hidden_state : list of (n_paticle, ndim_hidden) np.ndarray
+        list of particles
+    """
+
+    def __init__(self, init_particle, system, cov_system, nll, pdf=None):
+        """
+        construct state space model to perform particle filtering or smoothing
+
+        Parameters
+        ----------
+        init_particle : (n_particle, ndim_hidden) np.ndarray
+            initial hidden state
+        system : (ndim_hidden, ndim_hidden) np.ndarray
+            system matrix aka transition matrix
+        cov_system : (ndim_hidden, ndim_hidden) np.ndarray
+            covariance matrix of process noise
+        nll : callable (observation, particles)
+            function to compute negative log likelihood for each particle
+
+        Attribute
+        ---------
+        particle : list of (n_paticle, ndim_hidden) np.ndarray
+            list of particles at each step
+        weight : list of (n_particle,) np.ndarray
+            list of importance of each particle at each step
+        n_particle : int
+            number of particles at each step
+        """
+        self.particle = [init_particle]
+        self.n_particle, self.ndim_hidden = init_particle.shape
+        self.weight = [np.ones(self.n_particle) / self.n_particle]
+        self.system = system
+        self.cov_system = cov_system
+        self.nll = nll
+        self.smoothed_until = -1
+
+    def resample(self):
+        index = np.random.choice(self.n_particle, self.n_particle, p=self.weight[-1])
+        return self.particle[-1][index]
+
+    def predict(self):
+        predicted = self.resample() @ self.system.T
+        predicted += np.random.multivariate_normal(np.zeros(self.ndim_hidden), self.cov_system, self.n_particle)
+        self.particle.append(predicted)
+        self.weight.append(np.ones(self.n_particle) / self.n_particle)
+        return predicted, self.weight[-1]
+
+    def weigh(self, observed):
+        logit = -self.nll(observed, self.particle[-1])
+        logit -= logsumexp(logit)
+        self.weight[-1] = np.exp(logit)
+
+    def filter(self, observed):
+        self.weigh(observed)
+        return self.particle[-1], self.weight[-1]
+
+    def filtering(self, observed_sequence):
+        mean = []
+        cov = []
+        for obs in observed_sequence:
+            self.predict()
+            p, w = self.filter(obs)
+            mean.append(np.average(p, axis=0, weights=w))
+            cov.append(np.cov(p, rowvar=False, aweights=w))
+        return np.asarray(mean), np.asarray(cov)
+
+    def transition_probability(self, particle, particle_prev):
+        dist = cdist(
+            particle,
+            particle_prev @ self.system.T,
+            "mahalanobis",
+            VI=np.linalg.inv(self.cov_system))
+        matrix = np.exp(-0.5 * np.square(dist))
+        matrix /= np.sum(matrix, axis=1, keepdims=True)
+        matrix[np.isnan(matrix)] = 1 / self.n_particle
+        return matrix
+
+    def smooth(self):
+        particle_next = self.particle[self.smoothed_until]
+        weight_next = self.weight[self.smoothed_until]
+
+        self.smoothed_until -= 1
+        particle = self.particle[self.smoothed_until]
+        weight = self.weight[self.smoothed_until]
+        matrix = self.transition_probability(particle_next, particle).T
+        weight *= matrix @ weight_next / (weight @ matrix)
+        weight /= np.sum(weight, keepdims=True)
+
+    def smoothing(self, observed_sequence:np.ndarray=None):
+        if observed_sequence is not None:
+            self.filtering(observed_sequence)
+        while self.smoothed_until != -len(self.particle):
+            self.smooth()
+        mean = []
+        cov = []
+        for p, w in zip(self.particle, self.weight):
+            mean.append(np.average(p, axis=0, weights=w))
+            cov.append(np.cov(p, rowvar=False, aweights=w))
+        return np.asarray(mean), np.asarray(cov)
@@ -0,0 +1,5 @@
+class StateSpaceModel(object):
+    """
+    Base class for state-space models
+    """
+    pass
@@ -0,0 +1,34 @@
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.parameter import Parameter
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.array.flatten import flatten
+from prml.nn.array.reshape import reshape
+from prml.nn.array.split import split
+from prml.nn.array.transpose import transpose
+from prml.nn import linalg
+from prml.nn.image.convolve2d import convolve2d
+from prml.nn.image.max_pooling2d import max_pooling2d
+from prml.nn.math.abs import abs
+from prml.nn.math.exp import exp
+from prml.nn.math.gamma import gamma
+from prml.nn.math.log import log
+from prml.nn.math.mean import mean
+from prml.nn.math.power import power
+from prml.nn.math.product import prod
+from prml.nn.math.sqrt import sqrt
+from prml.nn.math.square import square
+from prml.nn.math.sum import sum
+from prml.nn.nonlinear.relu import relu
+from prml.nn.nonlinear.sigmoid import sigmoid
+from prml.nn.nonlinear.softmax import softmax
+from prml.nn.nonlinear.softplus import softplus
+from prml.nn.nonlinear.tanh import tanh
+from prml.nn import optimizer
+from prml.nn import random
+from prml.nn.network import Network
+
+
+__all__ = [
+    "optimizer",
+    "Network"
+]
@@ -0,0 +1,19 @@
+from prml.nn.array.broadcast import broadcast_to
+from prml.nn.array.flatten import flatten
+from prml.nn.array.reshape import reshape, reshape_method
+from prml.nn.array.split import split
+from prml.nn.array.transpose import transpose, transpose_method
+from prml.nn.tensor.tensor import Tensor
+
+
+Tensor.flatten = flatten
+Tensor.reshape = reshape_method
+Tensor.transpose = transpose_method
+
+__all__ = [
+    "broadcast_to",
+    "flatten",
+    "reshape",
+    "split",
+    "transpose"
+]
@@ -0,0 +1,36 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class BroadcastTo(Function):
+    """
+    Broadcast a tensor to an new shape
+    """
+
+    def forward(self, x, shape):
+        x = self._convert2tensor(x)
+        self.x = x
+        output = np.broadcast_to(x.value, shape)
+        if isinstance(self.x, Constant):
+            return Constant(output)
+        return Tensor(output, function=self)
+
+    def backward(self, delta):
+        dx = delta
+        if delta.ndim != self.x.ndim:
+            dx = dx.sum(axis=tuple(range(dx.ndim - self.x.ndim)))
+            if isinstance(dx, np.number):
+                dx = np.array(dx)
+        axis = tuple(i for i, len_ in enumerate(self.x.shape) if len_ == 1)
+        if axis:
+            dx = dx.sum(axis=axis, keepdims=True)
+        self.x.backward(dx)
+
+
+def broadcast_to(x, shape):
+    """
+    Broadcast a tensor to an new shape
+    """
+    return BroadcastTo().forward(x, shape)
@@ -0,0 +1,28 @@
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Flatten(Function):
+    """
+    flatten array
+    """
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self._atleast_ndim(x, 2)
+        self.x = x
+        if isinstance(self.x, Constant):
+            return Constant(x.value.flatten())
+        return Tensor(x.value.flatten(), function=self)
+
+    def backward(self, delta):
+        dx = delta.reshape(*self.x.shape)
+        self.x.backward(dx)
+
+
+def flatten(x):
+    """
+    flatten N-dimensional array (N >= 2)
+    """
+    return Flatten().forward(x)
@@ -0,0 +1,32 @@
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Reshape(Function):
+    """
+    reshape array
+    """
+
+    def forward(self, x, shape):
+        x = self._convert2tensor(x)
+        self._atleast_ndim(x, 1)
+        self.x = x
+        if isinstance(self.x, Constant):
+            return Constant(x.value.reshape(*shape))
+        return Tensor(x.value.reshape(*shape), function=self)
+
+    def backward(self, delta):
+        dx = delta.reshape(*self.x.shape)
+        self.x.backward(dx)
+
+
+def reshape(x, shape):
+    """
+    reshape N-dimensional array (N >= 1)
+    """
+    return Reshape().forward(x, shape)
+
+
+def reshape_method(x, *args):
+    return Reshape().forward(x, args)
@@ -0,0 +1,48 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Nth(Function):
+
+    def __init__(self, n):
+        self.n = n
+
+    def forward(self, x):
+        self.x = x
+        if isinstance(self.x, Constant):
+            return Constant(x.value)
+        return Tensor(x.value, function=self)
+
+    def backward(self, delta):
+        self.x.backward(delta, n=self.n)
+
+
+class Split(Function):
+
+    def __init__(self, indices_or_sections, axis=-1):
+        self.indices_or_sections = indices_or_sections
+        self.axis = axis
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self._atleast_ndim(x, 1)
+        self.x = x
+        output = np.split(x.value, self.indices_or_sections, self.axis)
+        if isinstance(self.x, Constant):
+            return tuple([Constant(out) for out in output])
+        self.n_output = len(output)
+        self.delta = [None for _ in output]
+        return tuple([Tensor(out, function=self) for out in output])
+
+    def backward(self, delta, n):
+        self.delta[n] = delta
+        if all([d is not None for d in self.delta]):
+            dx = np.concatenate(self.delta, axis=self.axis)
+            self.x.backward(dx)
+
+
+def split(x, indices_or_sections, axis=-1):
+    output = Split(indices_or_sections, axis).forward(x)
+    return tuple([Nth(i).forward(out) for i, out in enumerate(output)])
@@ -0,0 +1,36 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Transpose(Function):
+
+    def __init__(self, axes=None):
+        self.axes = axes
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        if self.axes is not None:
+            self._equal_ndim(x, len(self.axes))
+        self.x = x
+        if isinstance(self.x, Constant):
+            return Constant(np.transpose(x.value, self.axes))
+        return Tensor(np.transpose(x.value, self.axes), function=self)
+
+    def backward(self, delta):
+        if self.axes is None:
+            dx = np.transpose(delta)
+        else:
+            dx = np.transpose(delta, np.argsort(self.axes))
+        self.x.backward(dx)
+
+
+def transpose(x, axes=None):
+    return Transpose(axes).forward(x)
+
+
+def transpose_method(x, *args):
+    if args == ():
+        args = None
+    return Transpose(args).forward(x)
@@ -0,0 +1,33 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+
+
+class Function(object):
+    """
+    Base class for differentiable functions
+    """
+
+    def _convert2tensor(self, x):
+        if isinstance(x, (int, float, np.number, np.ndarray)):
+            x = Constant(x)
+        elif not isinstance(x, Tensor):
+            raise TypeError(
+                "Unsupported class for input: {}".format(type(x))
+            )
+        return x
+
+    def _equal_ndim(self, x, ndim):
+        if x.ndim != ndim:
+            raise ValueError(
+                "dimensionality of the input must be {}, not {}"
+                .format(ndim, x.ndim)
+            )
+
+    def _atleast_ndim(self, x, ndim):
+        if x.ndim < ndim:
+            raise ValueError(
+                "dimensionality of the input must be"
+                " larger or equal to {}, not {}"
+                .format(ndim, x.ndim)
+            )
@@ -0,0 +1,8 @@
+from prml.nn.image.convolve2d import convolve2d
+from prml.nn.image.max_pooling2d import max_pooling2d
+
+
+__all__ = [
+    "convolve2d",
+    "max_pooling2d"
+]
@@ -0,0 +1,100 @@
+import numpy as np
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+from prml.nn.image.util import img2patch, patch2img
+
+
+class Convolve2d(Function):
+
+    def __init__(self, stride, pad):
+        """
+        construct 2 dimensional convolution function
+
+        Parameters
+        ----------
+        stride : int or tuple of ints
+            stride of kernel application
+        pad : int or tuple of ints
+            padding image
+        """
+        self.stride = self._check_tuple(stride, "stride")
+        self.pad = self._check_tuple(pad, "pad")
+        self.pad = (0,) + self.pad + (0,)
+
+    def _check_tuple(self, tup, name):
+        if isinstance(tup, int):
+            tup = (tup,) * 2
+        if not isinstance(tup, tuple):
+            raise TypeError(
+                "Unsupported type for {}: {}".format(name, type(tup))
+            )
+        if len(tup) != 2:
+            raise ValueError(
+                "the length of {} must be 2, not {}".format(name, len(tup))
+            )
+        if not all([isinstance(n, int) for n in tup]):
+            raise TypeError(
+                "Unsuported type for {}".format(name)
+            )
+        if not all([n >= 0 for n in tup]):
+            raise ValueError(
+                "{} must be non-negative values".format(name)
+            )
+        return tup
+
+    def _check_input(self, x, y):
+        x = self._convert2tensor(x)
+        y = self._convert2tensor(y)
+        self._equal_ndim(x, 4)
+        self._equal_ndim(y, 4)
+        if x.shape[3] != y.shape[2]:
+            raise ValueError(
+                "shapes {} and {} not aligned: {} (dim 3) != {} (dim 2)"
+                .format(x.shape, y.shape, x.shape[3], y.shape[2])
+            )
+        return x, y
+
+    def forward(self, x, y):
+        x, y = self._check_input(x, y)
+        self.x = x
+        self.y = y
+        img = np.pad(x.value, [(p,) for p in self.pad], "constant")
+        self.shape = img.shape
+        self.patch = img2patch(img, y.shape[:2], self.stride)
+        return Tensor(np.tensordot(self.patch, y.value, axes=((3, 4, 5), (0, 1, 2))), function=self)
+
+    def backward(self, delta):
+        dx = patch2img(
+            np.tensordot(delta, self.y.value, (3, 3)),
+            self.stride,
+            self.shape
+        )
+        slices = [slice(p, len_ - p) for p, len_ in zip(self.pad, self.shape)]
+        dx = dx[slices]
+        dy = np.tensordot(self.patch, delta, axes=((0, 1, 2),) * 2)
+        self.x.backward(dx)
+        self.y.backward(dy)
+
+
+def convolve2d(x, y, stride=1, pad=0):
+    """
+    returns convolution of two tensors
+
+    Parameters
+    ----------
+    x : (n_batch, xlen, ylen, in_channel) Tensor
+        input tensor to be convolved
+    y : (kx, ky, in_channel, out_channel) Tensor
+        convolution kernel
+    stride : int or tuple of ints (sx, sy)
+        stride of kernel application
+    pad : int or tuple of ints (px, py)
+        padding image
+
+    Returns
+    -------
+    output : (n_batch, xlen', ylen', out_channel) Tensor
+        input convolved with kernel
+        len' = (len + p - k) // s + 1
+    """
+    return Convolve2d(stride, pad).forward(x, y)
@@ -0,0 +1,93 @@
+import numpy as np
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+from prml.nn.image.util import img2patch, patch2img
+
+
+class MaxPooling2d(Function):
+
+    def __init__(self, pool_size, stride, pad):
+        """
+        construct 2 dimensional max-pooling function
+
+        Parameters
+        ----------
+        pool_size : int or tuple of ints
+            pooling size
+        stride : int or tuple of ints
+            stride of kernel application
+        pad : int or tuple of ints
+            padding image
+        """
+        self.pool_size = self._check_tuple(pool_size, "pool_size")
+        self.stride = self._check_tuple(stride, "stride")
+        self.pad = self._check_tuple(pad, "pad")
+        self.pad = (0,) + self.pad + (0,)
+
+    def _check_tuple(self, tup, name):
+        if isinstance(tup, int):
+            tup = (tup,) * 2
+        if not isinstance(tup, tuple):
+            raise TypeError(
+                "Unsupported type for {}: {}".format(name, type(tup))
+            )
+        if len(tup) != 2:
+            raise ValueError(
+                "the length of {} must be 2, not {}".format(name, len(tup))
+            )
+        if not all([isinstance(n, int) for n in tup]):
+            raise TypeError(
+                "Unsuported type for {}".format(name)
+            )
+        if not all([n >= 0 for n in tup]):
+            raise ValueError(
+                "{} must be non-negative values".format(name)
+            )
+        return tup
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self._equal_ndim(x, 4)
+        self.x = x
+        img = np.pad(x.value, [(p,) for p in self.pad], "constant")
+        patch = img2patch(img, self.pool_size, self.stride)
+        n_batch, xlen_out, ylen_out, _, _, in_channels = patch.shape
+        patch = patch.reshape(n_batch, xlen_out, ylen_out, -1, in_channels)
+        self.shape = img.shape
+        self.index = patch.argmax(axis=3)
+        return Tensor(patch.max(axis=3), function=self)
+
+    def backward(self, delta):
+        delta_patch = np.zeros(delta.shape + (np.prod(self.pool_size),))
+        index = np.where(delta == delta) + (self.index.ravel(),)
+        delta_patch[index] = delta.ravel()
+        delta_patch = np.reshape(delta_patch, delta.shape + self.pool_size)
+        delta_patch = delta_patch.transpose(0, 1, 2, 4, 5, 3)
+        dx = patch2img(delta_patch, self.stride, self.shape)
+        slices = [slice(p, len_ - p) for p, len_ in zip(self.pad, self.shape)]
+        dx = dx[slices]
+        self.x.backward(dx)
+
+
+def max_pooling2d(x, pool_size, stride=1, pad=0):
+    """
+    spatial max pooling
+
+    Parameters
+    ----------
+    x : (n_batch, xlen, ylen, in_channel) Tensor
+        input tensor
+    pool_size : int or tuple of ints (kx, ky)
+        pooling size
+    stride : int or tuple of ints (sx, sy)
+        stride of pooling application
+    pad : int or tuple of ints (px, py)
+        padding input
+
+    Returns
+    -------
+    output : (n_batch, xlen', ylen', out_channel) Tensor
+        max pooled image
+        len' = (len + p - k) // s + 1
+    """
+    return MaxPooling2d(pool_size, stride, pad).forward(x)
@@ -0,0 +1,62 @@
+import itertools
+import numpy as np
+from numpy.lib.stride_tricks import as_strided
+
+
+def img2patch(img, size, step=1):
+    """
+    convert batch of image array into patches
+    Parameters
+    ----------
+    img : (n_batch, xlen_in, ylen_in, in_channels) ndarray
+        batch of images
+    size : tuple or int
+        patch size
+    step : tuple or int
+        stride of patches
+    Returns
+    -------
+    patch : (n_batch, xlen_out, ylen_out, size, size, in_channels) ndarray
+        batch of patches at all points
+        len_out = (len_in - size) // step + 1
+    """
+    ndim = img.ndim
+    if isinstance(size, int):
+        size = (size,) * (ndim - 2)
+    if isinstance(step, int):
+        step = (step,) * (ndim - 2)
+
+    slices = [slice(None, None, s) for s in step]
+    window_strides = img.strides[1:]
+    index_strides = img[[slice(None)] + slices].strides[:-1]
+
+    out_shape = tuple(
+        np.subtract(img.shape[1: -1], size) // np.array(step) + 1)
+    out_shape = (len(img),) + out_shape + size + (np.size(img, -1),)
+    strides = index_strides + window_strides
+    patch = as_strided(img, shape=out_shape, strides=strides)
+    return patch
+
+
+def patch2img(x, stride, shape):
+    """
+    sum up patches and form an image
+    Parameters
+    ----------
+    x : (n_batch, xlen_in, ylen_in, kx, ky, in_channels) ndarray
+        batch of patches at all points
+    stride : tuple
+        applied stride to take patches
+    shape : (n_batch, xlen_out, ylen_out, in_channels) tuple
+        output image shape
+    Returns
+    -------
+    img : (n_batch, len_out, ylen_out, in_channels) ndarray
+        image
+    """
+    img = np.zeros(shape, dtype=np.float32)
+    kx, ky = x.shape[3: 5]
+    for i, j in itertools.product(range(kx), range(ky)):
+        slices = [slice(b, b + s * len_, s) for b, s, len_ in zip([i, j], stride, x.shape[1: 3])]
+        img[[slice(None)] + slices] += x[..., i, j, :]
+    return img
@@ -0,0 +1,16 @@
+from prml.nn.linalg.cholesky import cholesky
+from prml.nn.linalg.det import det
+from prml.nn.linalg.inv import inv
+from prml.nn.linalg.logdet import logdet
+from prml.nn.linalg.solve import solve
+from prml.nn.linalg.trace import trace
+
+
+__all__ = [
+    "cholesky",
+    "det",
+    "inv",
+    "logdet",
+    "solve",
+    "trace"
+]
@@ -0,0 +1,45 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Cholesky(Function):
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self.output = np.linalg.cholesky(x.value)
+        if isinstance(self.x, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        delta_lower = np.tril(delta)
+        P = phi(self.output.T @ delta_lower)
+        S = np.linalg.solve(
+            self.output.T,
+            P @ np.linalg.inv(self.output)
+        )
+        dx = S + S.T + np.diag(np.diag(S))
+        self.x.backward(dx)
+
+
+def phi(x):
+    return 0.5 * (np.tril(x) + np.tril(x, -1))
+
+
+def cholesky(x):
+    """
+    cholesky decomposition of positive-definite matrix
+    x = LL^T
+    Parameters
+    ----------
+    x : (d, d) tensor_like
+        positive-definite matrix
+    Returns
+    -------
+    L : (d, d)
+        cholesky decomposition
+    """
+    return Cholesky().forward(x)
@@ -0,0 +1,35 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Determinant(Function):
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self._equal_ndim(x, 2)
+        self.output = np.linalg.det(x.value)
+        if isinstance(self.x, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        dx = delta * self.output * np.linalg.inv(self.x.value.T)
+        self.x.backward(dx)
+
+
+def det(x):
+    """
+    determinant of a matrix
+    Parameters
+    ----------
+    x : (d, d) tensor_like
+        a matrix to compute its determinant
+    Returns
+    -------
+    output : (d, d) tensor_like
+        determinant of the input matrix
+    """
+    return Determinant().forward(x)
@@ -0,0 +1,35 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Inverse(Function):
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self._equal_ndim(x, 2)
+        self.output = np.linalg.inv(x.value)
+        if isinstance(self.x, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        dx = -self.output.T @ delta @ self.output.T
+        self.x.backward(dx)
+
+
+def inv(x):
+    """
+    inverse of a matrix
+    Parameters
+    ----------
+    x : (d, d) tensor_like
+        a matrix to be inverted
+    Returns
+    -------
+    output : (d, d) tensor_like
+        inverse of the input
+    """
+    return Inverse().forward(x)
@@ -0,0 +1,37 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class LogDeterminant(Function):
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self._equal_ndim(x, 2)
+        sign, self.output = np.linalg.slogdet(x.value)
+        if sign != 1:
+            raise ValueError("matrix has to be positive-definite")
+        if isinstance(self.x, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        dx = delta * np.linalg.inv(self.x.value.T)
+        self.x.backward(dx)
+
+
+def logdet(x):
+    """
+    log determinant of a matrix
+    Parameters
+    ----------
+    x : (d, d) tensor_like
+        a matrix to compute its log determinant
+    Returns
+    -------
+    output : (d, d) tensor_like
+        determinant of the input matrix
+    """
+    return LogDeterminant().forward(x)
@@ -0,0 +1,43 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Solve(Function):
+
+    def forward(self, a, b):
+        a = self._convert2tensor(a)
+        b = self._convert2tensor(b)
+        self._equal_ndim(a, 2)
+        self._equal_ndim(b, 2)
+        self.a = a
+        self.b = b
+        self.output = np.linalg.solve(a.value, b.value)
+        if isinstance(self.a, Constant) and isinstance(self.b, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        db = np.linalg.solve(self.a.value.T, delta)
+        da = -db @ self.output.T
+        self.a.backward(da)
+        self.b.backward(db)
+
+
+def solve(a, b):
+    """
+    solve a linear matrix equation
+    ax = b
+    Parameters
+    ----------
+    a : (d, d) tensor_like
+        coefficient matrix
+    b : (d, k) tensor_like
+        dependent variable
+    Returns
+    -------
+    output : (d, k) tensor_like
+        solution of the equation
+    """
+    return Solve().forward(a, b)
@@ -0,0 +1,24 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Trace(Function):
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self._equal_ndim(x, 2)
+        self.x = x
+        output = np.trace(x.value)
+        if isinstance(self.x, Constant):
+            return Constant(output)
+        return Tensor(output, function=self)
+
+    def backward(self, delta):
+        dx = np.eye(self.x.shape[0], self.x.shape[1]) * delta
+        self.x.backward(dx)
+
+
+def trace(x):
+    return Trace().forward(x)
@@ -0,0 +1,50 @@
+from prml.nn.math.add import add
+from prml.nn.math.divide import divide, rdivide
+from prml.nn.math.exp import exp
+from prml.nn.math.log import log
+from prml.nn.math.matmul import matmul, rmatmul
+from prml.nn.math.mean import mean
+from prml.nn.math.multiply import multiply
+from prml.nn.math.negative import negative
+from prml.nn.math.power import power, rpower
+from prml.nn.math.product import prod
+from prml.nn.math.sqrt import sqrt
+from prml.nn.math.square import square
+from prml.nn.math.subtract import subtract, rsubtract
+from prml.nn.math.sum import sum
+
+
+from prml.nn.tensor.tensor import Tensor
+Tensor.__add__ = add
+Tensor.__radd__ = add
+Tensor.__truediv__ = divide
+Tensor.__rtruediv__ = rdivide
+Tensor.mean = mean
+Tensor.__matmul__ = matmul
+Tensor.__rmatmul__ = rmatmul
+Tensor.__mul__ = multiply
+Tensor.__rmul__ = multiply
+Tensor.__neg__ = negative
+Tensor.__pow__ = power
+Tensor.__rpow__ = rpower
+Tensor.prod = prod
+Tensor.__sub__ = subtract
+Tensor.__rsub__ = rsubtract
+Tensor.sum = sum
+
+
+__all__ = [
+    "add",
+    "divide",
+    "exp",
+    "log",
+    "matmul",
+    "mean",
+    "multiply",
+    "power",
+    "prod",
+    "sqrt",
+    "square",
+    "subtract",
+    "sum"
+]
@@ -0,0 +1,27 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Abs(Function):
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self.output = np.abs(x.value)
+        if isinstance(x, Constant):
+            return Constant(self.output)
+        self.sign = np.sign(x.value)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        dx = self.sign * delta
+        self.x.backward(dx)
+
+
+def abs(x):
+    """
+    element-wise absolute function
+    """
+    return Abs().forward(x)
@@ -0,0 +1,40 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+from prml.nn.array.broadcast import broadcast_to
+
+
+class Add(Function):
+    """
+    add arguments element-wise
+    """
+
+    def _check_input(self, x, y):
+        x = self._convert2tensor(x)
+        y = self._convert2tensor(y)
+        if x.shape != y.shape:
+            shape = np.broadcast(x.value, y.value).shape
+            if x.shape != shape:
+                x = broadcast_to(x, shape)
+            if y.shape != shape:
+                y = broadcast_to(y, shape)
+        return x, y
+
+    def forward(self, x, y):
+        x, y = self._check_input(x, y)
+        self.x = x
+        self.y = y
+        if isinstance(self.x, Constant) and isinstance(self.y, Constant):
+            return Constant(x.value + y.value)
+        return Tensor(x.value + y.value, function=self)
+
+    def backward(self, delta):
+        dx = delta
+        dy = delta
+        self.x.backward(dx)
+        self.y.backward(dy)
+
+
+def add(x, y):
+    return Add().forward(x, y)
@@ -0,0 +1,44 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+from prml.nn.array.broadcast import broadcast_to
+
+
+class Divide(Function):
+    """
+    divide arguments element-wise
+    """
+
+    def _check_input(self, x, y):
+        x = self._convert2tensor(x)
+        y = self._convert2tensor(y)
+        if x.shape != y.shape:
+            shape = np.broadcast(x.value, y.value).shape
+            if x.shape != shape:
+                x = broadcast_to(x, shape)
+            if y.shape != shape:
+                y = broadcast_to(y, shape)
+        return x, y
+
+    def forward(self, x, y):
+        x, y = self._check_input(x, y)
+        self.x = x
+        self.y = y
+        if isinstance(self.x, Constant) and isinstance(self.y, Constant):
+            return Constant(x.value / y.value)
+        return Tensor(x.value / y.value, function=self)
+
+    def backward(self, delta):
+        dx = delta / self.y.value
+        dy = - delta * self.x.value / self.y.value ** 2
+        self.x.backward(dx)
+        self.y.backward(dy)
+
+
+def divide(x, y):
+    return Divide().forward(x, y)
+
+
+def rdivide(x, y):
+    return Divide().forward(y, x)
@@ -0,0 +1,26 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Exp(Function):
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self.output = np.exp(x.value)
+        if isinstance(self.x, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        dx = self.output * delta
+        self.x.backward(dx)
+
+
+def exp(x):
+    """
+    element-wise exponential function
+    """
+    return Exp().forward(x)
@@ -0,0 +1,26 @@
+import scipy.special as sp
+from prml.nn.function import Function
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+
+
+class Gamma(Function):
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self.output = sp.gamma(x.value)
+        if isinstance(x, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        dx = delta * self.output * sp.digamma(self.x.value)
+        self.x.backward(dx)
+
+
+def gamma(x):
+    """
+    element-wise gamma function
+    """
+    return Gamma().forward(x)
@@ -0,0 +1,31 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Log(Function):
+    """
+    element-wise natural logarithm of the input
+    y = log(x)
+    """
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        output = np.log(self.x.value)
+        if isinstance(self.x, Constant):
+            return Constant(output)
+        return Tensor(output, function=self)
+
+    def backward(self, delta):
+        dx = delta / self.x.value
+        self.x.backward(dx)
+
+
+def log(x):
+    """
+    element-wise natural logarithm of the input
+    y = log(x)
+    """
+    return Log().forward(x)
@@ -0,0 +1,43 @@
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class MatMul(Function):
+    """
+    Matrix multiplication function
+    """
+
+    def _check_input(self, x, y):
+        x = self._convert2tensor(x)
+        y = self._convert2tensor(y)
+        self._equal_ndim(x, 2)
+        self._equal_ndim(y, 2)
+        if x.shape[1] != y.shape[0]:
+            raise ValueError(
+                "shapes {} and {} not aligned: {} (dim 1) != {} (dim 0)"
+                .format(x.shape, y.shape, x.shape[1], y.shape[0])
+            )
+        return x, y
+
+    def forward(self, x, y):
+        x, y = self._check_input(x, y)
+        self.x = x
+        self.y = y
+        if isinstance(self.x, Constant) and isinstance(self.y, Constant):
+            return Constant(x.value @ y.value)
+        return Tensor(x.value @ y.value, function=self)
+
+    def backward(self, delta):
+        dx = delta @ self.y.value.T
+        dy = self.x.value.T @ delta
+        self.x.backward(dx)
+        self.y.backward(dy)
+
+
+def matmul(x, y):
+    return MatMul().forward(x, y)
+
+
+def rmatmul(x, y):
+    return MatMul().forward(y, x)
@@ -0,0 +1,21 @@
+from prml.nn.math.sum import sum
+
+
+def mean(x, axis=None, keepdims=False):
+    """
+    returns arithmetic mean of the elements along given axis
+    """
+    if axis is None:
+        return sum(x, axis=None, keepdims=keepdims) / x.size
+    elif isinstance(axis, int):
+        N = x.shape[axis]
+        return sum(x, axis=axis, keepdims=keepdims) / N
+    elif isinstance(axis, tuple):
+        N = 1
+        for ax in axis:
+            N *= x.shape[ax]
+        return sum(x, axis=axis, keepdims=keepdims) / N
+    else:
+        raise TypeError(
+            "Unsupported type for axis: {}".format(type(axis))
+        )
@@ -0,0 +1,40 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+from prml.nn.array.broadcast import broadcast_to
+
+
+class Multiply(Function):
+    """
+    multiply arguments element-wise
+    """
+
+    def _check_input(self, x, y):
+        x = self._convert2tensor(x)
+        y = self._convert2tensor(y)
+        if x.shape != y.shape:
+            shape = np.broadcast(x.value, y.value).shape
+            if x.shape != shape:
+                x = broadcast_to(x, shape)
+            if y.shape != shape:
+                y = broadcast_to(y, shape)
+        return x, y
+
+    def forward(self, x, y):
+        x, y = self._check_input(x, y)
+        self.x = x
+        self.y = y
+        if isinstance(self.x, Constant) and isinstance(self.y, Constant):
+            return Constant(x.value * y.value)
+        return Tensor(x.value * y.value, function=self)
+
+    def backward(self, delta):
+        dx = self.y.value * delta
+        dy = self.x.value * delta
+        self.x.backward(dx)
+        self.y.backward(dy)
+
+
+def multiply(x, y):
+    return Multiply().forward(x, y)
@@ -0,0 +1,28 @@
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Negative(Function):
+    """
+    element-wise negative
+    y = -x
+    """
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        if isinstance(self.x, Constant):
+            return Constant(-x.value)
+        return Tensor(-x.value, function=self)
+
+    def backward(self, delta):
+        dx = -delta
+        self.x.backward(dx)
+
+
+def negative(x):
+    """
+    element-wise negative
+    """
+    return Negative().forward(x)
@@ -0,0 +1,57 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+from prml.nn.array.broadcast import broadcast_to
+
+
+class Power(Function):
+    """
+    First array elements raised to powers from second array
+    """
+
+    def _check_input(self, x, y):
+        x = self._convert2tensor(x)
+        y = self._convert2tensor(y)
+        if x.shape != y.shape:
+            shape = np.broadcast(x.value, y.value).shape
+            if x.shape != shape:
+                x = broadcast_to(x, shape)
+            if y.shape != shape:
+                y = broadcast_to(y, shape)
+        return x, y
+
+    def forward(self, x, y):
+        x, y = self._check_input(x, y)
+        self.x = x
+        self.y = y
+        self.output = np.power(x.value, y.value)
+        if isinstance(self.x, Constant) and isinstance(self.y, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        dx = self.y.value * np.power(self.x.value, self.y.value - 1) * delta
+        if self.x.size == 1:
+            if self.x.value > 0:
+                dy = self.output * np.log(self.x.value) * delta
+            else:
+                dy = None
+        else:
+            if (self.x.value > 0).all():
+                dy = self.output * np.log(self.x.value) * delta
+            else:
+                dy = None
+        self.x.backward(dx)
+        self.y.backward(dy)
+
+
+def power(x, y):
+    """
+    First array elements raised to powers from second array
+    """
+    return Power().forward(x, y)
+
+
+def rpower(x, y):
+    return Power().forward(y, x)
@@ -0,0 +1,55 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Product(Function):
+
+    def __init__(self, axis=None, keepdims=False):
+        if isinstance(axis, int):
+            axis = (axis,)
+        elif isinstance(axis, tuple):
+            axis = tuple(sorted(axis))
+        self.axis = axis
+        self.keepdims = keepdims
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self.output = np.prod(self.x.value, axis=self.axis, keepdims=True)
+        if not self.keepdims:
+            output = np.squeeze(self.output)
+            if output.size == 1:
+                output = output.item()
+        else:
+            output = self.output
+        if isinstance(self.x, Constant):
+            return Constant(output)
+        return Tensor(output, function=self)
+
+    def backward(self, delta):
+        if not self.keepdims and self.axis is not None:
+            for ax in self.axis:
+                delta = np.expand_dims(delta, ax)
+        dx = delta * self.output / self.x.value
+        self.x.backward(dx)
+
+
+def prod(x, axis=None, keepdims=False):
+    """
+    product of all element in the array
+    Parameters
+    ----------
+    x : tensor_like
+        input array
+    axis : int, tuple of ints
+        axis or axes along which a product is performed
+    keepdims : bool
+        keep dimensionality or not
+    Returns
+    -------
+    product : tensor_like
+        product of all element
+    """
+    return Product(axis=axis, keepdims=keepdims).forward(x)
@@ -0,0 +1,31 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Sqrt(Function):
+    """
+    element-wise square root of the input
+    y = sqrt(x)
+    """
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self.output = np.sqrt(x.value)
+        if isinstance(self.x, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        dx = 0.5 * delta / self.output
+        self.x.backward(dx)
+
+
+def sqrt(x):
+    """
+    element-wise square root of the input
+    y = sqrt(x)
+    """
+    return Sqrt().forward(x)
@@ -0,0 +1,30 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Square(Function):
+    """
+    element-wise square of the input
+    y = x * x
+    """
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        if isinstance(self.x, Constant):
+            return Constant(np.square(x.value))
+        return Tensor(np.square(x.value), function=self)
+
+    def backward(self, delta):
+        dx = 2 * self.x.value * delta
+        self.x.backward(dx)
+
+
+def square(x):
+    """
+    element-wise square of the input
+    y = x * x
+    """
+    return Square().forward(x)
@@ -0,0 +1,44 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+from prml.nn.array.broadcast import broadcast_to
+
+
+class Subtract(Function):
+    """
+    subtract arguments element-wise
+    """
+
+    def _check_input(self, x, y):
+        x = self._convert2tensor(x)
+        y = self._convert2tensor(y)
+        if x.shape != y.shape:
+            shape = np.broadcast(x.value, y.value).shape
+            if x.shape != shape:
+                x = broadcast_to(x, shape)
+            if y.shape != shape:
+                y = broadcast_to(y, shape)
+        return x, y
+
+    def forward(self, x, y):
+        x, y = self._check_input(x, y)
+        self.x = x
+        self.y = y
+        if isinstance(self.x, Constant) and isinstance(self.y, Constant):
+            return Constant(x.value - y.value)
+        return Tensor(x.value - y.value, function=self)
+
+    def backward(self, delta):
+        dx = delta
+        dy = -delta
+        self.x.backward(dx)
+        self.y.backward(dy)
+
+
+def subtract(x, y):
+    return Subtract().forward(x, y)
+
+
+def rsubtract(x, y):
+    return Subtract().forward(y, x)
@@ -0,0 +1,46 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Sum(Function):
+    """
+    summation along given axis
+    y = sum_i=1^N x_i
+    """
+
+    def __init__(self, axis=None, keepdims=False):
+        if isinstance(axis, int):
+            axis = (axis,)
+        self.axis = axis
+        self.keepdims = keepdims
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        output = x.value.sum(axis=self.axis, keepdims=self.keepdims)
+        if isinstance(self.x, Constant):
+            return Constant(output)
+        return Tensor(output, function=self)
+
+    def backward(self, delta):
+        if isinstance(delta, np.ndarray) and (not self.keepdims) and (self.axis is not None):
+            axis_positive = []
+            for axis in self.axis:
+                if axis < 0:
+                    axis_positive.append(self.x.ndim + axis)
+                else:
+                    axis_positive.append(axis)
+            for axis in sorted(axis_positive):
+                delta = np.expand_dims(delta, axis)
+        dx = np.broadcast_to(delta, self.x.shape)
+        self.x.backward(dx)
+
+
+def sum(x, axis=None, keepdims=False):
+    """
+    returns summation of the elements along given axis
+    y = sum_i=1^N x_i
+    """
+    return Sum(axis=axis, keepdims=keepdims).forward(x)
@@ -0,0 +1,90 @@
+from prml.nn.random.random import RandomVariable
+from prml.nn.tensor.parameter import Parameter
+
+
+class Network(object):
+    """
+    a base class for network building
+
+    Parameters
+    ----------
+    kwargs : tensor_like
+        parameters to be optimized
+
+    Attributes
+    ----------
+    parameter : dict
+        dictionary of parameters to be optimized
+    random_variable : dict
+        dictionary of random varibles
+    """
+
+    def __init__(self, **kwargs):
+        self.random_variable = {}
+        self.parameter = {}
+        for key, value in kwargs.items():
+            if isinstance(value, Parameter):
+                self.parameter[key] = value
+            else:
+                try:
+                    value = Parameter(value)
+                except TypeError:
+                    raise TypeError(f"invalid type argument: {type(value)}")
+                self.parameter[key] = value
+            object.__setattr__(self, key, value)
+
+    def __setattr__(self, key, value):
+        if isinstance(value, RandomVariable):
+            self.random_variable[key] = value
+        object.__setattr__(self, key, value)
+
+    def clear(self):
+        """
+        clear gradient and constructed bayesian network
+        """
+        for p in self.parameter.values():
+            p.cleargrad()
+        self.random_variable = {}
+
+    def log_pdf(self, coef=1.):
+        """
+        compute logarithm of probabilty density function
+        Parameters
+        ----------
+        coef : float
+            coefficient to balance likelihood and prior
+            assuming mini-batch size / whole data size for mini-batch training
+        Returns
+        -------
+        logp : tensor_like
+            logarithm of probability density function
+        """
+        logp = 0
+        for rv in self.random_variable.values():
+            if rv.observed:
+                logp += rv.log_pdf().sum()
+            else:
+                logp += coef * rv.log_pdf().sum()
+        return logp
+
+    def elbo(self, coef=1.):
+        """
+        compute evidence lower bound of this model
+        ln p(output) >= elbo
+        Parameters
+        ----------
+        coef : float
+            coefficient to balance likelihood and prior
+            assuming mini-batch size / whole data size for mini-batch training
+        Returns
+        -------
+        evidence : tensor_like
+            evidence lower bound
+        """
+        evidence = 0
+        for rv in self.random_variable.values():
+            if rv.observed:
+                evidence += rv.log_pdf().sum()
+            else:
+                evidence += -coef * rv.KLqp().sum()
+        return evidence
@@ -0,0 +1,32 @@
+import numpy as np
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class LogSoftmax(Function):
+
+    def __init__(self, axis=-1):
+        self.axis = axis
+
+    def _logsumexp(self, x):
+        x_max = np.max(x, axis=self.axis, keepdims=True)
+        y = x - x_max
+        np.exp(y, out=y)
+        np.log(y.sum(axis=self.axis, keepdims=True), out=y)
+        y += x_max
+        return y
+
+    def _forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self.output = x.value - self._logsumexp(x.value)
+        return Tensor(self.output, function=self)
+
+    def _backward(self, delta):
+        dx = delta
+        dx -= np.exp(self.output) * dx.sum(axis=self.axis, keepdims=True)
+        self.x.backward(dx)
+
+
+def log_softmax(x, axis=-1):
+    return LogSoftmax(axis=axis).forward(x)
@@ -0,0 +1,32 @@
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class ReLU(Function):
+    """
+    Rectified Linear Unit
+
+    y = max(x, 0)
+    """
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        output = x.value.clip(min=0)
+        if isinstance(x, Constant):
+            return Constant(output)
+        return Tensor(output, function=self)
+
+    def backward(self, delta):
+        dx = (self.x.value > 0) * delta
+        self.x.backward(dx)
+
+
+def relu(x):
+    """
+    Rectified Linear Unit
+
+    y = max(x, 0)
+    """
+    return ReLU().forward(x)
@@ -0,0 +1,31 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Sigmoid(Function):
+    """
+    logistic sigmoid function
+    y = 1 / (1 + exp(-x))
+    """
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self.output = np.tanh(x.value * 0.5) * 0.5 + 0.5
+        if isinstance(self.x, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        dx = self.output * (1 - self.output) * delta
+        self.x.backward(dx)
+
+
+def sigmoid(x):
+    """
+    logistic sigmoid function
+    y = 1 / (1 + exp(-x))
+    """
+    return Sigmoid().forward(x)
@@ -0,0 +1,39 @@
+import numpy as np
+from prml.nn.tensor.constant import Constant
+from prml.nn.tensor.tensor import Tensor
+from prml.nn.function import Function
+
+
+class Softmax(Function):
+
+    def __init__(self, axis=-1):
+        if not isinstance(axis, int):
+            raise TypeError("axis must be int")
+        self.axis = axis
+
+    def _softmax(self, array):
+        y = array - np.max(array, self.axis, keepdims=True)
+        np.exp(y, out=y)
+        y /= y.sum(self.axis, keepdims=True)
+        return y
+
+    def forward(self, x):
+        x = self._convert2tensor(x)
+        self.x = x
+        self.output = self._softmax(x.value)
+        if isinstance(x, Constant):
+            return Constant(self.output)
+        return Tensor(self.output, function=self)
+
+    def backward(self, delta):
+        dx = self.output * delta
+        dx -= self.output * dx.sum(self.axis, keepdims=True)
+        self.x.backward(dx)
+
+
+def softmax(x, axis=-1):
+    """
+    softmax function along specified axis
+    y_k = exp(x_k) / sum_i(exp(x_i))
+    """
+    return Softmax(axis=axis).forward(x)
--- a/Show More
+++ b/Show More