{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "
\n", "\n", "# Setting the Random Seed\n", "Author(s): Paul Miles | Date Created: June 19, 2018\n", "\n", "For the purpose of testing the MCMC simulation on a particular problem, it may be useful to check whether the results are repeatable. This can be accomplished by setting the seed for the random number generator used within [pymcmcstat](https://prmiles.wordpress.ncsu.edu/codes/python-packages/pymcmcstat/). This tutorial outlines how to accomplish this, and demonstrates the repeatability of the results." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we import the required packages, define our model/sum-of-squares functions, and define a data set." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from pymcmcstat.MCMC import MCMC\n", "np.seterr(over = 'ignore')\n", "\n", "# define test model function\n", "def test_modelfun(xdata, theta):\n", " m = theta[0]\n", " b = theta[1]\n", " nrow, ncol = xdata.shape\n", " y = np.zeros([nrow, 1])\n", " y[:,0] = m*xdata.reshape(nrow,) + b\n", " return y\n", "\n", "def test_ssfun(theta, data):\n", " xdata = data.xdata[0]\n", " ydata = data.ydata[0]\n", " # eval model\n", " ymodel = test_modelfun(xdata, theta)\n", " # calc sos\n", " ss = sum((ymodel[:, 0] - ydata[:, 0])**2)\n", " return ss\n", "\n", "# define data\n", "nds = 100\n", "x = np.linspace(2, 3, num=nds)\n", "y = 2.*x + 3. + 0.1*np.random.standard_normal(x.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Initialize MCMC object and set seed\n", "By default, no seed is specifield. To specify a seed simply define a numeric value for the keywork argument `rngseed`." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Initialize MCMC object\n", "mcstat = MCMC(rngseed=1234)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Setup and run rest of simulation\n", "We set the\n", "- data\n", "- simulation options\n", "- model settings\n", "- model parameters\n", "\n", "just like we would for any other simulation." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Sampling these parameters:\n", " name start [ min, max] N( mu, sigma^2)\n", " m: 2.00 [ -10.00, inf] N( 0.00e+00, inf)\n", " b: -5.00 [ -10.00, 100.00] N( 0.00e+00, inf)\n", " [-----------------100%-----------------] 5000 of 5000 complete in 1.5 sec" ] } ], "source": [ "mcstat.data.add_data_set(x, y)\n", "mcstat.simulation_options.define_simulation_options(\n", " nsimu=int(5.0e3), updatesigma=True, method='dram')\n", "# update model settings\n", "mcstat.model_settings.define_model_settings(sos_function=test_ssfun)\n", "mcstat.parameters.add_model_parameter(\n", " name='m',\n", " theta0=2.,\n", " minimum=-10,\n", " maximum=np.inf,\n", " sample=True)\n", "mcstat.parameters.add_model_parameter(\n", " name='b',\n", " theta0=-5.,\n", " minimum=-10,\n", " maximum=100,\n", " sample=True)\n", "# run mcmc\n", "mcstat.run_simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Extract Results and Display Chainstats\n", "In addition, we check the last row of the chain to see where the simulation ended." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "---------------------\n", "name : mean std MC_err tau geweke\n", "m : 2.0130 0.0320 0.0038 48.7245 0.9959\n", "b : 2.9673 0.0800 0.0096 49.0697 0.9921\n", "---------------------\n", "chain[-1, :] = [2.01337637 2.95431498]\n" ] } ], "source": [ "# Extract results\n", "results = mcstat.simulation_results.results\n", "chain = results['chain']\n", "s2chain = results['s2chain']\n", "sschain = results['sschain']\n", "names = results['names']\n", "# define burnin\n", "burnin = 2000\n", "# display chain statistics\n", "mcstat.chainstats(chain[burnin:, :], results)\n", "print('chain[-1, :] = {}'.format(chain[-1, :]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Check Repeatability\n", "To check the repeatability we simply create a new MCMC object with the same random seed and compare the result." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Sampling these parameters:\n", " name start [ min, max] N( mu, sigma^2)\n", " m: 2.00 [ -10.00, inf] N( 0.00e+00, inf)\n", " b: -5.00 [ -10.00, 100.00] N( 0.00e+00, inf)\n", " [-----------------100%-----------------] 5000 of 5000 complete in 1.4 sec\n", "---------------------\n", "name : mean std MC_err tau geweke\n", "m : 2.0130 0.0320 0.0038 48.7245 0.9959\n", "b : 2.9673 0.0800 0.0096 49.0697 0.9921\n", "---------------------\n", "chain[-1, :] = [2.01337637 2.95431498]\n" ] } ], "source": [ "# Initialize MCMC object\n", "mcstat = MCMC(rngseed=1234)\n", "mcstat.data.add_data_set(x, y)\n", "mcstat.simulation_options.define_simulation_options(\n", " nsimu=int(5.0e3), updatesigma=True, method='dram')\n", "# update model settings\n", "mcstat.model_settings.define_model_settings(sos_function=test_ssfun)\n", "mcstat.parameters.add_model_parameter(\n", " name='m',\n", " theta0=2.,\n", " minimum=-10,\n", " maximum=np.inf,\n", " sample=True)\n", "mcstat.parameters.add_model_parameter(\n", " name='b',\n", " theta0=-5.,\n", " minimum=-10,\n", " maximum=100,\n", " sample=True)\n", "# run mcmc\n", "mcstat.run_simulation()\n", "# Extract results\n", "results = mcstat.simulation_results.results\n", "chain = results['chain']\n", "s2chain = results['s2chain']\n", "sschain = results['sschain']\n", "names = results['names']\n", "# define burnin\n", "burnin = 2000\n", "# display chain statistics\n", "mcstat.chainstats(chain[burnin:, :], results)\n", "print('chain[-1, :] = {}'.format(chain[-1, :]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is clearly seen that the results are identical to the first simulation, so the random number process is repeatable." ] } ], "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.6.8" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": true, "user_envs_cfg": false }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }