{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Chapter 13: Poisson Processes\n", " \n", "This Jupyter notebook is the Python equivalent of the R code in section 13.5 R, pp. 534 - 536, [Introduction to Probability, 1st Edition](https://www.crcpress.com/Introduction-to-Probability/Blitzstein-Hwang/p/book/9781466575578), Blitzstein & Hwang.\n", "\n", "----" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1D Poisson process\n", "\n", "In Chapter 5, we discussed how to simulate a specified number of arrivals from a one-dimensional Poisson process by using the fact that the interarrival times are i.i.d. Exponentials. In this chapter, Story 13.2.3 tells us how to simulate a Poisson process within a specified interval $(0, \\left.L\\right]$. We first generate the number of arrivals $N(L)$, which is distributed $Pois(\\lambda L)$. Conditional on $N(L) = n$, the arrival times are distributed as the order statistics of $n$ i.i.d. $Unif(0, L)$ r.v.s. The following code simulates arrivals from a Poisson process with rate 10 in the interval $(0, \\left. 5\\right]$:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "np.random.seed(2971215073)\n", "\n", "from scipy.stats import poisson, uniform\n", "\n", "L = 5 \n", "lambd = 10 \n", "n = poisson.rvs(lambd*L, size=1)[0]\n", "t = np.sort(uniform.rvs(loc=0, scale=(L-0), size=n))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To visualize the Poisson process we have generated, we can plot the cumulative number of arrivals $N(t)$ as a function of $t$:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtQAAAEYCAYAAAB872GiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcZHV57/HPVxY1AsGRkQybg4oLGLdMQKMmKupFJUKiMYoaTVBMbkz0qlcwm6Im4tW4XTVXElR8yaIBDbhEwQW36MhAEGSLKPuMzMgwwhiVxef+cU5jTdtLdVd3naruz/v16ldXnTrLU1Wnqp566jm/k6pCkiRJ0vzcpesAJEmSpHFmQi1JkiQNwIRakiRJGoAJtSRJkjQAE2pJkiRpACbUkiRJ0gBMqKVFkuRDSd7U0baT5INJbkryrUXczsVJHr8A67kqyZPmuMzjklw+6LbnKsnWJPcd9na1tCT5f0n+dp7L/nuSFy50TIshSSW5/wy3vzTJOxdhu09qX6s/n+t7izQfJtRaNtqk7YYk9+iZ9uIk53QY1mJ5LPBkYK+qOnCxNlJVB1TVOYu1/lm2/dWqeuB8lk3yoiR3tB+4Nye5IMmhfW53p6r6/ny2q8WR5PVJPrII612V5Mwk69vEcPWk2++a5APtPvSDJK/sd91V9adV9cb5xFVVT62qE/uZN8k5SV48n+0stiQ7An8DvLVn2sOTnJfkv9v/D59h+XOS/LR9HW/t/YJdVZ+vqp2Aaxb1TkgtE2otN9sDL+86iLlKst0cF7kPcFVV/XgBtr19P9PG0DfaD9xdgROAjyVZ0XFMi2Zcn7OO4/458FngmdPc/npgP5rX2xOA1yQ5ZDihDcc83nvm4jDgsqq6vt3WjsAZwEeAewInAme006fzsvZL7k7z/YItLQQTai03bwVenWTXyTckWd1WobbvmXZndaetan49yTuSbEny/SS/1U6/NsnGKX6G3S3J2UluSfLlJPfpWfeD2ts2J7k8ybN7bvtQkn9K8pkkP6b5sJ4c7x5t9WxzkiuSvKSdfiTwL8Cj26rNsVMse78kX0xyY5IfJjmp9zFpq/lHJ7kQ+HGS7WeY9qQ2lp/0JqRJHtGue4fZtjcptgOTrGurfjckefs08z0+yXWTYn51kguT/CjJR5Pcbaple1XVz4EPAHcH7tuu6yXtY7q5fYz36NnOnT9hJ3lakkva5/f6JK9up++W5FPtfrI5yVeT3KW97cHtfrUlTcvMM3rW/aEk703y6Xada5Pcb5r7P7G/HpWmgrohyat6bn99ktOSfCTJzcCL0lRU39nOv769fNeeZQ5LU62/Ocn30iaHSX41yQntNq5P8qa0iVaS+7f79o/a5/aj7fSkea1sbG+7MMlDZns+eh7jP0/yXeC77bR3pXmd3Zymcvm4dvohwF8Bf9ju79+eLeZ+VdUNVfU+4NxpZvkj4I1VdVNVXQr8M/CiPu/jnS1hE/tykle1j9eGJH88w7KT35e+luRtaVq8rkzy1Pa2vwceB7ynfWze006fy3vPa9NU37frmef30rwPTLxev9HuzxuSvCczJ8C9ngp8uef642mKHu+sqp9V1buBAE/sc31SZ0yotdysA84BXj3P5Q8CLgTuBZwMnAr8JnB/4Pk0H1w79cz/POCNwG7ABcBJAGnaTs5u13Fv4LnA+5Ic0LPsEcDfAzsDX5sillOA64A9gGcB/5Dk4Ko6AfhT2gpsVb1uimUDvLld9sHA3jTVtl7PBZ4O7FpVt88wjapaD3yDbSt5RwCnVdVtfW5vwruAd1XVLsD9gI9NM99Ung0cAuwLPJQ+kps0X6BeDGwFvpvkiW2szwZWAVfTPM9TOQF4aVXtDDwE+GI7/VU0z81KYHeahK+S7AB8EjiL5nn/C+CkJL2VtecCx9JU6K6g2Qdm8gSaKulTgGOybb/oYcBpNFX4k4C/Bh4FPBx4GHAgzU/uJDkQ+DDwv9v5fxu4ql3PicDtNPv5I9ptTbQRvLG9P/cE9gL+bzv9Ke06HtCu7w+BG2e5L70Op3m97d9eP7eNewXN6+Zfk9ytqj4L/APw0XZ/f9hsMSd5bJsATvf32NmCS3JPmv352z2Tvw0cMPUSs/o14FeBPYEjgfe22+jHQcDlNO8z/wc4IUmq6q+Br/KLKu7L5vHe8zbgx2yb1B7RLg9wB/C/2m0/GjgY+J99xv3rbdwTDgAurKrqmXYhMz+mb26/yH09C3A8hzRfJtRajv4O+IskK+ex7JVV9cGqugP4KE1i+Ia2mnIWcCvNB/iET1fVV6rqZzTJzKOT7A0cStOS8cGqur2qzgdOp0mMJ5xRVV+vqp9X1U97g2jX8Vjg6Kr6aVVdQFOVfkE/d6Kqrqiqs9u4NwFvB35n0mzvrqprq+ons0ybcDLNhzNJAjynndbv9ibcBtw/yW5VtbWqvtnPfeqJb31VbaZJXKftvwQelWQL8IM27t+rqh/RfAn6QFWd3z5vr6V53lZPE+v+SXZpq5Tn90xfBdynqm5r+72LJpndCTiuqm6tqi8Cn2q3P+HjVfWt9gvLSbPcB4Bjq+rHVXUR8MFJ6/pGVf1buw/9pL1vb6iqje3zcCy/2GeObO/32e3811fVZUl2p6kkvqLdzkbgHTTP78R9vQ+wR7svfq1n+s7Ag4BU1aVVtWGW+9LrzVW1eWJfq6qPVNWN7evlH4G7AlP+xD9bzFX1taradYa/qb7ATjbxxflHPdN+1N7n+biN5rm5rao+Q/MFr98Whqur6p/b96UTafa93aeZdz7vPafwi9f2zsDT2mlU1XlV9c12XVcB72f61/ZkuwK39FzfiW0fT5j5MT2a5lelPYHjgU9mml90pMVmQq1lp6q+Q5PEHDOPxW/ouTzxQT95Wm+F+tqe7W4FNtNUte4DHNRbFaNJdn5tqmWnsAewuap6P4yupvlgmVWSeyc5tf0p/GaansXdJs021fZniuk0msRzD5rKZNFUx/rd3oQjaaqalyU5N30eLNj6Qc/l/2bb52Kyb7bJ025V9aiq+nw7fQ+axxK483m7kakf22fSJBdXt20Pj26nv5WmunxWmtagiX1tD+DaatpMJkx+3uZyH2Db5+TqdhtT3Tax/at7rvfOvzfwvSnWfx9gB2BDz776fprqJsBraH6B+FaaFpY/AWi/LLwHeC9wQ5Ljk+wyy32Z7n7RtkNc2raPbKGp5k63D80W80LY2v7vvU+7sG2COBc39v7qQ3/P/YQ795mq+u/24nTLzue952Tg99O0B/0+cH5VXQ2Q5AFp2pt+0L62/4Hpn5fJbmLbZHkr2z6eMMNjWlVrq+qW9ov6icDXaV6P0tCZUGu5eh3wErZNZCYO4PuVnmm9HzLzsffEhbYVZAWwnuYD68uTqmI7VdWf9SxbTG89sKKtFk3YB7i+z7je3K7/odW0VjyfJinqNdX2p42pqrbQ/PT/bJqfhE/p+em2n+1NrOe7VfVcmuTnLcBp6RmZZQjW0yQdwJ3tOfdiise2qs6tqsNoYv032vaU9kP+VVV1X+B3gVcmObhd995p+6lbc3neprJ3z+V92m3cGeKkebe5b5Pmv5amxWaya4GfAbv17Ku7VNUBAFX1g6p6SVXtAbyUpn3g/u1t766q36D5yf4BNO0k/boz9jT90kfT7Fv3rKpdaSqXmTxvPzGnGXJx6wx/j5s1uKqbgA00rTMTHgZcPIf7OAxTPTZzeu+pqktovnw9lW3bPQD+CbgM2K99bf8V07y2p3AhzX4x4WLgoe0vXBMeSv+Pac1h29KCMqHWslRVV9C0bPxlz7RNNInN85Ns11baBv358Gltv+aONL2ma6vqWpoK+QOSvCDNQXs7JPnNJA/uM/5rgf+g6R+8W5KH0lR2T+ozrp1pqkFbkuzJ3BKdmZxMc6DWM9n2Q7fv7SV5fpKVbRV3Szv5jgWKrx8nA3+cZviuu9JU3Na2P2f3xrljkucl+dVq+sRvnogzyaFpDtZLz/Q7gLU0X9xe0z7nj6dJuKfr0e7H3yb5lbYH9o9p9uvpnAL8TZKVSXajaX+aGG7uhPZ+H5zkLkn2TPKgtk3jLOAfk+zS3na/JL/T3tc/SLJXu46baJKaO9r9+aA0feM/Bn7a8/i8KMlVc7iPO9P0Q28Ctk/yd2xbybwBWD3xRWW2mNsWnJ1m+PvqxIrTHNg6ceDmXbPtga4fbh/PeyZ5EM2X9A/1LFvpvq/3BtqDbVvzfe85meb98reBf+2ZvjPNPr61fQz+bIplp/MZtm0POYdmH/nLNAfQvqyd/kXYdr9JsmuS/9G+/22f5HltbJ+bw/alBWNCreXsDcDkyudLaJK9G2mqav8x4DZOpqmGbwZ+g+anVdpWjafQ9HSup/nJ9i384oO7H88FVrfLfwJ4XVWd3eeyxwKPpKnyfRr4+By2O5MzaQ6Qu6Gqeg/Wmsv2DgEuTrKV5gDF59SkHvLFVFVfAP6Wpq90A82XqudMM/sLgKvan7r/lKbyDs1j8HmaLxHfAN5XVedU1a3AM2gqfT8E3gf8UVVdNkDIX6ZpL/kC8LZqevmn8yaaA3MvBC4Czm+nUVXfoknI30HzPH2ZX1Sz/wjYEbiEJmk+jaZPF5qDcte2z9eZwMur6kqahPef2/mvpnlNva1dZm+an+f79Tng34H/atf1U7ZtS5hI8G5MMtHHPlPMc/ETftHecVl7fcLraNpkrqZ5vN5azUGStF8yttI8zl16F/CsNCOAvHuA955TaEbh+GJV/bBn+qtpqta30DzfM32hm+yTwIPSjqLTvj4Op3nutgB/AhzeTodt95sdaPbdTTSvpb9o5x36yZ4kaA4U6ToGSdIcpTlI8kpgh0m9tyMvyVk0ifelXceyWJI8Hzigql7bdSyjLMlRwP5V9Yo+5u17v2lbrE6n+aLwtKr60sDBSjMwoZakMTTOCbUkLTW2fEiSJEkDsEItSZIkDcAKtSRJkjSA7bsOYK522223Wr16dddhSJIkaQk777zzflhVfZ1VeewS6tWrV7Nu3bquw5AkSdISluTq2edq2PIhSZIkDcCEWpIkSRqACbUkSZI0ABNqSZIkaQAm1JIkSdIATKglSZKkAZhQS5IkSQMwoZYkSdKydOwnL+bYT1488HrG7sQukiRJ0kK4ZP3NC7IeK9SSJEnSAEyoJUmSpAHY8iFJkqSxdvLaazjjguvnvNwlG25m/1W7DLx9K9SSJEkaa2dccD2XbJh7P/T+q3bhsIfvOfD2rVBLkiRp7O2/ahc++tJHd7JtK9SSJEnSAKxQS5Ikaez09k0vVC/0fFmhliRJ0tjp7ZteqF7o+bJCLUmSpLHUZd90LxNqSZIkDd18h7qb0HWbRy9bPiRJkjR08x3qbkLXbR69rFBLkiSpE6PSsjEoK9SSJEnSAKxQS5Ikad66Pu33KBhahTrJVUkuSnJBknXttBVJzk7y3fb/PYcVjyRJkgbX9Wm/R8GwK9RPqKof9lw/BvhCVR2X5Jj2+tFDjkmSJEkDWCq90PPVdQ/1YcCJ7eUTgcM7jEWSJEmas2FWqAs4K0kB76+q44Hdq2oDQFVtSHLvIcYjSZK0rA06FjQsrV7o+RpmQv2YqlrfJs1nJ7ms3wWTHAUcBbDPPvssVnySJEnLykT/8yAJ8VLqhZ6voSXUVbW+/b8xySeAA4Ebkqxqq9OrgI3TLHs8cDzAmjVralgxS5IkLXXLvf95IQwloU5yD+AuVXVLe/kpwBuAM4EXAse1/88YRjySJEnjZCFaM6Ziu8bCGFaFenfgE0kmtnlyVX02ybnAx5IcCVwD/MGQ4pEkSRobC9GaMRXbNRbGUBLqqvo+8LAppt8IHDyMGCRJksaZrRmjq+th8yRJkqSxZkItSZI0wk5eew1rr9zcdRiagQm1JEnSCJs4GNFe59FlQi1JkjTiDtp3BUcc5Lk4RpUJtSRJkjSAYZ4pUZIkacEt1hjNo8KxokefFWpJkjTWJsZoXqocK3r0WaGWJEljzzGa1SUr1JIkaWw5pJxGgQm1JEkaWw4pp1FgQi1JksaaQ8qpaybUkiRJ0gA8KFGSpGVmKQ0z55ByGgVWqCVJWmaW0jBzDimnUWCFWpKkZchh5qSFY4VakiRJGoAVakmSxsBC9j3bdywtLCvUkiSNgYXse7bvWFpYVqglSRoT9j1Lo8mEWpKkETbR6mGbhjS6bPmQJGmE9SbTtmlIo8kKtSRJI85WD2m0WaGWJEmSBmCFWpKkEdM7RJ6909Los0ItSdKI6R0iz95pafRZoZYkaQTZNy2NDxNqSZIWgGcylJYvWz4kSVoAnslQWr6GWqFOsh2wDri+qg5Nsi9wKrACOB94QVXdOsyYJElaKLZpSMvTsCvULwcu7bn+FuAdVbUfcBNw5JDjkSRJkgYytIQ6yV7A04F/aa8HeCJwWjvLicDhw4pHkqSFcvLaa1h75eauw5DUkWFWqN8JvAb4eXv9XsCWqrq9vX4dMGXDWJKjkqxLsm7Tpk2LH6kkSXMwcTCifc/S8jSUhDrJocDGqjqvd/IUs9ZUy1fV8VW1pqrWrFy5clFilCRpEAftu4IjDtqn6zAkdWBYByU+BnhGkqcBdwN2oalY75pk+7ZKvRewfkjxSJIkSQtiKAl1Vb0WeC1AkscDr66q5yX5V+BZNCN9vBA4YxjxSJI0m7mMK+240dLy1vU41EcDr0xyBU1P9QkdxyNJEjC3caUdN1pa3oZ+psSqOgc4p738feDAYccgSVI/HFdaUj+6rlBLkjRyHAZP0lyYUEuSNInD4EmaCxNqSZKm4DB4kvplQi1JkiQNYOgHJUqSNEqmGh7PYfAkzYUVaknSsjbV8HgOgydpLqxQS5KWPYfHkzQIK9SSJEnSAKxQS5KWpH5PHW6/tKRBWaGWJC1J/Z463H5pSYOyQi1JWrLsjZY0DCbUkqSxNl1rh60ckobFlg9J0librrXDVg5Jw2KFWpI09mztkNQlK9SSJEnSAKxQS5JGRr9D3fWyV1pS16xQS5JGRr9D3fWyV1pS16xQS5JGiv3QksaNFWpJkiRpAFaoJUmdm+idth9a0jiyQi1J6lxvMm0/tKRxY4VakjQS7J2WNK5MqCVJv2Q+w9cNwlYPSePMlg9J0i+Zz/B1g7DVQ9I4s0ItSZqSLRiS1B8r1JIkSdIArFBL0jLSb2+0Pc2S1D8r1JK0jPTbG21PsyT1bygV6iR3A74C3LXd5mlV9bok+wKnAiuA84EXVNWtw4hJkpYre6MlaWENq0L9M+CJVfUw4OHAIUkeBbwFeEdV7QfcBBw5pHgkSZKkBTHnCnWSewA/rao7+l2mqgrY2l7dof0r4InAEe30E4HXA/8015gkaSlZzDGg7Y2WpIU3a4U6yV2SHJHk00k2ApcBG5JcnOStSfbrZ0NJtktyAbAROBv4HrClqm5vZ7kOmLJhL8lRSdYlWbdp06Z+NidJY2sxx4C2N1qSFl4/FeovAZ8HXgt8p6p+DpBkBfAE4Lgkn6iqj8y0krai/fAkuwKfAB481WzTLHs8cDzAmjVrppxHkpYS+5wlaXz0k1A/qapuS3KfiWQaoKo2A6cDpyfZod8NVtWWJOcAjwJ2TbJ9W6XeC1g/t/AlaWk5ee01rL1yMwftu6LrUCRJfZq15aOqbmsvfmLybe2Bhb3zTCnJyrYyTZK7A08CLqWpfj+rne2FwBl9Ry5JS9BE77RtGZI0PmatUCd5NvBIYOckDwb+q+eAxOOBh/axnVXAiUm2o0niP1ZVn0pyCXBqkjcB/wmcMJ87IUlLyUH7ruCIg/bpOgxJUp/6afn4OnA34MXA24EHJtlC057xk342UlUXAo+YYvr3gQP7jlaSJEkaMbMm1FV1PfDhJN+rqq/DnQck7ksz4ockjZzFHHpuMTmsnSSNn36GzQvARDLdXt5cVedV1Y9755GkUbGYQ88tJoe1k6Tx09eweUlOB86oqmsmJibZEXgszcGEXwI+tCgRStI8OfScJGkY+kmoDwH+BDglyX1pThF+d5rq9lk0pw6/YPFClCRJkkZXPwn1vavqfcD72vGmdwN+UlVbFjc0SZq7id5pe5ElScPST0L92ST3phk3+iLgQuCiJBdV1fg1KEpa0nqTaXuRJUnD0M8oH/u3/dIHAL9OM+704cBDk/ysqvZd5BglaU7snZYkDdOso3wAVNWtVfWfNGdLXAv8gGYM6m8vYmySNCcTp+2WJGmY+jlT4gOBpwOHAiuBs4GTgKOq6tbFDU+S+udpuyVJXeinh/pSmtOCHwecWVU/W9yQJGn+PG23JGnY+kmo/4ymd/plwHuT3EhzcOJFwEVV9W+LGJ8kSZI00vo5KPH9vdeT7EVzYOKvA88ETKglzctCnx7cofIkSV3op0K9jaq6DrgO+MzChyNpOVno8aIdKk+S1IU5J9SStJAc4k6SNO76GjZPkiRJ0tSsUEsamsk90/Y8S5KWAivUkoZmomd6gj3PkqSlwAq1pKGyZ1qStNSYUEtacNMNh2eLhyRpKbLlQ9KCm9zaMcEWD0nSUmSFWtKisLVDkrRcWKGWJEmSBmCFWtK05ntqcHulJUnLiRVqSdOarhd6NvZKS5KWEyvUkmZkL7QkSTOzQi1JkiQNwAq1tIzMtSfaXmhJkmY3lAp1kr2TfCnJpUkuTvLydvqKJGcn+W77/57DiEdarubaE20vtCRJsxtWhfp24FVVdX6SnYHzkpwNvAj4QlUdl+QY4Bjg6CHFJC1L9kRLkrSwhpJQV9UGYEN7+ZYklwJ7AocBj29nOxE4BxNqaRvzHbpuKrZwSJK08IZ+UGKS1cAjgLXA7m2yPZF033uaZY5Ksi7Juk2bNg0rVGkkzHfouqnYwiFJ0sIb6kGJSXYCTgdeUVU3J+lruao6HjgeYM2aNbV4EUqjyTYNSZJG19Aq1El2oEmmT6qqj7eTb0iyqr19FbBxWPFIkiRJC2FYo3wEOAG4tKre3nPTmcAL28svBM4YRjzSuDh57TWsvXJz12FIkqQZDKvl4zHAC4CLklzQTvsr4DjgY0mOBK4B/mBI8UhjYeJgRPueJUkaXcMa5eNrwHQN0wcPIwZpXB207wqOOGifrsOQJEnT8NTjkiRJ0gA89bg0IqYab9pxoyVJGn1WqKURMdV4044bLUnS6LNCLY0Qx5uWJGn8mFBL87CQpwOfYHuHJEnjyZYPaR4W8nTgE2zvkCRpPFmhlubJ9gxJkgRWqCVJkqSBWKHWsrGQfc/2O0uSpAlWqLVsLGTfs/3OkiRpghVqLSv2PUuSpIVmhVqSJEkagBVqLSkz9Unb9yxJkhaDFWotKTP1Sdv3LEmSFoMVai059klLkqRhMqHW2LGtQ5IkjRJbPjR2bOuQJEmjxAq1xpJtHZIkaVRYoZYkSZIGYIVai24hT/kN9klLkqTRYoVai24hT/kN9klLkqTRYoVaQ2HPsyRJWqqsUEuSJEkDsEKtRdHbN23PsyRJWsqsUGtR9PZN2/MsSZKWMivUWjT2TUuSpOXAhFoDmW5IPNs8JEnScmHLhwYy3ZB4tnlIkqTlYigV6iQfAA4FNlbVQ9ppK4CPAquBq4BnV9VNw4hHC8vWDkmStJwNq0L9IeCQSdOOAb5QVfsBX2ivS5IkSWNlKAl1VX0F2Dxp8mHAie3lE4HDhxGLFs7Ja69h7ZWTn1ZJkqTlpcse6t2ragNA+//e082Y5Kgk65Ks27Rp09AC1MwmDka0V1qSJC1nY3FQYlUdX1VrqmrNypUruw5HPQ7adwVHHLRP12FIkiR1psuE+oYkqwDa/xs7jEWSJEmaly7HoT4TeCFwXPv/jA5jWdamG0t6No41LUmSNKQKdZJTgG8AD0xyXZIjaRLpJyf5LvDk9ro6MN1Y0rNxrGlJkqQhVair6rnT3HTwMLav2TmWtCRJ0vx46vExNN8WjenYuiFJkjR/YzHKh7Y13xaN6di6IUmSNH9WqMeULRqSJEmjwQq1JEmSNAAr1COmn/5oe54lSZJGhxXqEdNPf7Q9z5IkSaPDCvUIsj9akiRpfFihliRJkgZghboDM/VJ2x8tSZI0XqxQd2CmPmn7oyVJksaLFeqO2CctSZK0NJhQL6LpWjts65AkSVo6bPlYRNO1dtjWIUmStHRYoV5ktnZIkiQtbVaoJUmSpAFYoV4gU/VL2ystSZK09FmhXiBT9UvbKy1JkrT0WaFeQPZLS5IkLT9WqCVJkqQBLPsK9UynAZ8L+6UlSZKWp2VfoZ7pNOBzYb+0JEnS8rTsK9Rg77MkSZLmb2wTals1JEmSNArGtuXDVg1JkiSNgrGtUIOtGpIkSere2FaoJUmSpFFgQi1JkiQNoPOEOskhSS5PckWSY/pdbv89dmH/PTyYUJIkSd3qtIc6yXbAe4EnA9cB5yY5s6oumW3Z1/3uAYsdniRJkjSrrivUBwJXVNX3q+pW4FTgsI5jkiRJkvrWdUK9J3Btz/Xr2mmSJEnSWOg6oc4U0+qXZkqOSrIuybpNmzYNISxJkiSpP10n1NcBe/dc3wtYP3mmqjq+qtZU1ZqVK1cOLThJkiRpNl0n1OcC+yXZN8mOwHOAMzuOSZIkSepbp6N8VNXtSV4GfA7YDvhAVV3cZUySJEnSXHR+6vGq+gzwma7jkCRJkuaj65YPSZIkaayl6pcG1RhpSW4BLu86Do2c3YAfdh2ERo77habifqGpuF9osgdW1c79zNh5y8c8XF5Va7oOQqMlyTr3C03mfqGpuF9oKu4XmizJun7nteVDkiRJGoAJtSRJkjSAcUyoj+86AI0k9wtNxf1CU3G/0FTcLzRZ3/vE2B2UKEmSJI2ScaxQS5IkSSPDhFqSJEkawNgk1EkOSXJ5kiuSHNN1PBoNST6QZGOS73Qdi0ZDkr2TfCnJpUkuTvLyrmNS95LcLcm3kny73S+O7TomjY4k2yX5zySf6joWjYYkVyW5KMkF/QyfNxY91Em2A/4LeDJwHXAu8NyquqTTwNS5JL8NbAU+XFUP6ToedS/JKmBVVZ2fZGfgPOBw3y+WtyQB7lFVW5PsAHwNeHlVfbPj0DQCkrwSWAPsUlWHdh2PupfkKmBNVfV1sp9xqVAfCFxRVd+vqluBU4HDOo5JI6CqvgJs7jryg1YWAAACnUlEQVQOjY6q2lBV57eXbwEuBfbsNip1rRpb26s7tH+jX1HSokuyF/B04F+6jkXja1wS6j2Ba3uuX4cfkJJmkWQ18AhgbbeRaBS0P+tfAGwEzq4q9wsBvBN4DfDzrgPRSCngrCTnJTlqtpnHJaHOFNOsLEiaVpKdgNOBV1TVzV3Ho+5V1R1V9XBgL+DAJLaJLXNJDgU2VtV5XceikfOYqnok8FTgz9sW02mNS0J9HbB3z/W9gPUdxSJpxLU9sqcDJ1XVx7uOR6OlqrYA5wCHdByKuvcY4Bltv+ypwBOTfKTbkDQKqmp9+38j8Ama9uNpjUtCfS6wX5J9k+wIPAc4s+OYJI2g9uCzE4BLq+rtXcej0ZBkZZJd28t3B54EXNZtVOpaVb22qvaqqtU0ucUXq+r5HYeljiW5R3tQO0nuATwFmHE0sbFIqKvqduBlwOdoDjD6WFVd3G1UGgVJTgG+ATwwyXVJjuw6JnXuMcALaCpNF7R/T+s6KHVuFfClJBfSFGnOriqHSJM0ld2BryX5NvAt4NNV9dmZFhiLYfMkSZKkUTUWFWpJkiRpVJlQS5IkSQMwoZYkSZIGYEItSZIkDcCEWpIkSRqACbUkSZI0ABNqSZIkaQAm1JK0xCXZK8kfdh2HJC1VJtSStPQdDDyy6yAkaanyTImStIQleSxwBrAFuAX4vaq6stuoJGlpMaGWpCUuyWeBV1fVd7qORZKWIls+JGnpeyBweddBSNJSZUItSUtYknsBP6qq27qORZKWKhNqSVra9gXWdx2EJC1lJtSStLRdBuyW5DtJfqvrYCRpKfKgREmSJGkAVqglSZKkAZhQS5IkSQMwoZYkSZIGYEItSZIkDcCEWpIkSRqACbUkSZI0ABNqSZIkaQD/H/5fZfVvPr9+AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(12, 4))\n", "\n", "nt = np.arange(1, n+1)\n", "plt.plot(t, nt, linestyle='-', drawstyle='steps')\n", "\n", "plt.xlabel(r'$t$')\n", "plt.xlim([0,5])\n", "plt.ylabel(r'$N(t)$')\n", "plt.title(r'Number of arrivals in Poisson process, rate=10, in interval (0,5]')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This produces a staircase plot as in Figure 13.8." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Thinning\n", "\n", "The following code starts with an array of arrival times `t` and the corresponding number of arrivals `n`, generated according to the previous section. For each arrival, we flip a coin with probability `p` of Heads; these coin tosses are stored in the array `y`. Finally, the arrivals for which the coin landed Heads are labeled as type-1; the rest are labeled as type-2. The resulting arrays of arrival times $t_1$ and $t_2$ are realizations of independent Poisson processes, by Theorem 13.2.13. To perform this with $p = 0.3$, we can execute" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "np.random.seed(1)\n", "\n", "from scipy.stats import binom\n", "\n", "p = 0.3 \n", "y = binom.rvs(1, p, size=n)\n", "\n", "t1 = t[np.where(y==1)]\n", "t2 = t[np.where(y==0)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As before, we can plot the number of arrivals in each Poisson process, $N_{1}(t)$ and $N_{2}(t)$, as a function of $t$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtQAAAElCAYAAADEC3XRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmYHXWZ6PHvawwGQxQh4AQCJMPiIAwG7GFgcGEfdByVcXBglHEBo6OOKyp4VdTBKzwqMI7bjcp2VZArMKCgghpAdAQDBAQCwyq0iaRDCPsa3vtHVeOh08vpPt1VZ/l+nqefPqdOLW/Vqf6dt3/n91ZFZiJJkiRpYp5VdwCSJElSJzOhliRJklpgQi1JkiS1wIRakiRJaoEJtSRJktQCE2pJkiSpBSbUUo0i4pSIOKambUdEnBwR90bEFVO4nesjYs9JWM8dEbHvOJd5eUTc1Oq2p1JEfCMiPjnBZd8aEZdNdkytioiLI+LwuuMYTUS8KSIurDsOSd3BhFpqUCZtd0fEzIZph0fExTWGNVVeBuwHzM3MXadqI5m5Q2ZePFXrH2Pbv8zMF01k2aqS1cx8V2b++1RvpxM18x6UyfujEfFgRKyKiLMjYs5Y687M72bm/pMXraReZkItrevZwPvrDmK8ImLaOBfZCrgjMx+ahG0/u5lp0hR5b2ZuAGwHbAicUHM8U8a/K6k9mVBL6/oCcEREbDj0hYiYFxHZ+KHW+PV22aP2q4g4ISLWRMRtEfE35fS7ImJlRLxlyGpnR8RFEfFARFwSEVs1rPsvytdWR8RNEfHGhtdOiYivR8QFEfEQsNcw8W4WEeeVy98SEe8opx8GfAvYvezZ+8wwy24dEb+IiHvKnr/vNh6Tsjf/YxFxLfBQRDx7lGn7lrE8EhEbNaxj53Ld08fa3pDYdo2IJRFxf/mNwvEjzLdnRPQPifmIiLg2Iu6LiO9HxIxhltse+EbD8VkTEfPL388q5/lWRKxsWOY7EfGB0Y77CDE+PexnMN6I+HB5rqyIiLc1zLtxud77oxims/WQdY11vnyjhXPtqxFxfrns5RGxdcPr+0XEjeUx/QoQQ+J6e0Qsi2J40U+HbDcj4l0RcXP5+lejsM57MNIxHJSZq4GzgB3LdT8/Ik6LiIGI+H1EfKLh/Xu697vc3gnlMb+vPD8G1/HqiLih3O8/RMQRDbG/o3x/V5fvy2Zj7ddwcUfEpyPiB+X5+EBEXBURL2l4fbi/q+2jaHvWRDGs6rUN868fEV8q9/m+iLgsItYvX9stIn5dLndNNAzHKo/JbWUMt0fEm8rp25Tny31R/G1+v2GZ0c6bEY+d1HUy0x9//Cl/gDuAfYGzgWPKaYcDF5eP5wEJPLthmYuBw8vHbwWeBN4GTAOOAe4Evgo8B9gfeADYoJz/lPL5K8rX/wO4rHxtJnBXua5nA7sAq4AdGpa9D9iD4p/jGcPszyXA14AZwAJgANinIdbLRjkW21AMCXkOsAlwKXDikGO1FNgCWH+MafuWj38BvKNhHV8AvjGO7Q2u57+BQ8vHGwC7jbAPewL9Q9ZxBbAZsBGwDHjXCMuuc3zK9/Kl5eObgNuA7Rte23ms4z7Mdk7hT+fanhTnz2eB6cCrgYeBF5SvnwGcWZ4bOwJ/GOf50sq5thrYtXz9u8AZ5WuzgfuBfyxj/mC5D4N/E68HbgG2L5f9BPDrhv1P4EcUPctblsfqgGbO0WH+/mZTnGP/t3x+GnAuMIvib/d/gMOGrhv4W+DKMoYoY51TvrYCeHn5+AXALuXjvctjtEt5PP8TuLSZ/RpmHz4NPNFwDI8AbgemD/d3Vc5zC/BxYL0ylgeAF5Xzf7U8LptTtEN/U8a4OXAPxXn1LIq/t3so/t5mlu/j4DrmNLz/pwP/q1xmBvCyJs+bYY+dP/5044891NLwPgX8W0RsMoFlb8/MkzNzLfB9ig/Bz2bmY5l5IfA4RfI46PzMvDQzH6P40No9IrYAXkMxJOPkzHwyM6+i6H37x4Zlz83MX2XmU5n5aGMQ5TpeBnwsMx/NzKUUvdKHNrMTmXlLZl5Uxj0AHA+8cshsX87MuzLzkTGmDfoecEgZXwAHl9Oa3d6gJ4BtImJ2Zj6Ymb9pZp8a4lueRW/mDykS3mZdArwyIv6sfP6D8vl84HnANa0ed4p9+2xmPpGZFwAPAi+KYkjPG4BPZeZDmXkdcGrDcs2cL62ca2dn5hWZ+SRFQj143F4N3JCZP8jMJ4ATgT82LPdO4POZuaxc9n8DCxp7qYFjM3NNZt4JLGZ87wnAl8se7GsokrgPlcfrn4CjMvOBzLwD+BLDvw9PUCTdfwFEGeuKhtdeHBHPy8x7y2MD8CbgpMy8qjyeR1Ecz3kT3K8rG47h8RSJ626N+9jwd7UbxT+Sx2bm45n5C4rk/ZCyB/7twPsz8w+ZuTYzf13G+Gbggsy8oGwzLgKWULyHAE8BO0bE+pm5IjOvbzgGWwGblef04Lj2sc6bkY6d1HVMqKVhlMnKj4AjJ7D43Q2PHynXN3TaBg3P72rY7oMUPYGbUXyA/XX51eyaMmF4E/Bnwy07jM2A1Zn5QMO031P0Uo0pIjaNiDPKr2rvB75D0QPYaLjtjxbTDyiSjs0oekoT+OU4tjfoMIrxsjdGxG8j4jXN7FOpMdl7mGe+F2O5hKIX+RUUPegXUyT9rwR+mZlP0eJxB+4pE8+hMW5C0QvYeHx/3/B4XOfLBM61kY7bZkPWm0Ni3Ar4j4b1rqboBW48Hq28JwDvy8wNM3PzzHxT+Q/ZbIre28ZjNOz7UCakX6Ho2b07IhZFxPPKl99AkXD+vhz2sHvDfv++YR0PUvT2TnS/Go/hU0B/uY11Xi+n31XON3TfZlMk47cOs42tgIOGvM8vo+iNf4jiH5B3ASuiGN7zF+VyH6V4z64oh5e8vWF9o503Ix07qeuYUEsjOxp4B8/8gBws4Htuw7TGpGMithh8EBEbUAxFWE7xAXpJmSgM/myQmf/asGyOst7lwEYRMath2pYUwwSa8fly/Ttl5vMoereGjgEdbvsjxpSZa4ALgTcC/wycXiZgzW5vcD03Z+YhwKbAccAPouHKLJNkuP24BHg5RVJ9CXAZxZCbV5bPofXjPpIBiqEUWzRM27LhcTPnSyvn2khWDFlvDInxLuCdQ9a9fmb+uol1j3Z+j2UVf+pZHTTi+5CZX87MlwI7UPyz9pFy+m8z83UU59p/UQy5geK4NY4FnwlsPNL6m9B4DJ8FzC238XSIDY+XA1uU8w0a3LdVwKMMGV9fuotiOEzjezEzM48FyMyfZuZ+FMM9bgS+WU7/Y2a+IzM3o/jG4WsRsQ1jnDejHDup65hQSyPIzFsohmy8r2HaAMWH1psjYlrZUzPcB9d4vDoiXhYR6wH/DlyemXdR9JBvFxGHRlG0Nz0i/iqKYq1m4r8L+DXw+YiYERE7UfTsfrfJuGZRDDdYExGbUyYYk+B7wL9Q9F59byLbi4g3R8QmZQ/dYLHa2kmKb9DdwNzyfQGKRJ7iG4Y3U4yXvb+c7w2UCfUkHPdhZTGE6Gzg0xHx3Ih4MdBY4NrM+TIV59r5wA4R8Q9RFOu+j2f+k/kN4KiI2AGeLhQ8qMndXuc9aFZ5vM4EPhcRs8ohJh+i+ObjGcp9/euImE7xT/OjwNqIWC+K61U/vxyKcT9/Os++B7wtIhZExHMohrJcXg4tmYiXNhzDDwCPASMNZbq8jPOj5Xu1J/D3FOPanwJOAo6Pojh2WkTsXsb4HeDvI+Jvy+kzoiiEnRsRL4yI15b/GDxG8be4tjw+B0XE3HLb91Ik92sZ5bwZ49hJXceEWhrdZykKbxq9gyLZu4eiN6uZnrbRfI+iN3w18FKKr0wphwzsTzHOeDnF18fHURQXNesQimKs5cA5wNHluMlmfIaiyOg+iqTp7HFsdzTnAdsCd2fmNRPc3gHA9RHxIEVx3cE5ZAz5JPgFcD3wx4hY1TD9EophGXc2PA/g6oZ5Wjnuo3kvxbCBP1IUCp48+EKT58ukn2uZuQo4CDiW4m9iW+BXDa+fU67rjHIoz3XAq5rc35Heg2b9G0XieRvFtwnfo0g2h3oeRW/svRRDJ+4Bvli+dihwRxn7uyj+mSIzfw58kmLM8AqKf6wPnkCMg86lGHJxb7nNfygT0XVk5uPAaymO4yqKAth/ycwby1mOAH4H/JbivT4OeFb5z9PrKIoZByh6mD9CkQs8C/gwxfu/muJbl3eX6/sr4PLy7+08ivHZtzdx3gx77KRuFH/6tlWS1K0i4hSKK558ou5Y9EwR8Wlgm8w04ZQ6lD3UkiRJUgtMqCVJkqQWOORDkiRJaoE91JIkSVILTKglSZKkFphQS5IkSS0woZYkSZJaYEItSZIktcCEWpIkSWqBCbUkSZLUAhNqSZIkqQUm1JIkSVILTKglSZKkFphQS5IkSS0woZYkSZJaYEItSZIktcCEWpIkSWrBs+sOYCJmz56d8+bNqzsMSRq3K6+8clVmblJ3HFWyzZbUqZptszsyoZ43bx5LliypOwxJGreI+H3dMVTNNltSp2q2zXbIhyRJktQCE2pJkiSpBSbUkiRJUgs6cgy1pPbxxBNP0N/fz6OPPlp3KG1lxowZzJ07l+nTp9cdSlvyvBme543UmUyoJbWkv7+fWbNmMW/ePCKi7nDaQmZyzz330N/fz/z58+sOpy153qzL80bqXA75kNSSRx99lI033tikqEFEsPHGG9v7OgrPm3V53kidy4RaUstMitblMRmbx2hdHhOpM5lQS5IkScCZV/dz5tX9417OhFpSzzjxxBN5+OGHW17Pueeey0477cSCBQvo6+vjsssum4To1K4m67z5whe+wIIFC1iwYAE77rgj06ZNY/Xq1ZMQoaTJcteaR7hrzSPjXs6EWlJXyUyeeuqpYV+brMRon3324ZprrmHp0qWcdNJJHH744S2vU/Wq4rz5yEc+wtKlS1m6dCmf//zneeUrX8lGG23U8nol1a/ShDoipkXE1RHxo/L5/Ii4PCJujojvR8R6VcYjqTvccccdbL/99rz73e9ml1124bDDDqOvr48ddtiBo48+GoAvf/nLLF++nL322ou99toLgAsvvJDdd9+dXXbZhYMOOogHH3ywqe1tsMEGT491feihhxz32qGqPm8anX766RxyyCGTuj+S6lP1ZfPeDywDnlc+Pw44ITPPiIhvAIcBX684JkmT5OFfnM7alXdN6jqnbboFz9177MTjpptu4uSTT+ZrX/saq1evZqONNmLt2rXss88+XHvttbzvfe/j+OOPZ/HixcyePZtVq1ZxzDHH8LOf/YyZM2dy3HHHcfzxx/OpT32KD37wgyxevHidbRx88MEceeSRAJxzzjkcddRRrFy5kvPPP39S97nXnHl1/4S+Yh3NFhuuzxt3njvmfFWfNwAPP/wwP/nJT/jKV74yqfssqT6VJdQRMRf4O+BzwIei6NLZG/jncpZTgU9jQi1pArbaait22203AM4880wWLVrEk08+yYoVK7jhhhvYaaednjH/b37zG2644Qb22GMPAB5//HF23313AE444YQxt3fggQdy4IEHcumll/LJT36Sn/3sZ5O8R6pC1ecNwA9/+EP22GMPh3tINfrlrau44s5715nev+YR5m64/rjXV2UP9YnAR4FZ5fONgTWZ+WT5vB/YfKSFI2IhsBBgyy23nMIwJU1UMz3JU2XmzJkA3H777Xzxi1/kt7/9LS94wQt461vfOux1fTOT/fbbj9NPP32d15rtaQR4xStewa233sqqVauYPXv2JO1N5xtPm91MT/JUqeO8OeOMMxzuIdXsijvvHTZ5nrvh+uy65QvGvb5KEuqIeA2wMjOvjIg9BycPM2uOtI7MXAQsAujr6xtxPkm97f7772fmzJk8//nP5+677+bHP/4xe+65JwCzZs3igQceYPbs2ey222685z3v4ZZbbmGbbbbh4Ycfpr+/n+22227MnsZbbrmFrbfemojgqquu4vHHH2fjjTeuYO86R6e12VWcNwD33Xcfl1xyCd/5znemeI8kjWXuhuvz4b22nZR1VdVDvQfw2oh4NTCDYgz1icCGEfHsspd6LrC8ongkdamXvOQl7Lzzzuywww78+Z//+dNfzQMsXLiQV73qVcyZM4fFixdzyimncMghh/DYY48BcMwxx7DddtuNuY2zzjqL0047jenTp7P++uvz/e9/38LEDlfFeQPF2Pv999//6Z5xSd0hMqvtOCh7qI/IzNdExP8DzmooSrw2M7821jr6+vpyyZIlUx2qpCYsW7aM7bffvu4w2tJwxyYirszMvppCqsVwbbbnzcg8NtLU+9LimwHG7KFuts2u+zrUH6MoULyFYkz1t2uOR5IkSV3sl7eu4uaBhyZ1nVVfNo/MvBi4uHx8G7Br1TFIkiSpNw1e3WMixYcjqbuHWlIXqHroWCfwmIzNY7Quj4lUjW03mcnLt568KzOZUEtqyYwZM7jnnntMBBpkJvfccw8zZsyoO5S25XmzLs8bqXNVPuRDUneZO3cu/f39DAwM1B1KW5kxYwZz59Z3feV253kzPM8bqTOZUEtqyfTp05k/f37dYajDeN5Imioj3QVx0ETvhjgah3xIkiSpawzeBXEkE70b4mjsoZYkSVJXmcy7IDbDHmpJkiSpBSbUkiRJUgsc8iFJkqSOM1Lx4VQUHY7FHmpJkiR1nJGKD6ei6HAs9lBLkiSpI1VdfDgSe6glSZKkFphQS5IkSS1wyIckSZLaTh13PJwoe6glSZLUduq44+FE2UMtSZKkttQuRYdjqayHOiJmRMQVEXFNRFwfEZ8pp58SEbdHxNLyZ0FVMUmSJEmtqrKH+jFg78x8MCKmA5dFxI/L1z6SmT+oMBZJkiRpUlSWUGdmAg+WT6eXP1nV9iVJklS/sYoNB7VT0eFYKi1KjIhpEbEUWAlclJmXly99LiKujYgTIuI5Iyy7MCKWRMSSgYGBymKWJI2fbbakkYxVbDionYoOx1JpUWJmrgUWRMSGwDkRsSNwFPBHYD1gEfAx4LPDLLuofJ2+vj57tiWpjdlmSxpNpxQbNquWy+Zl5hrgYuCAzFyRhceAk4Fd64hJkiRJmogqr/KxSdkzTUSsD+wL3BgRc8ppAbweuK6qmCRJkqRWVTnkYw5wakRMo0jkz8zMH0XELyJiEyCApcC7KoxJkiRJNF8s2KpOKjZsVpVX+bgW2HmY6XtXFYMkSZKGN1gsONXJbicVGzbLOyVKkiQJ6L5iwarUUpQoSZIkdQsTakmSJKkFJtSSJElSC0yoJUmSpBaYUEuSJEktMKGWJEmSWmBCLUmSJLXA61BLkiS1garuVDiSbryDYVXsoZYkSWoDg3cqrEs33sGwKvZQS5IktQnvVNiZ7KGWJEmSWmBCLUmSJLXAIR+SJElTYLxFhhYFdi57qCVJkqbAeIsMLQrsXJX1UEfEDOBS4Dnldn+QmUdHxHzgDGAj4Crg0Mx8vKq4JEmSpopFhr2hyh7qx4C9M/MlwALggIjYDTgOOCEztwXuBQ6rMCZJkiSpJZUl1Fl4sHw6vfxJYG/gB+X0U4HXVxWTJEmS1KpKixIjYhpwJbAN8FXgVmBNZj5ZztIPbF5lTJIkSa0argDRIsPeUWlRYmauzcwFwFxgV2D74WYbbtmIWBgRSyJiycDAwFSGKUlqkW22es1wBYgWGfaOWi6bl5lrIuJiYDdgw4h4dtlLPRdYPsIyi4BFAH19fcMm3ZKk9mCbrV5kAWLvqqyHOiI2iYgNy8frA/sCy4DFwD+Ws70FOLeqmCRJkqRWVdlDPQc4tRxH/SzgzMz8UUTcAJwREccAVwPfrjAmSZIkqSWVJdSZeS2w8zDTb6MYTy1JktRRBosRLUDsbd4pUZIkaYIak2kLEHtXLUWJkiRJ3cJiRNlDLUmSJLXAHmpJkqQmDb2Bi2OnBfZQS5IkNW3oDVwcOy2wh1qSJGlcHDOtoeyhliRJklpgQi1JkiS1wCEfkiRJI7AIUc2wh1qSJGkEFiGqGfZQS5IkjcIiRI3FHmpJkiSpBSbUkiRJUgsc8iFJknre0OLDQRYhqhn2UEuSpJ43tPhwkEWIaoY91JIkSVh8qImrrIc6IraIiMURsSwiro+I95fTPx0Rf4iIpeXPq6uKSZIkSWpVlT3UTwIfzsyrImIWcGVEXFS+dkJmfrHCWCRJkqRJUVlCnZkrgBXl4wciYhmweVXblyRJvWGkAsPRWHyoVtRSlBgR84CdgcvLSe+NiGsj4qSIGHbkf0QsjIglEbFkYGCgokglSRNhm606jVRgOBqLD9WKyosSI2ID4CzgA5l5f0R8Hfh3IMvfXwLePnS5zFwELALo6+vL6iKWJI2XbbbqZoGhqlRpD3VETKdIpr+bmWcDZObdmbk2M58CvgnsWmVMkiRJUiuqvMpHAN8GlmXm8Q3T5zTMdiBwXVUxSZIkSa2qcsjHHsChwO8iYmk57ePAIRGxgGLIxx3AOyuMSZIktaGJFBYOssBQVavyKh+XATHMSxdUFYMkSeoMg4WFE0mMLTBU1bxToiRJaksWFqpT1HLZPEmSJKlbmFBLkiRJLRh3Qh0RMyNi2lQEI0lqb34GSNK6xkyoI+JZEfHPEXF+RKwEbgRWRMT1EfGFiHBwkyR1KT8DJGlszfRQLwa2Bo4C/iwzt8jMTYGXA78Bjo2IN09hjJKk+vgZIEljaOYqH/tm5hMRsVV5N0MAMnM1xV0PzyrvgChJ6j5+BkjSGMbsoc7MJ8qH5wx9LSJ2GzKPJKmL+BkgSWMbs4c6It4I7ALMiojtgf/JzLXly4uAnaYwPklSjfwMUDNauavhSLzboTpJM0M+fgXMAA4HjgdeFBFrgOXAI1MYmySpfn4GaEyt3NVwJN7tUJ1kzIQ6M/8AnBYRt2bmrwAiYiNgPkW1tySpS/kZoGZ5V0P1smaGfEQWfjU4rSxGWT10nimKUZJUEz8DJGlsTV02LyL+LSK2bJwYEetFxN4RcSrwlqkJT5JUMz8DJGkMzYyhPgB4O3B6RMwH1lCMp5sGXAickJlLpy5ESVKN/AzoUeMpNLSAUL2umTHUjwJfA75WXmt0NvBIZq6Z6uAkSfXyM6B3jafQ0AJC9bpmeqifVl5rdMUUxSJJamN+BvQeCw2l5jQzhnpUEfGxJufbIiIWR8SyiLg+It5fTt8oIi6KiJvL3/6LK0ltKiL2i4hvRsSC8vnCumOSpLqNq4caICLObHwKLACOa2LRJ4EPZ+ZVETELuDIiLgLeCvw8M4+NiCOBI4GmknRJUuXeDbwN+ER5+bwFNccjSbUbd0IN3J+Zhw8+iYivN7NQZq6g/KowMx+IiGXA5sDrgD3L2U4FLsaEWpLa1UA5fvqIiDgW+Ku6A9LkGixGtNBQat6YQz6GXioJ+NyQ5/9rvBuNiHnAzsDlwAvLZHsw6d50hGUWRsSSiFgyMDAw3k1KkiZgmM+A8wcfZOaRwGkjLGeb3aEak2kLDaXmNNND/ZOI2BRYBvwOuDYifgf8LjPvLy/w37SI2AA4C/hAZt4fEU0tl5mLgEUAfX193kBAkqox3GfAHvzpM+A/h1vINruzWYwojU8zl817cUSsB+wA/CWwE/B6YKeIeCwz5ze7sfKSS2cB383Ms8vJd0fEnMxcERFzgJXj3gtJ0pSYzM8ASepWTY2hzszHgasj4hbgEYrrkG5L0VvRlCi6or8NLMvM4xteOo/iLlvHlr/PbXadkqSpNxmfAZLUzcZMqCPiRcDfAa8BNgEuAr4LLCwb2WbtARwK/C4iBu+q9XGKRPrMiDgMuBM4aBzrlCRNoUn8DFDFxnOnw0YWI0rj10wP9TLgaorE97zMfGwiG8rMyyguszecfSayTknSlJuUzwBVb6JX6rAYURq/ZhLqf6UYN/de4KsRcQ/F13yDhYn/NYXxSZLq5WdAB7O4UKpGM0WJ/6fxeUTMpShK+UvgDYCNqSR1KT8DJGls476xS2b2A/3ABZMfjiSpnfkZIEnrmsidEiVJUo2aKTi0uFCqzph3SpQkSe1lsOBwNBYXStWxh1qSpA5kwaHUPuyhliRJklpgQi1JkiS1wIRakiRJaoEJtSRJktQCE2pJkiSpBSbUkiRJUgtMqCVJkqQWeB1qSZI6QOPdEb0LotRe7KGWJKkDNN4d0bsgSu2lsh7qiDgJeA2wMjN3LKd9GngHMFDO9vHMvKCqmCRJ6iTeHVFqT1X2UJ8CHDDM9BMyc0H5YzItSZKkjlJZQp2ZlwKrq9qeJEmSVIV2KEp8b0T8C7AE+HBm3lt3QJIkVamx4HAkFiJK7avuosSvA1sDC4AVwJdGmjEiFkbEkohYMjAwMNJskqQ2YJs9Po0FhyOxEFFqX7X2UGfm3YOPI+KbwI9GmXcRsAigr68vpz46SdJE2WaPnwWHUueqtYc6IuY0PD0QuK6uWCRJkqSJqPKyeacDewKzI6IfOBrYMyIWAAncAbyzqngkSZKkyVBZQp2Zhwwz+dtVbV/Suh675hIeX3Z53WFMimmbbsFz9x6umZHqZcGh1P3qLkqUVKPHl13O2pV31h2G1NUsOJS6XztcNk9SjaZtuiWzDv5o3WFIXc2CQ6m72UMtSZIktcCEWpIkSWqBQz6kLtVMweHalXcybdMtK4pI6m4jFR9acCh1P3uopS7VTMHhtE23ZL3t/7qiiKTuNlLxoQWHUvezh1rqYhYcStWy+FDqTfZQS5IkSS0woZYkSZJa4JAPqcsMFiNacChNHu92KGk09lBLXaYxmbbgUJoc3u1Q0mjsoZa6kMWI0uSz4FDSSOyhliRJklpgD7XUYca6YYtjp6XJ9ctbV3HzwENsu8nMukOR1KbsoZY6zFg3bHHstDS5BosRHR8taST2UEsdyDHSUrW23WQmL996dt1hSGpTlfVQR8RJEbEyIq5rmLZRRFwUETeXv/33X5IkSR2lyiEfpwAHDJl2JPDzzNwW+Hn5XJIkSeoYlQ35yMxLI2LekMmvA/YsH58KXAx8rKqYpHY2UvGhRYfS1Gu8kYs3bJE0lrqLEl+YmSsAyt+bjjRjRCyMiCURsWRgYKCyAKW6jFR8aNGhOkGnt9no8sHpAAAI+0lEQVSNN3Lxhi2SxtIxRYmZuQhYBNDX15c1hyNVwuJDdapuaLO9kYukZtXdQ313RMwBKH+vrDkeSZIkaVzqTqjPA95SPn4LcG6NsUiSJEnjVtmQj4g4naIAcXZE9ANHA8cCZ0bEYcCdwEFVxSNVYay7Go7G4kOpHt4ZUdJ4VXmVj0NGeGmfqmKQqjZYWDiRxNjiQ6ke3hlR0nh1TFGi1KksLJQ6j3dGlDQedY+hliRJkjqaCbUkSZLUAod8SEO0Ukg4lIWFUvtqvBtiI++MKGm87KGWhhjpDoUTYWGh1L4a74bYyDsjShove6ilYVhIKPUG74YoaTLYQy1JkiS1wIRakiRJaoFDPtRWJrMgcKIsJJQ630gFh40sPpQ0WeyhVluZzILAibKQUOp8IxUcNrL4UNJksYdabceCQEmTwYJDSVWxh1qSJElqgQm1JEmS1AKHfEiSOt7QIkQLDiVVyR5qSVLHG1qEaMGhpCrZQy1J6goWIUqqS1sk1BFxB/AAsBZ4MjP76o1IkiRJak5bJNSlvTJzVd1BSJIkSePRTgm1usxE7nroXQolNauxENEiREl1apeixAQujIgrI2LhcDNExMKIWBIRSwYGBioOTxMxkbseepdCqTtU0WY3FiJahCipTu3SQ71HZi6PiE2BiyLixsy8tHGGzFwELALo6+vLOoLU+HnXQ6k3VdVmW4goqR20RQ91Zi4vf68EzgF2rTciSZIkqTm1J9QRMTMiZg0+BvYHrqs3KkmSJKk57TDk44XAOREBRTzfy8yf1BuSmjVa4aEFhpImm4WIktpR7Ql1Zt4GvKTuODQxg4WHwyXOFhhKmmyDhYhzN1zfQkRJbaP2hFqdz8JDSVWyEFFSu6l9DLUkSZLUyUyoJUmSpBY45KNDTOSug1Ww8FDSVLMQUVK7s4e6Q0zkroNVsPBQ0lTzjoiS2p091B3E4j9JvcpCREntzB5qSZIkqQUm1JIkSVILunLIR7sW8LXC4j9JnerBx57kS4tvnvDyFiJKandd2UPdrgV8rbD4T1KneujxtU8XFU6EhYiS2l1X9lCDBXyS1E4sKpTUzbqyh1qSJEmqigm1JEmS1AITakmSJKkFXTmGetqmW9QdgiSptN60YAuv0iGpi7VFQh0RBwD/AUwDvpWZx7ayvufufcikxCVJat0Lnrseb9x5bt1hSNKUqX3IR0RMA74KvAp4MXBIRLy43qgkSZKk5tSeUAO7Ardk5m2Z+ThwBvC6mmOSJEmSmtIOCfXmwF0Nz/vLaZIkSVLba4eEOoaZluvMFLEwIpZExJKBgYEKwpIkTZRttqRe0g4JdT/QeFmOucDyoTNl5qLM7MvMvk022aSy4CRJ42ebLamXtENC/Vtg24iYHxHrAQcD59UckyRJktSU2i+bl5lPRsR7gZ9SXDbvpMy8vuawJEmSpKbUnlADZOYFwAV1xyFJkiSNV2SuU//X9iLiAeCmuuOo0GxgVd1BVMx97g29uM8vysxZdQdRpR5ss6E3z+1e2+de21/ozX1uqs1uix7qCbgpM/vqDqIqEbGkl/YX3Ode0av7XHcMNeipNht699zupX3utf2F3t3nZuZrh6JESZIkqWOZUEuSJEkt6NSEelHdAVSs1/YX3Ode4T73Bve5N/TaPvfa/oL7PKKOLEqUJEmS2kWn9lBLkiRJbaGjEuqIOCAiboqIWyLiyLrjmWoRcVJErIyI6+qOpSoRsUVELI6IZRFxfUS8v+6YplpEzIiIKyLimnKfP1N3TFWIiGkRcXVE/KjuWKoSEXdExO8iYmkvXO2j19ps6L122za7d9ps6L12ezxtdscM+YiIacD/APsB/RS3LD8kM2+oNbApFBGvAB4ETsvMHeuOpwoRMQeYk5lXRcQs4Erg9V3+PgcwMzMfjIjpwGXA+zPzNzWHNqUi4kNAH/C8zHxN3fFUISLuAPoys+uv49qLbTb0Xrttm907bTb0Xrs9nja7k3qodwVuyczbMvNx4AzgdTXHNKUy81Jgdd1xVCkzV2TmVeXjB4BlwOb1RjW1svBg+XR6+dMZ/+lOUETMBf4O+FbdsWjK9FybDb3Xbttm90abDbbbY+mkhHpz4K6G5/10+R9tr4uIecDOwOX1RjL1yq/RlgIrgYsys9v3+UTgo8BTdQdSsQQujIgrI2Jh3cFMMdvsHmOb3fV6sd1uus3upIQ6hpnW9f8R9qqI2AA4C/hAZt5fdzxTLTPXZuYCYC6wa0R07VfFEfEaYGVmXll3LDXYIzN3AV4FvKccHtCtbLN7iG1297bZ0NPtdtNtdicl1P3AFg3P5wLLa4pFU6gck3YW8N3MPLvueKqUmWuAi4EDag5lKu0BvLYcm3YGsHdEfKfekKqRmcvL3yuBcyiGRXQr2+weYZvd9W029Gi7PZ42u5MS6t8C20bE/IhYDzgYOK/mmDTJymKPbwPLMvP4uuOpQkRsEhEblo/XB/YFbqw3qqmTmUdl5tzMnEfxd/yLzHxzzWFNuYiYWRZtEREzgf2Bbr4ShG12D7DN7v42G3qz3R5vm90xCXVmPgm8F/gpRdHDmZl5fb1RTa2IOB34b+BFEdEfEYfVHVMF9gAOpfjvd2n58+q6g5pic4DFEXEtRRJyUWb2xCWJeswLgcsi4hrgCuD8zPxJzTFNmV5ss6En223bbNvsbjWuNrtjLpsnSZIktaOO6aGWJEmS2pEJtSRJktQCE2pJkiSpBSbUkiRJUgtMqCVJkqQWmFBLkiRJLTChliRJklpgQi01iIi5EfFPdcchSRqbbbbahQm19Ez7ALvUHYQkqSm22WoL3ilRKkXEy4BzgTXAA8CBmXl7vVFJkoZjm612YkItNYiInwBHZOZ1dcciSRqdbbbahUM+pGd6EXBT3UFIkppim622YEItlSJiY+C+zHyi7lgkSaOzzVY7MaGW/mQ+sLzuICRJTbHNVtswoZb+5EZgdkRcFxF/U3cwkqRR2WarbViUKEmSJLXAHmpJkiSpBSbUkiRJUgtMqCVJkqQWmFBLkiRJLTChliRJklpgQi1JkiS1wIRakiRJaoEJtSRJktSC/w8H91UgDjTXOwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "_, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12, 4))\n", "\n", "# graph for t1\n", "nt1 = np.arange(1, t1.size+1)\n", "ax1.plot(t1, nt1, linestyle='-', drawstyle='steps', color='#ef8a62', label='rate=3')\n", "ax1.set_xlabel(r'$t$')\n", "ax1.set_xlim([0,5])\n", "ax1.set_ylabel(r'$N_{1}(t)$')\n", "ax1.legend()\n", "\n", "# graph for t2\n", "nt2 = np.arange(1, t2.size+1)\n", "ax2.plot(t2, nt2, linestyle='-', drawstyle='steps', color='#67a9cf', label='rate=7')\n", "ax2.set_xlabel(r'$t$')\n", "ax2.set_xlim([0,5])\n", "ax2.set_ylabel(r'$N_{2}(t)$')\n", "ax2.legend()\n", "\n", "plt.suptitle('Number of arrivals in two independent Poisson processes')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2D Poisson process \n", "\n", "Simulating a 2D Poisson process is nearly as straightforward as simulating a 1D Poisson process. In the square $(0, \\left. L\\right] \\times (0, \\left. L \\right]$, the number of arrivals is distributed $Pois(\\lambda L^2)$. Conditional on the number of arrivals, the locations of the arrivals are Uniform over the square, so by Example 7.1.22, their horizontal and vertical coordinates can be generated independently:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "np.random.seed(1)\n", "\n", "L = 5 \n", "lambd = 10 \n", "n = poisson.rvs(lambd*(L**2), size=1)[0]\n", "x = uniform.rvs(loc=0, scale=(L-0), size=n)\n", "y = uniform.rvs(loc=0, scale=(L-0), size=n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then easily plot the locations of the arrivals as in Figure 13.7:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAGHCAYAAAC3RzktAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnX+4HVV5778vCa0KCJqTIiFCmsMvNQ2cJsUYi0mjgihYwac/pNDSSwvmFjFFW8CHlpa20PZpqJdrbyRU4RrE2Apya4AktpgEEgXyg5yChHhOiBiCcE7VFFAUkvf+MbMP++wze++ZPTPr5/fzPPOcs3/NvLNmzXe9613vWiOqCkIIIfFwkG0DCCGEmIXCTwghkUHhJ4SQyKDwE0JIZFD4CSEkMij8hBASGRR+QgiJDAo/IYYRkXeIyDdFZL2IfElEDrZtE4kLCj8h5vkugEWqugDALgC/btkeEhkU/jaIyGMisrDkPnaLyHtsHNvEPtP99nSOPR7rVhH56xr2KyLyXRHpz/Hd60VkSZnjqepeVf1J+vIVAAea9v+QiLytzP59oIpy7LL/+0TkJRF5oKL91WpveoxKbe5EsMIvIutE5Ici8vO9/F5V36aq6yo2y9qxbZ5Pg7obCRF5g4ioiLwgIj8Wkb15blZNOFZVh7vsfyqA3wVwU9N7j4vIniyxFpG+1J6j2+zvFwGcCWBV09v/AODabja37KetDb0iIpeKyGYR+amI3Jrx+RtF5Ksi8mLaaJ5XYN9Z5Vh4f53KV1UXAfhoXptM2GvS5m4EKfwiMgPAaQAUwAc7fG9ynveIN5wCYERVD1XV1wFYDOAfRWR6Rfu/EMA9Td46AMwCsBPAhzO+PwBgVFWfbv1ARF4P4P8CuEBVf9b00b8B+DUROaqAXZ1s6JW9AP4awOfbfP5PAH4G4EgAvwNgWYGG50JMLMde9te2fDuR9vAGMt4/WUQm1WhvzzZXTZDCj6R1/haAWwH8XvMHqdd5hYgMAnhRRCZ3eO89InKliHylZR//S0RuTP+/UkSGReR5Efm2iJzTzqj0GE+n331CRN7d5ntjnnH6/ydFZFBE9onIl0XkNUWPkbHPP0n3+aKIfE5EjhSRe9Pf/buIvKFpnyoixzW9bhtyaVceIrICwDEAvpZ65H+avj9NRO4QkREReVJELmva14CIbE339WUAbc875RQADze9fjD9+3PpzX5F6p39SET+RUQOT4/zByJyV9NxZ4rIKhEZTcv86+lHZwJY33xAVd0P4AEAJ7ex55GMMpoM4EsA/kJVn2jZ30sAtgA4vcu55rWhJ1T1TlW9C8B/tX4mIocgaWT+TFVfUNUHkDRYF+Tc/bhyLLG/zPLNwQwAa0XkfU02nAbgPwC8pUZ7y9hcKSEL/xfT7QwRObLl848A+ACAI1T1lQ7vAckN+v7UQ0PqEfwmgNvTz4eR9C4OB/CXAG7L8tZE5EQAlwL4FVU9DMAZAHbnPJ/fBPA+AL8IYDYSD2QCBY/xYQDvBXACgLMB3AvgUwD6kNSLy9r8rhuZ5aGqFwB4CsDZqUf+9yJyEICvAdgO4GgA7wawRETOEJGfA3AXgBUA3gjgX9Hdox0A8BAAiMgRAP4GiYg+CeCvkNzA8wC8CcDPA/jz9HezUxsafAFJeRyZbn+Rvv9LAMYJtYi8FsBvp/topd1N/hEAbwfw55KEJH+r5fPHUUDEu9jQ+M6qtMHL2la1+10bTgCwX1V3Nr23HUBej7+1HHvdX08iqqpPIqlLXxSRXxORUwHcCeB8VX20Rnt7trlqghN+EflVAMcC+BdV3YJEiFrjbzeq6vdaum5Z70FVvwtgK4APpW8tAvBjVf1W+vm/poN1B1T1ywC+A+DUDNP2IxGbt4rIwaq6u1tMucW2var6AyRCeUqb7xU5xv9W1WfTLuf9AB5U1W2q+lMAX0UiooUpUB4A8CsApqrqtar6M1XdBeBmJCI2D8DBAD6tqi+r6lcw3pvP4hQAfyIiP0Ai+IqkUfsFAB8DcJ6qPpN61V8BMDf93ckYL/z9ACYBmKSqL6nqxvT9IwA833LMvwHwNIB+ETk0w54JN7mqrlDVPlVdmG5fbvnK8+mx8tLJhsYxz1LVI9psZxU4FgAcCmBfy3v7AByW8/et5djr/noWUVXdgEQXvoJkjOViVV3d5utV2QtQ+Gvj9wCsVdXR9PXtaAn3APhexu+y3mtwOxIvDUgqS8Pbh4j8rog80vCekMRb+1p3oKpDAJYg8R6fE5GVIjItx/kAwPeb/v8xkoo3gYLHeLbp/59kvM48RjfylkfKsQCmNXufSHodRwKYBuBpHf/AiO92OO7PI+mmn6Sqb1TVflX9A1V9BkkP5D9VdW/TT/oAPJP+3+rx/w6SFMu9aRjsjen7P0TTzS0i70DSG/swkht/VtNnrwFwIoBt7WzuwGEAfpTni51sqJEXALy+5b3XY2Kj2I5x5djL/kqWb4OnkGRVCTr3vkvbC1RmcyUEJfxpl/c3ASwQke+LyPcB/DGAk0Wkueuc9fSZTk+k+VcACyUZJDwHqfCLyLFIPNRLAUxR1SMAPIqkIk08gOrtqtrokSiAvytyfnmo6Rg/BvC6ptdvyvpSjvJoLePvAXiyxfs8TFXfj0SUjxaR5rI8poONswC8qKp7Mj6bioke2q8DeCC1eTKSfPrESNX7VPXdAN6KpDdwYfrRIJJufuMm/jyAj6Y9se0YH575JSSDf+NCQzl5C8Y3RJnksKH5u/emYytZ270F7dsJYLKIHN/03skAHsv5+7FyLLG/MuULSVJ3vw7gCiSZNPdI+8HZKuwtbXOVBCX8SMIx+5HcsKek21uQhDJ+t9edquoIgHUAbkEiVI+nHx2CRMxGAEBEfh9tPC4ROVFEFqWe6UtIvOr9vdpk+BiPADhPRCalA2IL2nyvW3k8C2Bm0+uHAPx3Ouj62nT/s0TkVwB8E4k3dpkkg+3non3ICEhCU+1uvIcBvENE+kXkUBG5Fkmv4vNIbtjBRs9CRM4VkePTBucwAG/Aq13ze5rO/VoA31TVRnz8EYyPsZ8C4NsADhaR16Rb19Ti9DtzkIhSYyD91jZf72bDGKp6Zjq2krWdmWHH5LRhmQRgUmr/5HRfLyKJiV8rIoeIyDuRNKQrctgMjC/HXvfXU/mm+5uGZCD3b1T1VlW9A8AnkQz4zsz4SSF767C5akIT/t8DcIuqPqWq329sAD4D4HekXKrm7QDeg6Ywj6p+G8BSJCL1LJIWfWPmr5PY+98CGEUSuvkFJGGNKqnrGB9HEiv/EZIwyF1ZX8pRHtcDuDoN63wyzUY5G8kN8WRq9z8DOFyTFMdzkXjbPwTwW0hutnacgqR3kWXXZiRx8AcA7EHiDCxS1R8jFf6mr/8qkgyO55Hc8H+rqveln30ByUD/AgC/gaQ32eARjPe2T0Yi4D9p2r7dwf4GHwSwriks9WZk1Kl0QLKbDWW4GonNVwI4P/3/6qbP/yeA1wJ4DkkCxGJVbTS8mTY30SjH15bYX6/lCySZSp9Q1WWNN1T1i6kNz1Vgbx02V4oon7lLSG5E5DoAz6nqp2va/4MALlLVRyXJbNoOYLaqvlzH8aomr815y7HXMpAkBXcegIfSsF0pilx3V2zueCwKPyGExEVooR5CCCFdoPATQkhkUPgJISQyKPyEEBIZFH5CCIkMJ5cg7uvr0xkzZtg2gxBCvGHLli2jqjo1z3edFP4ZM2Zg8+bNts0ghBBvEJG2a1m1wlAPIYREBoWfEEIig8JPCCGRQeEnhJDIoPATQkhkUPgJISQyKPyEEBIZRvL4RWQ3kgdb7AfwiqrO7fwLQgghdWFyAtevNT0AnRBCiCWCDvV8dv0wNg2Pb2s2DY/is+uHLVlECCH2MSX8iuRBxltE5GJDx8Ts6Yfj0tu3jYn/puFRXHr7NsyefrgpEwjxDjpM4WNK+N+pqr8M4EwAfyQi72r9gohcLCKbRWTzyMhIJQed39+Hz5w3gEtv34Yb1j6BS2/fhs+cN4D5/X2V7J+Qbvgooq44TD6WnS8YEX5V3Zv+fQ7AVwGcmvGd5ao6V1XnTp2aa4G5XMzv78P5bz8GN943hPPffgxFnxjFFREtgisOk49l5wu1C7+IHCIihzX+B3A6gEfrPm6DTcOjuO3Bp3DZouNw24NPTfAgXCavx0PPyF1cEdGiuOAw+Vp2PmDC4z8SwAMish3AQwDuVtXVBo475iF85rwBXH76iWOVyBfxz+vx0DNyGxdEtCiuOEw+lp0XqKpz25w5c7QKlq0b0o1DI+Pe2zg0osvWDTm536z9bxwa0YFr1+qSldv0pKvvnXDc5uMPXLtWl67ZoQPXrm37PVvUXWYu4/q1aaVhb3MdPOnqe3X5hqEJ36v7+vlWdjYBsFlzaqx1kc/aqhL+usi6MaqslK37W7Jymx57xSpdsnJbx98tXbNDj71ilS5ds6MSO6qk7jJzFR/PO6uRXr5haJzjYeI8fCw7mxQRfkm+7xZz585V15/A1QinnP/2Y3Dbg09VHnts7H/BCX24a9tefGjgaKzfOdL2OHXbUwU+2Fg1n10/jNnTDx93npuGRzG4Zx8+uqDfomXFMX39Qio7E4jIFs27KkLeFsLk5rrH36BuD3vJyq2pp79VVdt7PD55Ri73Skh3eP3cBQU8/qBn7tZJ3YNfm4ZHsfrRZ3HOwDSs3zmKTcOjY1kOg3v2jfvu4J5947yvdt+zjSsDhqQ3eP0CIm8LYXJz3eM3HeOvw4M3PdjqU6+ETITXz31Aj79e6vawTXjwplNAfemVkGx4/cKCg7sRE+Ngawhw0JNkUWRwlx5/BLSb2Tu4Zx8nx3gIJ+xlwxns+aHwR0A7oZh0EDhY5yFcyiAbNoj5YagnElrDOosXzsSydbvGBKN5eYvYBcQXblj7BG68bwiXLToOl59+om1znCDm8CVDPWQCrWue7D8ADtZ5DFMrJ9II6TTX8+b3yauYfPQisUirUGR5QvP7+6LxjnymtXc2r38Ke2tIQj2XrNgCALhs0XG4ZdNu3LJpN266YI5ly9yDHn8E+L5KKRkPUytJWRjjjwCm/5EYaNTzbw3/19jYx7z+KdHU8yIxfgo/ISQYOLjLwV1CSEQwpJkfCj8pDCfKEBfh2Ed+ohd+ilhxOFGmN6qoa6yv7fnogv7MTLUY4vtFiV74KWLF6WXmKAWrmrrG+kqqIHrh5/T33ij6EOwyguVTo9HJ1irqGusrqYLohR8oLmKh0YuwNk8Iu/n+J3Hz/cMTPm/+fRnB8snL7WZrL3Wt9frM7+/DghOmRltfXcQn5wSg8APg9PeiwtqaPXH56cfjurt3jIl/u9/32sD65OV2s7WXutZ6fW6+fxh3bXsa5wwc7XR99U0My+CTcwKAT+Dik4USGue9dM2Oruef9fSu5RuG9KSr7+34+yLHyMKn571m2VqmrjW+u2TlVp1xxSpdvmGo8D5ME9u9VbZ+lwUFnsBlXeSzNpPCb/oRhC5TVlg7/b6sCNi+qYrQztayda1RvktWbu15H6bx6bpVgU3nhMJPClP2Bu32+zKi55Pn2GrblXds11nXrB5nay9C7bOA+tRTK4Pta0ThJ4WoyhuvS5h96pW12rpxaERnXbNar7xj+9jromXjU8PXim0xNIUL14jCTwpRVlh9EmYblBU/X8vXBTE0hQvXqIjwc5E2Qgzgy9OyqlzJlavCmoWLtBHiED6lC1eZlsglFNyFwo+48o2JWXxbMdKnOROkdyj88HDyBfEGH1eMjH0me6/45EBS+EEvxzd8usF8DHf4FJpyCZ8cSAp/iikvxyfRchWfbjDf8C005RI+OZAU/hRTXg5Fqzw+3WBFse0Y+BiacglfwmQUfpj1ckIWLZP4coMVpdUxuOrOQVyyYss4x6BsQ9CpcfExNOUSvoTJKPww7+WEKlom8eUG60bWksuLF87ERbduxg1rn8CqwWfGfb+KHqKJXqftnosNvAqT5Z3pZXILfeZuLNPY6yKkGaHtzmXJym1j69vUUV/qroNVXCMXZsMWwba94JIN7hKSaNnC9g1WNa0ivHzD0ARRrmOhs7oXT6tq4T/eK/koIvxBLtng8lRxl20rg+vn5bp9jSUdzhmYhvU7R8dCj5uGR3HJii0AgN+fPwO3PfhUJWNCjbDE+W8/prJ9ZlF2qQpTdoZA9Es2uJw5E+rgmctlDrhtX/N4xepHn8XihTMn1JGzZh9VWdzYVCy6inEYjofVRN6ugcmtilAP4+jmcb3MXbSvWzijjrCWiVBZVWEaF6+Zq4Ax/oRYHgDhEq6XuWv2hTZe0aCK82KMvxhFhD/IGD/A2KANXC9z1+0j43F9XMY1isT4rXv3WVtZj5+egnlcL3OT9oXqxceAz9cOBTz+IAd3Oe3cPK6XuUn7ehlIjnHCk4u4nARQJcGGegixSdGwUnOmTSONk8t52MHXkGD06ZyE2KZoGiLXcMpP3b2jGFJIKfyE1EAvOewxCE4V1B2OCWUdqE5Mtm0AIa7SmlXy2fXDmHQQsP8AxrJKsrJMWsM08/qn5PLgWwVnXv8Uin8Gzb2jqsMxvV4736DHT0gbWj3LSQcB1929A5PSu6adp9nLQLJXKzs6QF29I9eTFKqCg7ttYA4xASYO9C1eOBPL1u0a52kO7tlXuq6EUt9MnYevA7B14uTgrohMEpFtIrLK1DHLEEtaF+lMq2f5h6f1T/A0q6groazhNHv64bjo1s24+f5koLVRFpMOQmWDr1X1jmJOoTUZ6vk4gMcNHq8URbMsqqxEMVdI12iNu998//CEgb9GXbno1s344y8/MiEtM6brNr+/D5effjyuu3vHWFk0eklVOU1VhWOidu7yzvQqswGYDuA/ACwCsKrb911ajz/v2i5Vzgx1fRZsLLSW+/INQzrjilW6fMNQ5ueNh6csWbkt8/OYaJTFbyzb6HQZVLkInO1Zvygwc9dUVs+nAfwpgMPafUFELgZwMQAcc8wxhszqTJEsiyozDerMWiD5afUs9x8APvWBk7D/QPJ5q6e5fucIzhk4GndtexqAjltXPyY2DY9i/c4RnDrjDXho9w9xzsA0Z8ugOZR32aLjStnZ6EFkTcJzjrwtRK8bgLMA/J/0/4XwxOPv1euucvVH11aSJNlM9Py3pp7/VsuWmadRFo2niC1ZuXVcL8k1ql722eYy0nBsrZ53AvigiOwGsBLAIhG5zcBxS9FrSl5VEz9imEQSCs11JfF2R3HOwDSsfvRZp66bibGjwT37xmL6nzlvAP/4WwP41AdOwg1rv+NUWQC9DxJ3KkdvJuHlbSGq2OCRx18UxviJ69fNlH22Y9156dXOTuXoi8dvNI9fRBYC+KSqntXpey7k8RelyvzlUHK6Y8OH68b892rIKkcAVhfaK5LHzwlchERG2Qegt+JDg1cHreVouxycnMBFiC04L+JV6hg7CjEfvludySpHnybhUfhJ8IQgTFU0XnWtBxTiktKd6kwI6ypR+EnwhCBMVTRedS5A5k02S0461ZkQFnJjjJ9EQ9WxbdO4PDDrsm1l8KnOMMZPSAshzItw1asOIfSRRQh1ph0U/gCIefAyz7mHIkyuClEIoY9WQqkz7YhC+EMXxhAGL3slz7m7Ikxl6qHLQuRTNks7Wq9NYwZyo470Wmec1Z68M71MblXP3C06Y9GXmYfN2JwxaBtfzr3MzFkf66RP1DWr2eRsbhSYuWtd5LO2OpZsKCIOrk+9b0fMi7r5cu6+NFK+UUXDWNe1MXXNiwh/FKEeoNjAmI/pf67Gf03g07m7OkDrO1WEO+u6Ni5e82iEv6g4uHix2uFy/LdufDt3nxopn6jCWavr2jh5zfN2DUxutmP8zd/xoUsec/zXp3P3NYRogqquY68hP8b4HdiqFv6ilYo3KKmDqhspnxq9blRxz5Vx1jqVZZlyNnmNKPwlCemGIuESmoNSRrirKouse3/5hiE96ep7nS/nIsLPJRsI8ZjQlkrodYmEqpZEbl1Dv/G68VQxl8u5yJINph62TgipgSofFm6b1kHQef1Tcp9PlrjP7+8rXB7Ng8StIv/8T14JopyBiLJ6CGng7GzKHnAyY6QHXMrOysroC6WcG1D4SXSEssSFS2JZFleW1QAmNqY33z8cTDk3YIyfRIlrsfFeYtS2H/UXIlkx/otu3YzLTz8e+w9grLwb5Tx7+uHOlDeXZSakC65N0OulFxLC4miukdXz+NyFc8dEv3GN5vf3jb32racI0OMnFeOLF+qax++qTWQ8Ll8jevzEGj7Ez12NjbvWCyETCeUaUfhJpfiwwF23gcS6sn667Te0zJEQCeYa5Z3pZXKzPXOXlMeXZZKzsLFuS2izcEOk12tkaiUAcFlmYhPfvaK6ei3z+/twxtuOxCUrtozbLwAs37Br3DFanwAF+DvXIBR6TTl1MfzJwV1SKe2mvLsW7slDr8sHdGLT8Cj+x60P46WXD+CyRcdh5IWfYtXgM7jpgjlj5bNpeBRf274Xax57NohyJGYGhblkA7FGJ6/IJ8Eqs3xANw6elHS0l9+/CweJYNJBMu64DYE/++RpzmaQkGK4trRGEKGekKbg+04IueV1Zf009nvTBXNw8Wkz8dLLB3BAFZe9+7jMsFIoGSTEvfBnEMLvYgyN+Etdywc09gtgTAQOnnQQdo28mCnwrolFKJh2FJ1MH847Cmxy6yWrx6cnZpF4ycoMmXXNap11zepxdZdZPvVhumxdzOoJanC3jsE4QqqkdWbzpuFRXLJiC86afRSuP3f2mHd4xtuOxNknT3N+BnQdmJj97fIM3F6JcuYuu8XEB1rHQAb37MNNF8zB9efOBvBqWOnYKYd4P1bSKyZCt7GPnwTh8YeUQkiIDVxbY6luj5wefwC4tJY3IT7iWoJEnR65k4OthgnC4yeElMclL7hOW1zr3VRFdB4/8QPOt3AbV+LedXvkIcw1KQuFnxjDtXCCr9TVgLqSIMHQbf0w1EOM4lI4wVfqSGZggoT/MNRDnMWVcILP1LF6KL3suOAibcQodS5+FhNVL/qVFd9urBdEwoMef41wMHM8TKOrDlfi8cRPghV+F0SXg5njYTihGupsQF24b0j9BCv8LoiuD8+fNQnT6KqhzgbUhfuG1E/QWT2uZJBw8bhihDrBxmWay7xx3yw4YSpWP/p9fO7CuVbuG9aDYjCrJ6XqDJJeusGMxRbHttcZY7ijuczn9/dhwQlT8dVtT+N9s95krYdqux6ETNDCX7XoFq2IHMzsDdshMl8Ep8oGqrnM//jL23DXtqdxzsDRWL9zxFp9tV0Pgibvwv0mt14exNJKXQ9bKPLAl6IPYDD1wAZfWLpmhx57xSpdumaH8WP78GCfOur4kpVb9dgrVumSlVsr22dZbNQDH+9FFHgQS7Aef10DYEXCR0UHM33xNE1gO0Tmw0Szqj3iTcOjWP3oszhnYBrW7xwdC/vYzLyyVQ+CvxfzthAmtyo8/rqo2xP0wdOsGxceO1j0Otj0EKvwiF0oc9ds8u1eRAGP37rIZ22uCr+pimgzxOECtrvZvVxnWyJVlTjZLvMsXLDJp3vRKeEH8BoADwHYDuAxAH/Z7TeuCr+Jiuibl2EKkyLQ67FMXzvbHnHo+HYvuib8AuDQ9P+DATwIYF6n37gq/HXDG7k9vpSNSQ/RBY84VHypb80UEX6jE7hE5HUAHgCwWFUfbPe9WJdl5oSVzrgyIa8drttH8uPjvVhkApcR4ReRSQC2ADgOwD+p6hWdvh+r8JPuuDoLOu969j4Iig82kok4N3NXVfer6ikApgM4VURmtX5HRC4Wkc0isnlkZMSEWcQzbKd4diJv+rAPaYI+2EjKYXytHhG5BsCLqvoP7b5Dj5+0EtITonwICflgY6+E2qNxyuMXkakickT6/2sBvAfAjrqPS8IipCWdfZkc5rqNvcIejZlQz1EAviEigwAeBvB1VV1l4LgE4Sw4FtKSzi6HrBr4YGOvcA0gA8KvqoOqOqCqs1V1lqpeW/cxyavQu3GLrIX7Lrp1M26+f3jC92w1zjEsLhhyjyYPwa7VQxJ69W586in4ZGtWyOry04/HDWu/40zjHFJYrR0h92hykTfh3+QWwwQu05Nvik4s8mkCi0+2tsO3WaI+k6e++Dg5Dlyd031MhmB68W58ioP6ZGs7Yg89mCRPjyb4EGneFsLkFoPHr2rGyyvrDfu0SJVPtrZCj989fLsmiN3j9yXma8LLKxOv9SkO6pOtrcQwmOojIffCghR+X7ppJsSq1zRIn8TIJ1uziGEw1Ud8dia6krdrYHKr8tGLrnbTXB+Q9GlwyydbiR+4fn9mAVdX58xLVUs2uLqgFxDutHFCQsDH+9O51TmLUoXwN681cvP9T+Ly04/HH57WP+5zly8iIYQUwam1emzQGvO9/PTjcd3dO8ZmR7oa8yckFHxJsMhLaOcTpPAv37ALixfOHOum/eFp/Tjv7W/G3937hLd53lmEVhljJcTr6EuCRV5CO58ghf/id83EsnW7xl2kex99FmeffFRQqVmhVcZYCfE6hjCprpnQzieKGP9tDz6FxQuTxsDG+uJ1DhSFvG56TIR6HV1OsOgFl88n+hg/MH7yxYITpmLZul3W8rzr9OhCnmQSEyFex9Dy4EM6n2CFv/kirRrcizNnHTlugszihTOxfMMuI7bU2U0MqTLGTGjX0fdJda3UcT5Wx3byJvyb3MpO4GqdbLF8w5DOuGKVLt8wlPm5KapeS8bEJBNOjqofHycLdSO0elPH+VR93VFgApd1kc/aygp/1kVavmFIT7r63spn8uatEHXMJDZxc4UoSq4RmkiShDzXtUpdiF7421HH6o15hNF38XR9+QviHkUas1Abvrz3fVW6VET4g43xt1JXDDVP/N73RbhCHHgk9VIkocF0Oqup2HoebbA2tpO3hTC5Ve3xm/C4W1vtkLwYevzVE1L9aEeRemOyjpnugbfz6G3G+KPw+Ov2uLNa7VAm5YSWneEKodSPThTpKZrsVZqcjNXJo7caCcjbQpjcfHoCV6dWOwRPOQbP1BYh1I9OuOrxN6j7iW2mexbg4K45ugmjz48DJPUTaoiwiOjZSH4w0dCYvpYUJ8RtAAAXeUlEQVQUfkcI3aMj5ciqH75ngDVwOasnlDJupYjwd12rR0T+HcAnVHV7/YGnhKoexGKT5tj4/P6+Ca9J3HSqHwCCXLfHFXx8yEoeKn0Qi4j8MoB/APBdAJ9S1WfKm9iZEIQ/1MpFqqFb/XB5MTDiJrU8gUtEPgzgzwHcCeDvVfUnvZvYmRCEn5BeCXWlTlIvla/OKSIC4AkAywB8DMB3ROSC3k0khGTB9Fligq7CLyIPAHgawD8COBrAhQAWAjhVRJbXaRyZSIhPayKv4vssb+IHeTz+jwI4WlXfq6p/pqqrVHVIVT8G4LSa7SMtxDDxJ2Y+uqB/Qlhnfn+fc+NCdED8pqvwq+qj2n4g4AMV20O6ENoj4Ei1mBJkOiB+U2rJBlU18yQTMg4umkbaYUqQ6YD4TRRr9YSGjRX92LX3g6oEOc/1pgOSDxfvHQq/Z9jK+mDX3h+qEOQ81zu0x0XWhZP3Tt4pviY320s2uLxeik3buASFH1R1nTrtJ9RlD+rCxL0DrtVTjjoqtcuNSRG46JzbVF13213vUOqzSeq+d4oIP0M9GdQxcOVkdw/F4o++dO1djKmaosp5AJ2uty9pp67g3L2Tt4Uwudn2+BtU3UK7GCrJ6yH61LX3yVZXsV2GofQolq0b0uUbhsaV3fINQ3rS1fdWXpagx1+eOlrowT37sOCEqeMG3Wx7onl7Nz7NKGWqYXlsX29Xe8hFmT39cNyw9jtYvHDm2P2+bN0uXH768VbvndyLtJnE9iJtdS2pfPP9w7ju7h340MDRWL9zBIsXzsSydbucEKUQV4MM8ZxioLFyKfDq8tS3bNqNs2YfhevPnW3ZuuKYWnSv8kXaYqMOb6fR0n/qAydh/c4RLDihD9fdvWPME7CJc/HHCgjxnGKh4e0DGEtLfXn/AZx98jTLlvWGi/MdJts2wEWyBqjm9/eVumDNjcnzP3kFN943hHMGpmH/gTKWlqe1NzOvf4r3oZEQzykmGo7WJSu24OX9B/Cagw/CwZP89VFbnZB5/VOs10N/S9MzGlkQzZVg/c5R6zFL27HcOvDxnGLORGrHy/sP4KWXD+Di02bipgvmeLk8tavLbEcp/LZuMhcrQYhpeT6eUyiDmVXxte17cfCkg8a8ZADjGm9fGsp2TsjyDbvs2p83/cfkVnc6Z2tq2pV3bNdZ16yeMDOx6tSxUFLUXCOUcnUx3dcGeVJJbaeblqUO+8GZu91pvslmXbN6nPD7Volix3cRaIYzo/M35L43lFXbX0T4o07nbE73awwA8jmnfhLCc2pDOAfT+J6yW6X9TOfMQetIOwDnUq7Iq3SL6bqYMlcEF8d/XMf3lF2b9kcp/Fk32SUrtuCWTbu9rUSh023w03cRcCUTyZdBU98bSuv2540JmdzqjvG3xhA3Do3orGtW65V3bB977WPMMHTaxURDivHbxpey9H1Avw77wRh/MRpTxJvDA5uGRzG4Z5/TKYAxkhUT5fWrFo41+IlTMX4RebOIfENEHheRx0Tk43Ufsyg+5n3HSLtwDq9ftfg+XtKJMqGsor91OWxmIsb/CoBPqOpbAMwD8Eci8lYDxyUBYT0mGhG+j5d0osxEuaK/dXpSXt6YUFUbgP8H4L2dvuPKevzEHXyP6fpClTF+V69Zmfz5or81OdcArq7HLyIzAAwAeNDkcYn/MJxjhiqzi1z1eMuEsor+1tWwmTHhF5FDAdwBYImq/nfG5xeLyGYR2TwyMmLKLEJIE1U2sK4+EKdMKKvob50Nm+XtGpTZABwMYA2Ay/N8n6EeQsLBpWUoyoSyiv7WdGosXAr1iIgA+ByAx1X1hrqPVyUuj8q7AsuIdMI1j7dMKKvob12ZlJdJ3hai1w3ArwJQAIMAHkm393f6TVUef9nBJV8ms9iEZUTawbphFrjk8avqA6oqqjpbVU9Jt3vqPi5QfnDJpRilq561S2VE3MJpjzdygl6rpwpRcmVU3tUMCcCdMiJuwUwsdwla+IHyouRKjNJlz9qVMiKE5CN44S+buuXSbFGTnnXe0JJrZUSISVwNwXYjaOEvK0quxShNetZ5Q0uulREhJnE5BNuJoFfnDGnVxuZGbH5/34TXdR6TqzQS0h5X7hOnVue0SUiDSzY8aw7aEp8xFYbx8T4JWvhDwkYjxkFb4jOmwjA+3ieTbRtA3KQ1lNR4GD3DPcQXmjPh6grD+Hqf0OMvia+j+t3goC0JgbrDMFXcJzY0hMJfEl9H9bvRHFpqVMzm0FIIjZtNQnUYXKPuMExWCHZwz74J93+na2tDQyj8JXF5YlVVhNq42YRlWj+25pgUvbY2NCTodE6TZD0EPCRcSVkLCZZpvdhM5+7l2pbVEKZzGsbHUf2i+Jiy5jos03qxmc5d9Nqa1hAKf0liWbIghsbNNCzTcClyba1oSN71m01uPj2Bq3XN/2XrhnT5hqFxa/678IDpMsS4rnrdDwqPsUxjoei1raquwaX1+EOntTs5e/rhWLZu19hATh2DdqYzQmJM7ax78DXGMo2FotfWRkiKg7s1UPegnY11e2KEg6/xEMK6XhzctUzdg3YxpJC6AAdf7WOqdxtbem2Qwm97ckzVg3ZZ5wMAb3nTYVGLUt3XmYOv9jElyLE5U0EKv83Wu44R+qzzuWTFFgw+vS9qUarzOseSreU6JgU5ph5esDF+W/HZumKFzedzy6bdAICbLpgTfYy/ruscQsw3JExMkPR9TKdIjN966mbWVlU659I1O/TYK1bp0jU7Cv2u7lS+Xmmcz0eWf9NJ+2zR63UmftBIh1y6ZkdtKa8hpNeC6Zzl4rMuDvQ0n8+O7z8/4fPW9C/b4xymYBw+bEyF3KJLr83bQpjcynr8VbTeJryMorYUOZ8QPJhuxHCOseNq79tFUMDjDzLGX1V81pWF13o9H99jlt1gHJ6QVykS4w9S+KsgFNE00XhRgAmxDydwlcSXVL5ucXxT8W8Xx0QIIe2h8Gfgy0BPJ8E12XjFNvmFEN9hqMdz2oWkbIRfXBkTISRGioR6JtdtDKmX5tmGly06bkzos8R9fn9fbV54a1hpXv8UevyEOApDPZ7jQh67L2MiJCximatSBxR+j3FFcH0ZEykKhcVtmFTQOxR+j3FFcLMeJDG4Z9+EG9A30aSwuI1LSQW+OQkUfo+x+TDpboQgmi4Ji0+YFEFXVtT0rb4HK/y+tcChEYpouiIsVRDiQ01Mj3G1K8NG79uX+h6s8PvWAodICKLpwuB5VYT2UBMbY1ydytCn+h6s8A/u2YfFC2eOq3yLF870fsDRJ3wXTVcGz6sitIea2Bjj6lSGPtX3YIV/9vTDsWzdLiw4YSpuvG8IC06YimXrdln1+It0tX0PVZkUzbrKypXB8yox5ZWaEME8Y1x11I2sMvTNSQhW+Of392Hxwpm4a9vTOHXGG3DXtqexeOFMq92vIl1t30NVJkWzrrJyefC8V0wIsksiWEfdyCpD75yEvOs3m9yqeAJXY232JSu36rFXrNIlK7c6sVZ7kXX+XXomgOuwrLpj6vkFrq2hX2XdcPkZEOATuF6N8a/fOYrLFh2H9TtHnYjxF+lq+zRYZBuWVXdMeaWu9ZSqrBveefZtCFb4GzH+5u6m7Rg/UKyr7dNgkW1YVt1xTZBNwbqRQd6ugcmtilCPa93NxvHzdhNd7lK6BsuKtKPquuFyXUPsj150lSJLJfOpVvlhWZF21FE3XH06Hx+9SKKAgk9s4eKzJ/joRTIB3+cFZNGaqnfVnYO4ZMWWceM4vp8jcY8Qxgwo/JHg+7yALFpnUa4afGbc5yGcI3ELl+YolIGhnohwNTZZluZu97z+KUGeI3EDl8OLDPXkIMTQRzdCzHVv7XYDCO4ciTuEkhIbrfDbCH3YbmxCiE02k9XtvmTFFtyyaXcw50hIHUQr/DbWi7cZZ/cxNtmtoWydRdngrNlHeXOOhNigduEXkc+LyHMi8mjdxyqK6dCHzYeT+DjVvFtD2drtHtyzDzddMAfXnzsbgB/nSIgNah/cFZF3AXgBwBdUdVae35ga3C072NnrQI+LOcCuEuqANCFV49TgrqpuAPCDuo9TlCpCH72EbkKLs9dNiAPShNjGmRi/iFwsIptFZPPIyEjtx6si9FE0dONjnN02bCiJ69hO2ugFZ4RfVZer6lxVnTt16tTS++t2MapKyyrikfoYZ7fJpuFRXHTrZixeOHNcQ3nz/cNO31QkLnycHOmM8FeNqYtRxCMNJQfYFIN79uHy04/HsnW7sGl4dOypajes/Y7TNxWJC5tJG70y2bYBddF8MeoaGGwO3czv7xubNer6Rc/CxRmJjeO+bdrh467j5y6c6135krBp7vlftug45+uniXTOLwH4JoATRWSPiFxU9zEb1D0wGFLoxuXuKgd4iev4NhZVu8evqh+p+xjtaL0Y8/qnVCoaWZ7w/P6+sWO46EW3w0QPqVfqvo6ElMHHnn+wMX4XMmhc9qKzcNGzduE6EtIJH3v+wa7O6Yq37dMEJBdtdeU6EuI6fAKXY/gwU7e1u9r6mhDiNk7N3I0dXwZ9fOyu+kLRCT4+TggifkHhrxGf4tOcY1AfRcd6yo4NseEg3aDw1wi9aAIUn+BTdkKQb0kFxDyM8RNiiKJjPWXGhlwcqHeJEJMGGONPYZeXuELRsZ6yY0Mupua6ROy9oqCFP/aLWxQ2lPVw1Z2DuGTFlgmPiLzqzsHM71cxNuRCUoHL9SlvOM3lcyhD0MLv4+JJNmFD6QZlx4ZcSSpwvT7l6RW5fg69EkWM34c8eldgbLgeTJarS/Frl+tTXttcPodmGONvwoUur0t067oyNlwPJsvVpdRcV+tTkV6Rq+dQhqCF35Uur0t067qyoayHWMvV1fMuEk5z9RxKoarObXPmzNEqWLZuSDcOjYx7b+PQiC5bN1TJ/n1l49CIDly7Vpeu2aED164dK6PG++1ek96ItVxDOG+fzgHAZs2psVHE+MlEssY9XIoNh0Ss5Zp13o1MpuvPnT32nstl4dO14yJtDuJSBfJlsIqEBxcDrA8O7jqIK2lhHPcgNmGKtRtQ+A3hSoV3df2gUCfKkImEmCXjGxR+g7hQ4V1K9WvGlR5RXthQ9U6QWTKeQeE3SLsKTxFxp0eUF98aKldgqNENKPyG6FThKSIJLvSI8uJbQ+UKroYaY4NZPYboltUTQ6ZNiGXA5UCIK0SZ1eN6uKRbbN0nb7dXOvVsfAwBMFZNfCUY4fc9XBKDiHQKj/gWAvCxoSKkQVChHh9DBUB8k1pCCI+4NCGPECDSUA/gb7jEN2+3DKH0bFxNiyUkD5NtG1AlraIyr3+KF+KfJRbz+/u8sL0IrT2Zef1Tgu7ZEOIqwXj8jLm6T0w9m7qpKpnB9aQIUg/BCD9FxX0YHqmOqpIZqk6KYEPiB0EN7hISE1UlM1SZFBFbooJLFBncDSrGT0hMNCczXLbouJ6Ftar9NPbVCLP6ll0XE8GEelzFl66vL3aSV6kqQ6rqTCtfs+tigsJfM75MLPPFTpJQVTJDHUkRoaTshgxj/AbwZWKZL3aS6iaQVT0RjTF+e/DRiz1S52xMX2ar+mIncRPOaLZHtDN3y1JXuMOXrq8vdhJ3YcquHzCrp4k6MhJ8ma3qi52EkPLQ42+h6owEXyaW+WInIaQ8jPG3wAFOQoiPMMbfI1zvhxASA4zxN9EIdzTCG63hDmYmEEJCgB5/E42MhObsnubXLkxmCmWGbSjnQYiPUPgz6PSIwCLUIW6hzLAN5TxIeMTglFD421BFdk8d4lZVo2SbUM4jC9PCYeJ4Js/JtvDG4JRQ+NtQxWSmusQtlEWwQjmPVkwJR0Mgm4+3aXgUV905WPnxTIqhbeEN2SkZQ1Wd2+bMmaM22Tg0ogPXrtWNQyOZr4uydM0OPfaKVbp0zY5K7Vu6Zkcpu2wTynlkYeLcmuvlxqERnXXNaj3x6nt01jWraz2eievlQt2o+r6tGwCbNafG0uPPoMrJTFUvgxBKymko59EOE72ZZs/0W8P/hZf3H8BLLx/A78+f0fF4vYZSTPbQbPcGQ1++hMKfQVXrjdQhbqHMsA3lPNphSjiaBRJAruP1GkoxKYY2hTd0pwQAQz11smzd0IQu6sahEV22bsiSRcQEVYcKux2rNcST53hFQymmz8nUsVQn3qfL1g3p8g1D4+7T5vvW1fsaBUI91kU+awtF+EmcmBKGhiBeecf2CYKf53hFYtgmxc60sBZtaEw3THkpIvxcqydwuD56uJS5tlyTajxFy8PF8nNurR4ReZ+IPCEiQyJypYljuo6pXGXbqXGkPnodi4oihl2QooPJtgefy1K78IvIJAD/BOBMAG8F8BEReWvdx3UdU4IcRU4yKUToA+u9UHQw2fusn7wxoV43AO8AsKbp9VUArur0m1hi/CZzlX3LSSbEFDHG+E2Eeo4G8L2m13vS98YhIheLyGYR2TwyMmLALPuY6i56750QUiNFe0Ah9JhqH9wVkd8AcIaq/kH6+gIAp6rqx9r9JpbBXRMDRK2PVGx9TQgJA9cGd/cAeHPT6+kA9ho4rtOYGmALwTshhFSLCY9/MoCdAN4N4GkADwM4T1Ufa/ebGDx+plkSQqqkiMdf+xO4VPUVEbkUwBoAkwB8vpPox0KWuM/v72P4hRBSO0Yevaiq9wC4x8SxCCGEdIaLtBFCSGRQ+AkhJDIo/IQQEhkUfkIIiQwKPyGERAaFnxBCIoPCTwghkUHhJ4SQyKDwE0JIZDj56EURGQHw3R5/3gcgtnWHec5xwHOOg17P+VhVnZrni04KfxlEZHPehYpCgeccBzznODBxzgz1EEJIZFD4CSEkMkIU/uW2DbAAzzkOeM5xUPs5BxfjJ4QQ0pkQPX5CCCEdCEb4ReR9IvKEiAyJyJW27TGBiHxeRJ4TkUdt22ICEXmziHxDRB4XkcdE5OO2baobEXmNiDwkItvTc/5L2zaZQkQmicg2EVll2xZTiMhuEflPEXlERGp7/mwQoR4RmYTkub7vRfJw94cBfERVv23VsJoRkXcBeAHAF1R1lm176kZEjgJwlKpuFZHDAGwB8KGQr7OICIBDVPUFETkYwAMAPq6q37JsWu2IyOUA5gJ4vaqeZdseE4jIbgBzVbXWuQuhePynAhhS1V2q+jMAKwH8umWbakdVNwD4gW07TKGqz6jq1vT/5wE8DuBou1bViya8kL48ON3899a6ICLTAXwAwD/btiVEQhH+owF8r+n1HgQuCLEjIjMADAB40K4l9ZOGPB4B8ByAr6tq8OcM4NMA/hTAAduGGEYBrBWRLSJycV0HCUX4JeO94L2iWBGRQwHcAWCJqv63bXvqRlX3q+opAKYDOFVEgg7richZAJ5T1S22bbHAO1X1lwGcCeCP0nBu5YQi/HsAvLnp9XQAey3ZQmokjXPfAeCLqnqnbXtMoqo/ArAOwPssm1I37wTwwTTevRLAIhG5za5JZlDVvenf5wB8FUkYu3JCEf6HARwvIr8oIj8H4LcB/Jtlm0jFpAOdnwPwuKreYNseE4jIVBE5Iv3/tQDeA2CHXavqRVWvUtXpqjoDyb18n6qeb9ms2hGRQ9KkBYjIIQBOB1BLxl4Qwq+qrwC4FMAaJAN+/6Kqj9m1qn5E5EsAvgngRBHZIyIX2bapZt4J4AIkHuAj6fZ+20bVzFEAviEig0gcnK+rajTpjZFxJIAHRGQ7gIcA3K2qq+s4UBDpnIQQQvIThMdPCCEkPxR+QgiJDAo/IYREBoWfEEIig8JPCCGRQeEnhJDIoPATQkhkUPgJyUH6HID3pv//tYjcaNsmQnplsm0DCPGEawBcKyK/gGRV0A9atoeQnuHMXUJyIiLrARwKYGH6PABCvIShHkJyICK/hGTdnJ9S9InvUPgJ6UL6yMcvInmq24sicoZlkwgpBYWfkA6IyOsA3AngE6r6OIC/AvAXVo0ipCSM8RNCSGTQ4yeEkMig8BNCSGRQ+AkhJDIo/IQQEhkUfkIIiQwKPyGERAaFnxBCIoPCTwghkfH/AWoXatTU+ODMAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(6, 6))\n", "\n", "plt.plot(x, y, 'x')\n", "\n", "plt.xlabel(r'$x$')\n", "plt.ylabel(r'$y$')\n", "plt.title(r'Arrivals in simulated $Pois(\\lambda L^2)$, $\\lambda={}$, $(0,L] \\times (0,L]$'.format(lambd))\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "\n", "© Blitzstein, Joseph K.; Hwang, Jessica. Introduction to Probability (Chapman & Hall/CRC Texts in Statistical Science)." ] } ], "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.1" } }, "nbformat": 4, "nbformat_minor": 2 }