RNG and RVG for Simulation Engineering (Teknik Simulasi)

RNG-RVG.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm, expon, erlang, laplace, kstest\n",
    "\n",
    "import timeit\n",
    "from datetime import timedelta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RNG_MiddleSquares(n, seed = 4321, digit = 4):\n",
    "    \"\"\"\n",
    "    Generates n random numbers with distribution Uniform(0,1) using middle squares method.\n",
    "    Prints information if the random numbers have repeated.\n",
    "    \"\"\"\n",
    "    if not isinstance(n, int):\n",
    "        raise Exception(\"n must be an integer\")\n",
    "    if n <  0:\n",
    "        raise Exception(\"n must be nonnegative\")\n",
    "    if n == 0 :\n",
    "        return []\n",
    "    if len(str(seed)) > digit:\n",
    "        raise Exception(\"seed cannot be longer than digit\")\n",
    "    \n",
    "    result = []\n",
    "    temp = seed\n",
    "    \n",
    "    begin = int(digit/2)\n",
    "    end = begin + digit\n",
    "    \n",
    "    repeats_at = None\n",
    "    for i in range(n):\n",
    "        square = temp*temp\n",
    "        padded_square = str(square).zfill(2*digit)\n",
    "        middle = int(padded_square[begin:end])\n",
    "        \n",
    "        if middle in result and repeats_at is None:\n",
    "            repeats_at = i + 1\n",
    "        \n",
    "        result.append(middle)\n",
    "        temp = middle\n",
    "        print(i + 1, \"Middle:\", middle, \"Padded squares:\", padded_square)\n",
    "    \n",
    "    if repeats_at is not None:\n",
    "        print(\"RNG repeated itself at iteration\", repeats_at)\n",
    "    else:\n",
    "        print(\"RNG hasn't repeated itself.\")\n",
    "    \n",
    "    result = [i/(10**digit) for i in result]\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 Middle: 6710 Padded squares: 18671041\n",
      "2 Middle: 241 Padded squares: 45024100\n",
      "3 Middle: 580 Padded squares: 00058081\n",
      "4 Middle: 3364 Padded squares: 00336400\n",
      "5 Middle: 3164 Padded squares: 11316496\n",
      "6 Middle: 108 Padded squares: 10010896\n",
      "7 Middle: 116 Padded squares: 00011664\n",
      "8 Middle: 134 Padded squares: 00013456\n",
      "9 Middle: 179 Padded squares: 00017956\n",
      "10 Middle: 320 Padded squares: 00032041\n",
      "RNG hasn't repeated itself.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[0.671, 0.0241, 0.058, 0.3364, 0.3164, 0.0108, 0.0116, 0.0134, 0.0179, 0.032]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RNG_MiddleSquares(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RNG_LCM(n, a = 7**5, c = 0, m = 2**31-1, seed = 123457):\n",
    "    \"\"\"\n",
    "    Generates n random numbers with distribution Uniform(0,1) using Linear Congruential Method.\n",
    "    \"\"\"\n",
    "    if not isinstance(n, int):\n",
    "        raise Exception(\"n must be integer\")\n",
    "    if n <  0:\n",
    "        raise Exception(\"n must be nonnegative\")\n",
    "    if n == 0 :\n",
    "        return []\n",
    "    \n",
    "    result = [(a*seed)+c%m]\n",
    "    for i in range(n-1):\n",
    "        result.append((a*result[-1]+c) % m)\n",
    "    \n",
    "    result = [i/(2**31) for i in result]\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.9662200692109764,\n",
       " 0.26071079075336456,\n",
       " 0.7662622318603098,\n",
       " 0.5693368730135262,\n",
       " 0.8448291937820613,\n",
       " 0.044266507029533386,\n",
       " 0.9871839913539588,\n",
       " 0.601350411772728,\n",
       " 0.8963753702118993,\n",
       " 0.38085416657850146]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RNG_LCM(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RVG_InverseTransform(n):\n",
    "    \"\"\"\n",
    "    Generates n random numbers with distribution Exponential(1)\n",
    "    using Inverse Transform method\n",
    "    \"\"\"\n",
    "    start = timeit.default_timer()\n",
    "    \n",
    "    U = np.random.uniform(size = n)\n",
    "    result = [expon.ppf(i) for i in U]\n",
    "    \n",
    "    end = timeit.default_timer()\n",
    "    print(\"Time elapsed:\", timedelta(seconds = end - start))\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time elapsed: 0:00:00.009182\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[1.2373579392393725,\n",
       " 3.2536210594560564,\n",
       " 3.3370889929753074,\n",
       " 0.2502478616122124,\n",
       " 1.725568133588623,\n",
       " 3.435813250960463,\n",
       " 1.9588112265840711,\n",
       " 0.07342301847320508,\n",
       " 0.4684263100155705,\n",
       " 0.030347852939645835]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RVG_InverseTransform(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time elapsed: 0:00:00.841655\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "KstestResult(statistic=0.033025168962058926, pvalue=0.22056078296859538)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kstest(RVG_InverseTransform(1000), \"expon\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RVG_Convolution(n):\n",
    "    \"\"\"\n",
    "    Generates n random numbers with distribution Erlang(10,1)\n",
    "    by convolution over 10 independent Exponential(1) distributions\n",
    "    \"\"\"\n",
    "    start = timeit.default_timer()\n",
    "    \n",
    "    arr = np.random.exponential(size = (n, 10))\n",
    "    result = np.sum(arr, 1).tolist()\n",
    "    \n",
    "    end = timeit.default_timer()\n",
    "    print(\"Time elapsed:\", timedelta(seconds = end - start))\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time elapsed: 0:00:00.000305\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[13.370118617627085,\n",
       " 14.933038189330732,\n",
       " 5.039015105297263,\n",
       " 8.812694574851518,\n",
       " 9.071790504715352,\n",
       " 8.481998567223117,\n",
       " 10.650739699409161,\n",
       " 5.861713407564656,\n",
       " 8.782132159215605,\n",
       " 8.511690380862975]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RVG_Convolution(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time elapsed: 0:00:00.001029\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "KstestResult(statistic=0.03140020822326928, pvalue=0.27191818598626316)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kstest(RVG_Convolution(1000), 'erlang', args = (10,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RVG_Composition(n):\n",
    "    \"\"\"\n",
    "    Generates n random numbers with distribution Laplace(0,1)\n",
    "    by compositing two back-to-back Exponential(1) distributions\n",
    "    \"\"\"\n",
    "    start = timeit.default_timer()\n",
    "    \n",
    "    result = []\n",
    "    \n",
    "    for i in range(n):\n",
    "        if np.random.uniform() < 0.5:\n",
    "            result.append(-np.log(np.random.uniform()))\n",
    "        else:\n",
    "            result.append(np.log(np.random.uniform()))\n",
    "    \n",
    "    end = timeit.default_timer()\n",
    "    print(\"Time elapsed:\", timedelta(seconds = end - start))\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time elapsed: 0:00:00.000920\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[2.472549607799876,\n",
       " -3.257675005162992,\n",
       " -2.7038705358248922,\n",
       " -0.3559187463309637,\n",
       " -0.3988664081454697,\n",
       " 1.296661596855381,\n",
       " 2.128365294270027,\n",
       " 3.8595992690624925,\n",
       " -1.2657119840418871,\n",
       " -0.6422545909684113]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RVG_Composition(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time elapsed: 0:00:00.035158\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "KstestResult(statistic=0.04089328092765532, pvalue=0.06861893005040165)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kstest(RVG_Composition(1000), \"laplace\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RVG_Acceptance(n):\n",
    "    \"\"\"\n",
    "    Generates n random numbers with distribution Normal(0,1)\n",
    "    by Acceptance-Rejection approach. The result is limited along\n",
    "    (-5,5) interval.\n",
    "    Note: P(X < -5 | X ~ Normal(0,1)) = 2e-7\n",
    "    \"\"\"\n",
    "    start = timeit.default_timer()\n",
    "    \n",
    "    x_lower_limit = -5\n",
    "    x_upper_limit = 5\n",
    "    \n",
    "    Tx = norm.pdf(0)\n",
    "    \n",
    "    result = []\n",
    "    rejected = 0\n",
    "    \n",
    "    while len(result) < n:\n",
    "        Ux = np.random.uniform(x_lower_limit, x_upper_limit)\n",
    "        fx = norm.pdf(Ux)\n",
    "        if np.random.random() <= fx/Tx:\n",
    "            result.append(Ux)\n",
    "        else:\n",
    "            rejected +=1\n",
    "    \n",
    "    end = timeit.default_timer()\n",
    "    print(\"Time elapsed:\", timedelta(seconds = end - start))\n",
    "    print(\"n rejected:\", rejected)\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time elapsed: 0:00:00.033655\n",
      "n rejected: 40\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[-1.220596895843201,\n",
       " 0.6589347766067108,\n",
       " 0.24104797834240443,\n",
       " -0.8276294293413313,\n",
       " 1.3867220437840801,\n",
       " -0.3257231450333622,\n",
       " 0.8466631311650676,\n",
       " -2.062863122145325,\n",
       " 1.6278220305466897,\n",
       " -1.5355890834542563]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RVG_Acceptance(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time elapsed: 0:00:02.533878\n",
      "n rejected: 2942\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "KstestResult(statistic=0.015988002999985262, pvalue=0.956666058216043)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kstest(RVG_Acceptance(1000), 'norm')"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}