{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Programming Assignment 2: Wiener and Poisson Processes"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"print \"Modules Imported!\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Wiener Process"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will now use Python to construct, plot, and manipulate (discrete approximations to) sample paths of Wiener processes. In the lectures, we have constructed the Wiener process as a continuous-time limit of random walks: Let $(U_t)_{t \\in {\\mathbb Z}_+}$ be an i.i.d. sequence of Rademacher random variables, that is, ${\\mathbf P}[U_t = \\pm 1] = \\frac{1}{2}$ for all $t$. Fix two positive parameters $h > 0$ and $\\tau > 0$ and define a continuous-time stochastic signal $W = (W_t)_{t \\ge 0}$ with the initial condition $W_0 = 0$ as a linear interpolation of the random points \n",
"\n",
"$$\n",
"W_{n\\tau} = h (U_0 + U_1 + \\ldots + U_{n-1}), \\qquad n = 1,2,3,\\ldots.\n",
"$$\n",
"\n",
"The Wiener process with diffusion parameter $D$ arises when we take the continuous-time limit $\\tau \\to 0$, while ensuring that $h^2/\\tau = D$.\n",
"\n",
"In this assignment, we will use an alternative construction of the Wiener process, using standard normal increments. To that end, let $(Z_t)_{t \\in {\\mathbb Z}_+}$ be a sequence of i.i.d. $N(0,1)$ random variables. Fix a small parameter $h > 0$ and define a continuous-time stochastic signal $W = (W_t)_{t \\ge 0}$ with zero initial condition $W_0 = 0$ as a linear interpolation of the random points\n",
"\n",
"$$\n",
"W_{nh} = \\sqrt{Dh} (Z_0 + Z_1 + \\ldots Z_{n-1}), \\qquad n = 1,2,3,\\ldots.\n",
"$$\n",
"\n",
"The code below constructs such a signal on a given interval $[0,t_{{\\max}}]$ for given values of $D$ and $h$ and plots it."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# define WienerProcess: this function will produce random output\n",
"# D: diffusion constant (set D=1 for standard Wiener process)\n",
"# h: timescale\n",
"# tmax: length of time interval\n",
"def WienerProcess(D,h,tmax):\n",
" z = np.random.normal(size = int(tmax/h)-1) # generate a sequence of iid N(0,1) random increments\n",
" zsums = np.append([0],np.cumsum(np.sqrt(D*h) * z)) # compute scaled running sums of the increments\n",
" tvals = np.linspace(0,tmax,tmax*int(1/h)) # define the interpolation times\n",
" t = np.linspace(0,tmax,10*int(1/h))\n",
" w = np.interp(t,tvals,zsums) # construct linear interpolation of w\n",
" return (t,w) \n",
"\n",
"# generate four successively better approximations of standard Wiener process\n",
"(t1,W1) = WienerProcess(1,0.1,5)\n",
"(t2,W2) = WienerProcess(1,0.01,5)\n",
"(t3,W3) = WienerProcess(1,0.001,5)\n",
"(t4,W4) = WienerProcess(1,0.0001,5)\n",
"\n",
"# plot these approximations\n",
"# use plt.figure() to initialize a new plot\n",
"plt.figure() \n",
"plt.plot(t1,W1)\n",
"plt.title('Wiener Process (h=0.1)')\n",
"plt.xlabel('Time')\n",
"plt.ylabel('Position')\n",
"\n",
"plt.figure()\n",
"plt.plot(t2,W2)\n",
"plt.title('Wiener Process (h=0.01)')\n",
"plt.xlabel('Time')\n",
"plt.ylabel('Position')\n",
"\n",
"plt.figure()\n",
"plt.plot(t3,W3)\n",
"plt.title('Wiener Process (h=0.01)')\n",
"plt.xlabel('Time')\n",
"plt.ylabel('Position')\n",
"\n",
"plt.figure()\n",
"plt.plot(t4,W4)\n",
"plt.title('Wiener Process (h=0.001)')\n",
"plt.xlabel('Time')\n",
"plt.ylabel('Position')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Problem 1. By now you know that if $W = (W_t)_{t \\ge 0}$ is a Wiener process with diffusion constant $D$, its time-rescaled version $V = (V_t)_{t \\ge 0}$ with \n",
"\n",
"$$\n",
"V_t = \\frac{1}{\\sqrt{c}} W_{ct}\n",
"$$\n",
"\n",
"for some $c > 0$ is also a Wiener process with the same diffusion constant. \n",
"\n",
"Use the code cell below to generate sample paths of time-rescaled Wiener processes. Your function should take the following inputs: diffusion constant $D$, time scale $h$, maximum time $t_{{\\max}}$, rescaling parameter $c$. Generate and plot sample paths of time-rescaled Wiener processes with $D = 2$, $h=0.001$, $t_{{\\max}} = 5$ and $c = 0.01, 0.1, 1, 10, 100$. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"###### student code for Problem 1 goes here ######"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Problem 2. As a sanity check, modify the above code for WienerProcess to implement the construction from the lecture notes, where the random increments $U_0,U_1,\\ldots$ are i.i.d. Rademacher random variables, i.e., ${\\bf P}[U_k = \\pm 1] = \\frac{1}{2}$. To generate the Rademacher samples, you will need the following Python code that outputs an array containing $n$ i.i.d.\\ ${\\rm Bern}(p)$ ranodm variables:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"b = np.random.binomial(1,p,size=n) # p is the bias between 0 and 1; n is the desired number of samples"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the code cell below to write your modified code and to plot four sample paths with the same parameters as above. You should be seeing roughly the same behavior as in the case of i.i.d. standard normal increments."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"###### student code for Problem 2 goes here ######"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Poisson Process"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our second set of problems concerns the Poisson process. Recall that we can construct a Poisson process $N = (N_t)_{t \\ge 0}$ with arrival rate $\\lambda$ as follows. Let $Z_0,Z_1,\\ldots$ be a sequence of i.i.d. ${\\rm Exp}(\\lambda)$ random variables. Define the random arrival times $T_1,T_2,\\ldots$ by\n",
"\n",
"$$\n",
"T_i = Z_0 + \\ldots + Z_{i-1}, \\qquad i = 1,2,\\ldots\n",
"$$\n",
"\n",
"and for every $t \\ge 0$ take\n",
"\n",
"$$\n",
"N_t = \\sum^\\infty_{i=1} u(t-T_i),\n",
"$$\n",
"\n",
"where $u(\\cdot)$ is the unit step function.\n",
"\n",
"The Python code below generates a sample path of a Poisson process with a given rate up to a given number of arrivals."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# define UnitStep\n",
"def UnitStep(t): # the input to UnitStep is an array of reals\n",
" return np.where(t<0, 0., 1.) # returns an array of 0's (where entries of t are negative) and 1's (otherwise)\n",
"\n",
"# define PoissonProcess: this function will produce random output\n",
"# rate: arrival rate (mean number of arrivals per unit time)\n",
"# nmax: total number of arrivals\n",
"def PoissonProcess(rate,nmax):\n",
" z = np.random.exponential(1/rate,size=nmax) # generate i.i.d. sequence of interarrival times\n",
" arrivals = np.cumsum(z) # compute arrival times T1, ..., Tn\n",
" t = np.linspace(0,int(arrivals[-1])+1,1000) # define time interval from t=0 to t=T_n\n",
" n = [np.sum([UnitStep(time-arrival) for arrival in arrivals]) for time in t] \n",
" return (t,n)\n",
"\n",
"# generate a sample path of unit Poisson process with 10 arrivals and plot it\n",
"(a,N) = PoissonProcess(1,10)\n",
"\n",
"plt.figure()\n",
"plt.plot(a,N)\n",
"plt.title('Poisson process (lambda = 1)')\n",
"plt.xlabel('Time')\n",
"plt.ylabel('Number of arrivals')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Problem 3: The code above produces a segment of a Poisson process on the random interval $[0,T_n]$, where $T_n$ is the $n$th arrival time for a given value of $n$, whereas ourr code for the Wiener process generated sample paths on a given interval $[0,t_{{\\max}}]$. \n",
"\n",
"In the code cell below, construct a new function that would generate a sample path of a Poisson process with a given rate on a given interval $[0,t_{\\max}]$. To do this, use the discrete-time approximation discussed in class, where you chop the interval $[0,t_{\\max}]$ into $n$ subintervals, for a very large $n$, and then let the arrivals be represented by $n$ i.i.d. Bernoulli random variables with bias $p = \\lambda t_{\\max}/n$. Your function should take the following inputs: the arrival rate $\\lambda$ and the maximum time $t_{{\\max}}$. You will need to set the number of subintervals $n$ large enough, so that $\\lambda t_{\\max} < n$. You may find the UnitStep function useful. \n",
"\n",
"Plot the resulting sample paths for $\\lambda = 1,2,4,10$ and for $t_{\\max} = 5$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"###### student code for Problem 3 goes here ######"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Problem 4: Poisson processes are used to model situations where discrete events happen at random times. For example, a Poisson process with rate $\\lambda$ can be used to model the number of customers arriving at a ticket counter in the airport, where $\\lambda$ is the average number of new customer arrivals per unit time. \n",
"\n",
"In this problem, we will consider the situation when there are several independent queues of customers arriving at the counter. Formally, let $m$ be the number of queues, and for each $k \\in \\{ 1,\\ldots,m \\}$ let $N^{(k)} = (N^{(k)}_t)_{t \\ge 0}$ be a Poisson process with rate $\\lambda_k$. Thus, $\\lambda_k$ is the average number of customers arrivals per unit time via the $k$th queue. We assume that these Poisson processes are mutually independent. The total number of arrivals at the counter at time $t$ is then\n",
"\n",
"$$\n",
"N_t = \\sum^m_{k=1} N^{(k)}_t.\n",
"$$\n",
"\n",
"In the cell below, adapt your code from Problem 3 to generate two independent Poisson processes with rates $\\lambda_1 = 3$ and $\\lambda_2 = 5$ on the interval $[0,10]$ and plot the sample path of their sum."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"###### student code for Problem 4 goes here ######"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the markdown cell below, comment on what you see. Do you think that the sum process $N$ is also a Poisson process? If so, why do you think so, and what do you think its arrival rate is? \n",
"\n",
"Hint: Consider two independent Poisson random variables $X \\sim {\\rm Pois}(\\lambda_1)$ and $Y \\sim {\\rm Pois}(\\lambda_2)$. Prove that their sum $X + Y$ is also a Poisson random variable with parameter $\\lambda_1 + \\lambda_2$. Now let $N^{(1)}$ and $N^{(2)}$ be two independent Poisson processes with rates $\\lambda_1$ and $\\lambda_2$, respectively, and analyze the increments $N_t - N_s$ of the sum $N_t = N^{(1)}_t + N^{(2)}_t$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----- student answer for Problem 4 goes here -----"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Problem 5: Now we consider the opposite situation: customers arrive at a rate of $\\lambda$, but each new customer independently decides to join one of two queues with respective probabilities $p$ and $1-p$. That is, we have a Poisson process $N = (N_t)_{t \\ge 0}$ with rate $\\lambda$, and then we form two counting processes $N^{(1)} = (N_t)_{t \\ge 0}$ and $N^{(2)} = (N^{(2)}_t)_{t \\ge 0}$ as follows:\n",
"\n",
"Let $T_1,T_2,\\ldots$ be the arrival times of $N$, and let $U_1,U_2,\\ldots$ be i.i.d. ${\\rm Bern}(p)$ random variables that are also independent of all the arrival times. Then\n",
"\n",
"$$\n",
"N^{(1)}_t = \\sum^{N_t}_{i=1} U_{N_{T_i}}\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
"N^{(2)}_t = \\sum^{N_t}_{i=1} (1-U_{N_{T_i}}).\n",
"$$\n",
"\n",
"Use your code from Problem 3 to model such a situation. Plot the resulting sample paths of $N$, $N^{(1)}$, and $N^{(2)}$ with the following parameters: $\\lambda = 2$, $p = 0.4$, $t_{\\max} = 10$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"###### student code for Problem 5 goes here ######"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the markdown cell below, comment on what you see. Do you think that the two counting processes $N^{(1)}$ and $N^{(2)}$ are also Poisson processes? If so, why do you think so, and what do you think their arrival rates are?\n",
"\n",
"Hint: Compute the distribution of $N^{(1)}_t$ using the law of total probability as follows:\n",
"\n",
"\\begin{align*}\n",
"{\\bf P}[N^{(1)}_t = m] &= \\sum^\\infty_{k=0} {\\bf P}[N^{(1)}_t = m\\, |\\, N_t = k] {\\bf P}[N_t = k] \\\\\n",
"&= \\sum^\\infty_{k=0} {\\bf P}[U_1 + \\ldots + U_k = m\\, |\\, N_t = k] {\\bf P}[N_t = k],\n",
"\\end{align*}\n",
"\n",
"then use the fact that $U_1 + \\ldots + U_n \\sim {\\rm Binom}(n,p)$ and that the $U_i$'s are independent of $N$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----- student answer for Problem 5 goes here -----"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [conda root]",
"language": "python",
"name": "conda-root-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.12"
}
},
"nbformat": 4,
"nbformat_minor": 0
}