Chapter 1. Monte Carlo Methods.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d02c17d0",
   "metadata": {},
   "source": [
    "**Essential imports:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b88afb0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random as r\n",
    "import numpy as n\n",
    "import math as m\n",
    "from numpy.linalg import eig\n",
    "import matplotlib.pyplot as p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d65839f",
   "metadata": {},
   "source": [
    "**Calculation of $\\pi$ using direct sampling:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0dac259f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Value of pi = 3.1627\n"
     ]
    }
   ],
   "source": [
    "# Assume the radius = 1, and the center of the cirlce is the origin of the Cartesian coordinates.\n",
    "\n",
    "n_trials = 40000\n",
    "n_hits = 0\n",
    "\n",
    "for i in range(n_trials):\n",
    "    x, y = r.uniform(-1.0, 1.0), r.uniform(-1, 1)\n",
    "    if x**2 + y**2 < 1: n_hits += 1 # r^2 = 1\n",
    "        \n",
    "print(\"Value of pi =\", 4*n_hits/n_trials)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ecb0d9d",
   "metadata": {},
   "source": [
    "**Calculation of $\\pi$ using direct sampling - multi-runs:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4fe3ddf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def direct_pi(N):\n",
    "    n_hits = 0\n",
    "    for i in range(N):\n",
    "        x, y = r.uniform(-1, 1), r.uniform(-1, 1)\n",
    "        if x**2 + y**2 < 1: n_hits += 1\n",
    "    return n_hits\n",
    "\n",
    "n_runs = 1000\n",
    "n_trials = 4000\n",
    "\n",
    "for run in range(n_runs): obs = direct_pi(n_trials)/n_trials"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05d778d9",
   "metadata": {},
   "source": [
    "**Calculation of $\\pi$ using Markov chain sampling:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cae9303",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = 1, 1 # Starting at the upper right corner of the square and assuming that the center of the cicle is the origin\n",
    "delta = 0.1 # There are lower and upper limits to the extent of the throw.\n",
    "\n",
    "# Does the convergence of the markov-chain Monte Carlo algorithm depend on the starting point and the value of delta?\n",
    "\n",
    "n_trials = int(input(\"Number of trials: \"))\n",
    "n_hits = 0\n",
    "for i in range(n_trials):\n",
    "    del_x, del_y = r.uniform(-delta, delta), r.uniform(-delta, delta)\n",
    "    if abs(x + del_x) < 1 and abs(y + del_y) < 1: # That is, if the pebble is thrown inside the square... Note that the trial is counted anyway.\n",
    "        x, y = x + del_x, y + del_y\n",
    "    if x**2+y**2 < 1: n_hits += 1 # That is, if the pebble is inside the circle... Rejections inside the circle are still counted as hits.\n",
    "print(\"Value of pi:\", 4*n_hits/n_trials)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19caa4ba",
   "metadata": {},
   "source": [
    "**Calculation of $\\pi$ using Markov chain sampling - multi-runs:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e75eb554",
   "metadata": {},
   "outputs": [],
   "source": [
    "def markov_pi(N, delta):\n",
    "    x, y = 1, 1\n",
    "    n_hits = 0\n",
    "    for i in range(N):\n",
    "        del_x, del_y = r.uniform(-delta, delta), r.uniform(-delta, delta)\n",
    "        if abs(x + del_x) < 1 and abs(y + del_y) < 1:\n",
    "            x, y = x + del_x, y + del_y\n",
    "        if x**2 + y**2 < 1.0: n_hits += 1\n",
    "    return n_hits\n",
    "\n",
    "n_runs = int(input(\"Num of runs = \"))\n",
    "n_trials = int(input(\"Num of trials = \"))\n",
    "delta = float(input(\"Value of delta = \"))\n",
    "for run in range(n_runs): print(4*markov_pi(n_trials, delta)/n_trials)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78ee6986",
   "metadata": {},
   "source": [
    "**Pebble game, discrete Markov-chain sampling and the use of the neighbor table:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd014654",
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbor = [[1, 3, 0, 0], # The current site is just the raw number.\n",
    "            [2, 4, 0, 1],\n",
    "            [2, 5, 1, 2],\n",
    "            [4, 6, 3, 0],\n",
    "            [5, 7, 3, 1],\n",
    "            [5, 8, 4, 2],\n",
    "            [7, 6, 6, 3],\n",
    "            [8, 7, 6, 4],\n",
    "            [8, 8, 7, 5]]\n",
    "\n",
    "# Numbering scheme of the pepple game box:\n",
    "\n",
    "#   6   7   8\n",
    "#   3   4   5\n",
    "#   0   1   2\n",
    "\n",
    "t_max = int(input(\"Number of desired moves = \"))\n",
    "site = int(input(\"Current site = \"))\n",
    "t = 0 # Actually, this is discretized time!\n",
    "while t < t_max:\n",
    "    t = t + 1\n",
    "    site = neighbor[site][r.randint(0, 3)] # The next site is chosen randomly. randint(a, b) generates a random integer bet. a & b inclusive.\n",
    "    print(site)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e42fcb90",
   "metadata": {},
   "source": [
    "**Pebble game with inhomogenous transition probability:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2617f3bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "hist = [0, 0, 0, 0, 0, 0, 0, 0, 0]\n",
    "\n",
    "neighbor =  [[1, 3, 0, 0],\n",
    "             [2, 4, 0, 1],\n",
    "             [2, 5, 1, 2],\n",
    "             [4, 6, 3, 0],\n",
    "             [5, 7, 3, 1],\n",
    "             [5, 8, 4, 2],\n",
    "             [7, 6, 6, 3],\n",
    "             [8, 7, 6, 4],\n",
    "             [8, 8, 7, 5]]\n",
    "\n",
    "weight = [3.0, 0.5, 1.0, 0.5, 1.0, 0.5, 2.0, 0.5, 1.0] # Statistical weights of configs, not probs.\n",
    "pos = 8\n",
    "n_iterations = int(input(\"Number of iterations: \"))\n",
    "\n",
    "for i in range(n_iterations):\n",
    "    new_pos = neighbor[pos][r.randint(0, 3)]\n",
    "    if r.random() < weight[new_pos]/weight[pos]: pos = new_pos # MOST IMPORTANT LINE IN THE CODE!\n",
    "    hist[pos] += 1\n",
    "\n",
    "norm = sum(weight)\n",
    "for k in range(9): print('Site:\\t', k, '\\tWeight:\\t', weight[k], '\\tHistogram:\\t', norm*hist[k]/n_iterations)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d388c03",
   "metadata": {},
   "source": [
    "**Performing data bunching on a Markov-chain generated dataset (experimental code):**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cd33565",
   "metadata": {},
   "outputs": [],
   "source": [
    "def markov_pi(N, d):\n",
    "    x, y = 1, 1\n",
    "    data = []\n",
    "    for i in range(N):\n",
    "        dx, dy = r.uniform(-d, d), r.uniform(-d, d)\n",
    "        if abs(x + dx) < 1 and abs(y + dy) < 1: x, y = x + dx, y + dy\n",
    "        if x**2+y**2 < 1: data.append(4)\n",
    "        else: data.append(0)\n",
    "    return data\n",
    "\n",
    "n_trials = 1000\n",
    "delta = 0.1\n",
    "data = markov_pi(n_trials, delta)\n",
    "\n",
    "sum1 = 0\n",
    "sum2 = 0\n",
    "N = int(len(data)/2)\n",
    "y = []\n",
    "for i in range(N):\n",
    "    # print(\"data[2i]\\t\", str('% 0.4f' % data[2*i]).replace(\"'\", \"\"), \"\\tdata[2i+1]\\t\", str('% 0.4f' % data[2*i+1]).replace(\"'\", \"\"))\n",
    "    sum1 += data[2*i] + data[2*i+1]\n",
    "    sum2 += data[2*i]**2 + data[2*i+1]**2\n",
    "    y.append((data[2*i] + data[2*i+1])/2)\n",
    "\n",
    "error = m.sqrt(sum2/(2*N) - (sum1/(2*N))**2)/m.sqrt(2*N)\n",
    "print(\"Sum/2N:\", \"% 0.6f\" % (sum1/(2*N)), \"\\tError:\", \"% 0.6f\" % error)\n",
    "# print(\"\\n\"+str([\"% 0.4f\" % i for i in y]).replace(\"'\", \"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5d1fdef",
   "metadata": {},
   "source": [
    "**Generating the transfer matrix from the neighbor table, and then generating the probability vectors from the transfer matrix:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b06d62fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbor =  [[1, 3, 0, 0],\n",
    "             [2, 4, 0, 1],\n",
    "             [2, 5, 1, 2],\n",
    "             [4, 6, 3, 0],\n",
    "             [5, 7, 3, 1],\n",
    "             [5, 8, 4, 2],\n",
    "             [7, 6, 6, 3],\n",
    "             [8, 7, 6, 4],\n",
    "             [8, 8, 7, 5]]\n",
    "\n",
    "transfer = n.zeros((9, 9))\n",
    "\n",
    "# Filling the transfer matrix:\n",
    "for i in range(9):\n",
    "    for j in range(4): transfer[neighbor[i][j], i] += 0.25\n",
    "\n",
    "n.set_printoptions(formatter={'float': '{:0.2f}'.format})\n",
    "\n",
    "print(\"Transfer matrix:\\n\")\n",
    "for row in range(transfer.shape[0]): print(transfer[row, :])\n",
    "\n",
    "print(\"\\nProbabiliy vectors:\\n\")\n",
    "\n",
    "n.set_printoptions(formatter={'float': '{:0.5f}'.format})\n",
    "\n",
    "# Constructing the probability vectors:\n",
    "prob_vec = n.zeros(9)\n",
    "prob_vec[8] = 1.0\n",
    "\n",
    "for t in range(100):\n",
    "    print(t, \"\\t\", prob_vec, \"\\n\")\n",
    "    prob_vec = n.matmul(transfer, prob_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de311d61",
   "metadata": {},
   "source": [
    "**Finding eigenvalues and eigenvectors of the transfer matrix**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11a32ef1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.linalg import eig\n",
    "\n",
    "eigenvalues, eigenvectors = eig(transfer)\n",
    "n.set_printoptions(formatter={'float': '{: 0.4f}'.format})\n",
    "for i in range(9): print(\"% 0.2f\" % (eigenvalues[i]), eigenvectors[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc1e293b",
   "metadata": {},
   "source": [
    "**Escaping reducibility in Markov-chain Monte Carlo**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c81ded7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "transfer = n.zeros((18, 18))\n",
    "\n",
    "for i in range(9):\n",
    "    for j in range(4): \n",
    "        transfer[neighbor[i][j], i] += 0.25     # Red pebble game\n",
    "        transfer[neighbor[i][j]+9, i+9] += 0.25 # Blue pebble game   \n",
    "\n",
    "# Small transition epsilon between red 2 and blue 6\n",
    "epsilon = 0.04\n",
    "\n",
    "transfer[6+9, 2] = epsilon\n",
    "transfer[2, 6+9] = epsilon\n",
    "\n",
    "# This prob. has to be taken from the not-moving prob. so that prob. sum = 1.\n",
    "transfer[2, 2] -= epsilon    \n",
    "transfer[6+9, 6+9] -= epsilon\n",
    "\n",
    "eigenvalues, eigenvectors = eig(transfer)\n",
    "\n",
    "n.set_printoptions(formatter={'complexfloat': '{: .4f}'.format})\n",
    "for i in eigenvalues: print('{: .4f}'.format(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34f64f6f",
   "metadata": {},
   "source": [
    "**Escaping periodicity (i.e., recurrence) in Markov-chain Monte Carlo**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7215e462",
   "metadata": {},
   "outputs": [],
   "source": [
    "eps = 0.01\n",
    "transfer = n.array([[eps,   1-eps],\n",
    "                    [1-eps, eps  ]])\n",
    "eigenvalues, eigenvectors = eig(transfer)\n",
    "for i in eigenvalues: print('{: .4f}'.format(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4506df77",
   "metadata": {},
   "source": [
    "**The following code is wrong because it is sampling two random numbers:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "262d28b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 20\n",
    "position = 0\n",
    "for t in range(10):\n",
    "    if r.uniform(0, 1) < 0.5: position = (position + 1) % N\n",
    "    elif r.uniform(0, 1) > 0.5: position = (position - 1) % N\n",
    "    print(position)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "077ad610",
   "metadata": {},
   "source": [
    "**In direct sampling, the standard deviation decreases as the number of trials is increased:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d61e3ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "def direct_pi(N):\n",
    "    n_hits = 0\n",
    "    for i in range(N):\n",
    "        x, y = r.uniform(-1, 1), r.uniform(-1, 1)\n",
    "        if x**2+y**2 < 1.0:\n",
    "            n_hits += 1\n",
    "    return n_hits\n",
    "\n",
    "n_runs = 500\n",
    "n_trials_list = []\n",
    "sigmasq_list = []\n",
    "\n",
    "for i in range(4, 13):\n",
    "    n_trials = 2**i\n",
    "    sigmasq = 0.0\n",
    "    for run in range(n_runs):\n",
    "        pi_est = 4*direct_pi(n_trials)/n_trials\n",
    "        sigmasq += (pi_est - m.pi)**2\n",
    "    sigmasq_list.append(m.sqrt(sigmasq/(n_runs)))\n",
    "    n_trials_list.append(n_trials)\n",
    "\n",
    "p.plot(n_trials_list, sigmasq_list, 'o')\n",
    "p.xscale('log')\n",
    "p.yscale('log')\n",
    "p.xlabel('Number of trials')\n",
    "p.ylabel('Root mean square deviation')\n",
    "p.title('Direct sampling of $\\pi$: root mean square deviation vs. number of trials')\n",
    "p.savefig('direct_sampling_rms_deviation.png')\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ada8e6be",
   "metadata": {},
   "source": [
    "**Choosing $\\delta$ that gives an acceptance rate approx. equal to 0.5:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b637da7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = 1.0, 1.0\n",
    "delta = [0.062, 0.125, 0.25, 0.5, 1.0, 2.0, 4.0]\n",
    "\n",
    "n_trials = 2**12\n",
    "n_hits = 0\n",
    "accept_rate = 0\n",
    "for d in delta:\n",
    "    for i in range(n_trials):\n",
    "        del_x, del_y = r.uniform(-d, d), r.uniform(-d, d)\n",
    "        if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:\n",
    "            x, y = x + del_x, y + del_y\n",
    "            accept_rate += 1\n",
    "        if x**2+y**2 < 1.0: n_hits += 1\n",
    "    print(\"delta\\t\", d,\"\\tAcceptance rate:\\t\", n_hits/n_trials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41788d2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = 1.0, 1.0\n",
    "delta = [0.062, 0.125, 0.25, 0.5, 1.0, 2.0, 4.0]\n",
    "\n",
    "n_trials = 2**12\n",
    "n_hits = 0\n",
    "\n",
    "print(\"delta\\t|\\t Acceptance rate\")\n",
    "for d in delta:\n",
    "    accept_rate = 0\n",
    "    for i in range(n_trials):\n",
    "        del_x, del_y = r.uniform(-d, d), r.uniform(-d, d)\n",
    "        if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:\n",
    "            x, y = x + del_x, y + del_y\n",
    "            accept_rate += 1\n",
    "    print(d, \"\\t|\\t\", accept_rate/n_trials)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aecc517",
   "metadata": {},
   "source": [
    "**Defining `obs` to generate a 1-D Markov-chain generated dataset that is suitable for data bunching:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2613072",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_trials = 400000\n",
    "n_hits = 0\n",
    "obs_sum = 0\n",
    "obs2_sum = 0\n",
    "\n",
    "for i in range(n_trials):\n",
    "    x, y = r.uniform(-1, 1), r.uniform(-1, 1)\n",
    "    obs = 0\n",
    "    if x**2+y**2 < 1:\n",
    "        n_hits += 1\n",
    "        obs = 4 # \"obs\" is a Bernoulli variable that takes either 0 or 4.\n",
    "    obs_sum += obs\n",
    "    obs2_sum += obs**2\n",
    "\n",
    "obs_av = obs_sum/n_trials\n",
    "obs2_av = obs2_sum/n_trials\n",
    "var = obs2_av - obs_av**2\n",
    "sdev = m.sqrt(var)\n",
    "\n",
    "print('Obs av:\\t\\t{0:.4f}\\nObs^2 av:\\t{1:.4f}\\nVariance:\\t{2:.4f}\\nStandard dev:\\t{3:2.4f}'.format(obs_av, obs2_av, var, sdev))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c3a5a10",
   "metadata": {},
   "source": [
    "**Data bunching for a Markov chain:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75d6ef4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generating 1-D dataset for a 2-D problem: \n",
    "\n",
    "def markov_pi(N, d):\n",
    "    x, y = 1, 1\n",
    "    data = []\n",
    "    for i in range(N):\n",
    "        dx, dy = r.uniform(-d, d), r.uniform(-d, d)\n",
    "        if abs(x + dx) < 1 and abs(y + dy) < 1: x, y = x + dx, y + dy\n",
    "        if x**2+y**2 < 1: data.append(4)\n",
    "        else: data.append(0)\n",
    "    return data\n",
    "\n",
    "k = 20\n",
    "n_trials = 2**k\n",
    "delta = 0.1\n",
    "data = markov_pi(n_trials, delta)\n",
    "errors  = []\n",
    "bunches = []\n",
    "\n",
    "for i in range(k):\n",
    "    new_data = []\n",
    "    mean = 0\n",
    "    mean_sq = 0\n",
    "    N = len(data)\n",
    "    while data != []:\n",
    "        x = data.pop()\n",
    "        y = data.pop()\n",
    "        mean += x + y\n",
    "        mean_sq += x**2 + y**2\n",
    "        new_data.append((x + y)/2)\n",
    "    errors.append(m.sqrt(mean_sq/N - (mean/N)**2)/m.sqrt(N))\n",
    "    bunches.append(i)\n",
    "    data = new_data[:]\n",
    "    print(abs(m.pi - mean/N))\n",
    "\n",
    "p.plot(bunches, errors, '--s')\n",
    "p.xlabel('Iteration')\n",
    "p.ylabel('Apparent error')\n",
    "p.title('Bunching: naive error vs. iteration number')\n",
    "p.grid(\"on\")\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "740a3a47",
   "metadata": {},
   "source": [
    "**Bayesian Monte Carlo:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec01f661",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_hits = 0\n",
    "\n",
    "for k in range(10):\n",
    "    while n_hits != 3156: \n",
    "        pi_test = r.uniform(0, 4)\n",
    "        n_hits = 0\n",
    "        for i in range(4000):\n",
    "            if r.uniform(0, 1) < pi_test/4: n_hits += 1\n",
    "    print(pi_test)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}