{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 1" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m registry at `/opt/julia/registries/General`\n", "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n", "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Qt5Base_jll ─── v5.15.3+1\n", "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Ogg_jll ─────── v1.3.5+1\n", "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m libvorbis_jll ─ v1.3.7+1\n", "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m MKL_jll ─────── v2021.1.1+2\n", "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m GLFW_jll ────── v3.3.5+1\n", "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Glib_jll ────── v2.68.3+2\n", "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Libffi_jll ──── v3.2.2+1\n", "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m StatsBase ───── v0.33.14\n", "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m LERC_jll ────── v3.0.0+1\n", "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Libtiff_jll ─── v4.3.0+1\n", "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.6/Project.toml`\n", " \u001b[90m [2913bbd2] \u001b[39m\u001b[92m+ StatsBase v0.33.14\u001b[39m\n", "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.6/Manifest.toml`\n", " \u001b[90m [2913bbd2] \u001b[39m\u001b[93m↑ StatsBase v0.33.10 ⇒ v0.33.14\u001b[39m\n", " \u001b[90m [0656b61e] \u001b[39m\u001b[93m↑ GLFW_jll v3.3.5+0 ⇒ v3.3.5+1\u001b[39m\n", " \u001b[90m [7746bdde] \u001b[39m\u001b[93m↑ Glib_jll v2.68.3+0 ⇒ v2.68.3+2\u001b[39m\n", " \u001b[90m [88015f11] \u001b[39m\u001b[92m+ LERC_jll v3.0.0+1\u001b[39m\n", " \u001b[90m [e9f186c6] \u001b[39m\u001b[93m↑ Libffi_jll v3.2.2+0 ⇒ v3.2.2+1\u001b[39m\n", " \u001b[90m [89763e89] \u001b[39m\u001b[93m↑ Libtiff_jll v4.3.0+0 ⇒ v4.3.0+1\u001b[39m\n", " \u001b[90m [856f044c] \u001b[39m\u001b[93m↑ MKL_jll v2021.1.1+1 ⇒ v2021.1.1+2\u001b[39m\n", " \u001b[90m [e7412a2a] \u001b[39m\u001b[93m↑ Ogg_jll v1.3.5+0 ⇒ v1.3.5+1\u001b[39m\n", " \u001b[90m [ea2cea3b] \u001b[39m\u001b[93m↑ Qt5Base_jll v5.15.3+0 ⇒ v5.15.3+1\u001b[39m\n", " \u001b[90m [f27f6e37] \u001b[39m\u001b[93m↑ libvorbis_jll v1.3.7+0 ⇒ v1.3.7+1\u001b[39m\n", "\u001b[32m\u001b[1mPrecompiling\u001b[22m\u001b[39m project...\n", "\u001b[32m ✓ \u001b[39m\u001b[90mLibffi_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mLERC_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mMKL_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mOgg_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mGLFW_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mGlib_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mWayland_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mLibtiff_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mlibvorbis_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39mStatsBase\n", "\u001b[32m ✓ \u001b[39m\u001b[90mWayland_protocols_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mCoupledFields\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mCairo_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mFFMPEG_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mxkbcommon_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mFFMPEG\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mQt5Base_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mGR_jll\u001b[39m\n", "\u001b[32m ✓ \u001b[39mDistributions\n", "\u001b[32m ✓ \u001b[39m\u001b[90mFFTW\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mGR\u001b[39m\n", "\u001b[32m ✓ \u001b[39m\u001b[90mKernelDensity\u001b[39m\n", "\u001b[32m ✓ \u001b[39mGadfly\n", "\u001b[32m ✓ \u001b[39mPlots\n", " 24 dependencies successfully precompiled in 42 seconds (190 already precompiled)\n" ] } ], "source": [ "using Pkg\n", "Pkg.add(\"StatsBase\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Generating exponentially distributed random variables\n", "\n", "Let us simulate $10^5$ realizations of a random variable $X$ that follows an Exponential distribution of parameter $\\lambda = 10$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- with the uniform distribution" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "λ = 10\n", "g = y -> -log(y)/λ # write what your function should return here\n", "u = rand(100000) # draw 10^5 realizations from a uniform distribution with rand()e\n", "g_eval = g.(u); # evaluate your function g at those realizations" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Plots\n", "histogram(g_eval, normalize = true, xlim = (0, 1), label = \"from Uniform\")\n", "#histogram(u, normalize = true, xlim = (0, 1), label = \"from Uniform\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- with the exponential distribution\n", "\n", "[1]https://github.com/JuliaStats/Distributions.jl/tree/ebeab79ff28f144506f6aa51b284b67486283ba0/src/univariate/continuous\n", "[2]https://github.com/JuliaLang/julia/blob/master/stdlib/Random/src/normal.jl" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Distributions\n", "pdf_evaluation = x -> pdf.(Exponential(1/λ), x)\n", "plot!(0:0.01:2.0, pdf_evaluation, \n", " linewidth = 3.0, \n", " label = \"from Exponential\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Gillespie algorithm (naive implementation)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "κ = 0.01 # rate constant\n", "x₀ = 20 # initial count of molecules\n", "Δt = 0.001; " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "gillespie_alg_naive = function(x₀, κ, Δt, stoptime)\n", " \n", "\n", " tt = zeros(Integer(round(stoptime/Δt))+1)\n", " xx = Integer.(zeros(length(tt)))\n", " \n", " t = 0\n", " k = 0\n", " \n", " xx[1] = x₀ # initial count\n", " tt[1] = t # initial time\n", " \n", " k = 1\n", " while (t < stoptime)\n", " x = xx[k]\n", " \n", " # step 1\n", " ν = x*κ*Δt\n", " # step 2\n", " u = rand()\n", " # step 3\n", " u < ν ? xx[k+1] = x - 1 : xx[k+1] = x\n", " \n", " t = k*Δt\n", " tt[k+1] = t\n", " k = k + 1 \n", " end\n", " \n", " return tt, xx\n", "end;" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "tt, xx = gillespie_alg_naive(x₀, κ, Δt, 30);" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(tt, xx, linetype=:steppre, xlabel = \"time (sec)\", ylabel = \"# molecules of A\", legend=false)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Gillespie algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Preliminary steps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us consider the reactions\n", "\n", "$A \\xrightarrow{\\kappa_1} \\emptyset$\n", "\n", "$\\emptyset \\xrightarrow{\\kappa_2}A $\n", "\n", "\\begin{align}\n", "(\\text{Reactants})\\,\\underline{S} = \\begin{bmatrix}\n", "n_{A,\\,r_1}^{(\\text{reac})}\\\\\n", "n_{A,\\,r_2}^{(\\text{reac})}\\\\\n", "\\end{bmatrix}_{\\text{#reactions }\\times \\text{ #species}} &= \\begin{bmatrix}\n", "1\\\\\n", "0\\\\\n", "\\end{bmatrix}\\\\\n", "\\hspace{2cm}\n", "(\\text{Products})\\,\\overline{S} = \\begin{bmatrix}\n", "n_{A,\\,r_1}^{(\\text{prod})}\\\\\n", "n_{A,\\,r_2}^{(\\text{prod})}\\\\\n", "\\end{bmatrix}_{\\text{#reactions }\\times \\text{ #species}} &= \\begin{bmatrix}\n", "0\\\\\n", "1\\\\\n", "\\end{bmatrix}\n", "\\end{align}\n", "\n", "\\begin{align}\n", "(\\text{Stoichiometric Matrix})\\, S = (\\overline{S} - \\underline{S})^\\top = \\begin{bmatrix}\n", "-1 & 1\n", "\\end{bmatrix}\n", "\\end{align}\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×2 adjoint(::Vector{Int64}) with eltype Int64:\n", " -1 1" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Number of reactions\n", "Nᵣ = 2\n", "\n", "# Compute the stoichiometric matrix \n", "Products = [0; 1]\n", "Reactants = [1; 0]\n", "SM = (Products - Reactants)'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(step 1) Compute the propensity function $\\nu_{n_r}(t)$ for each of the $N_r$ reactions. Then compute the total propensity $\\alpha(t)$." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "#9 (generic function with 1 method)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Compute the propensity function for each of the reactions\n", "\n", "V = 0.001 # Volume \n", "ν₁ = function(κ₁, x) return κ₁*x[1] end # Propensity function for degradation\n", "ν₂ = function(κ₂, x) return κ₂*V end # Propensity function for production" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "#11 (generic function with 1 method)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Compute the total propensity\n", "propensities = function(k, x, ν, Nᵣ)\n", " ν_eval = zeros(Nᵣ)\n", " for r in 1:Nᵣ\n", " ν_eval[r] = ν[r](k[r], x) \n", " end\n", " return ν_eval\n", "end" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 0.0\n", " 1.0" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Example : \n", "k = [0.1, 1000]\n", "x = [0]\n", "ν = [ν₁, ν₂]\n", "\n", "ν_eval = propensities(k, x, ν, 2) # Evaluate propensity functions " ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.0" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "α = sum(ν_eval)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(step 2) Compute the time when the first chemical reaction takes place as $\\tau$, where $\\tau$ is sampled from an exponential distribution with parameter $\\alpha(t)$" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.4397891544092207" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "τ = rand(Exponential(1/α)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(step 3) Find out which reaction occurs at time $\\tau$" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Reaction n.2 takes place!\n" ] } ], "source": [ "using StatsBase\n", "index_j = StatsBase.sample(1:2, Weights(ν_eval./α));\n", "println(\"Reaction n.\" .* string(index_j) .* \" takes place!\")" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Update number of molecules A with the stoichiometric vector: [1]\n" ] } ], "source": [ "println(\"Update number of molecules A with the stoichiometric vector: \" .* string(SM[:, index_j]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gillespie algorithm (main algorithm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Combine all these steps in a loop." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "#13 (generic function with 1 method)" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gillespie_alg = function(SM, κ, ν, x₀, stoptime)\n", " t = 0.0\n", " x = x₀\n", " \n", " Nₛ, Nᵣ = size(SM)\n", " \n", " tt = [t]\n", " xx = copy(x) \n", " \n", " k=1\n", " \n", " while (t <= stoptime)\n", " \n", " # step 1\n", " ν_eval = propensities(κ, x, ν, Nᵣ) \n", " α = sum(ν_eval) \n", "\n", " # step 2\n", " τ = rand(Exponential(1/α)) \n", " t += τ\n", "\n", " # step 3\n", " index_j = StatsBase.sample(1:Nᵣ, Weights(ν_eval./α))\n", " \n", " x = x + SM[:, index_j]\n", " append!(tt, t)\n", " xx = cat(xx, x, dims=2)\n", " end\n", " \n", " return tt, xx'\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (a) Stochastic simulation of degradation" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×1 adjoint(::Vector{Int64}) with eltype Int64:\n", " -1" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Number of reactions Nᵣ = 1\n", "# Compute the stoichiometric matrix \n", "Products = [0]\n", "Reactants = [1]\n", "SM = (Products-Reactants)'" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "#15 (generic function with 1 method)" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Compute the propensity function for each of the reactions\n", "V = 0.001\n", "ν₁ = function(κ₁, x) return κ₁*x[1] end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1 realization" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tt, xx = gillespie_alg(SM, [0.1], [ν₁], [20], 100.0)\n", "plot(tt[1:end-1], xx[1:end-1], linetype=:steppre, xlabel=\"time (sec)\", ylabel =\"# molecules of A\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "20 realizations" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Nₛ = 20\n", "p = plot(legend = false, xlabel = \"time (sec)\", ylabel = \"# molecules of A\")\n", "for nₛ in 1:Nₛ \n", " tt, xx = gillespie_alg(SM, [0.1], [ν₁], [20], 100.0)\n", " plot!(tt[1:end-1], xx[1:end-1], linetype=:steppre)\n", "end\n", "p" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (b) stochastic simulation of production and degradation" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×2 adjoint(::Vector{Int64}) with eltype Int64:\n", " -1 1" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Number of reactions Nᵣ = 2\n", "# Compute the stoichiometric matrix \n", "Products = [0; 1]\n", "Reactants = [1; 0]\n", "SM = (Products-Reactants)'" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "#19 (generic function with 1 method)" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Compute the propensity function for each of the reactions\n", "V = 0.001\n", "ν₁ = function(κ₁, x) return κ₁*x[1] end\n", "ν₂ = function(κ₂, x) return κ₂*V end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1 realization" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "κ₁ = 0.1\n", "κ₂ = 1000.0\n", "\n", "tt, xx = gillespie_alg(SM, [κ₁, κ₂], [ν₁, ν₂], [0], 30.0)\n", "plot(tt, xx, linetype=:steppre, label = \"1 realization of Production + Degradation\", legend=:bottomright)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "10 realizations" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Nₛ = 10\n", "p = plot(legend=false, xlabel=\"time (sec)\", ylabel =\"# molecules of A\")\n", "for nₛ in 1:Nₛ \n", " tt, xx = gillespie_alg(SM, [κ₁, κ₂], [ν₁, ν₂], [0], 30.0)\n", " plot!(tt[1:end-1], xx[1:end-1], linetype=:steppre)\n", "end\n", "p" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Evolution of the stochastic mean $M(t)$" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M₀ = 0\n", "M = t -> κ₂*V/κ₁ + (M₀ - κ₂*V/κ₁)*exp(-κ₁*t)\n", "plot!(tt[1:end-1], M, linestyle = :dash, colour = :black, linewidth = 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (c) Steady-state version of the chemical master equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Number of molecules at $t=100$sec, from $10^5$ realizations" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "N = 1e5 #number of realisations\n", "n_upperbound = 24\n", "p = zeros(n_upperbound+1) #vector of probabilities\n", "\n", "\n", "for n in 1:N \n", " tt, xx = gillespie_alg(SM, [κ₁, κ₂], [ν₁, ν₂], [0], 50.0)\n", " A = xx[end] \n", " if (A <= n_upperbound) \n", " p[xx[end]+1] += 1\n", " end\n", "end\n", "\n", "p /= N;" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bar(collect(0:24), p, \n", " colour = :lightgrey,\n", " label = \"Gillespie SSA\",\n", " xlabel = \"number of molecules\",\n", " ylabel = \"stationary distribution\")" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "#25 (generic function with 1 method)" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "compute_ϕ = function(n_upperbound, κ₁, κ₂)\n", " \n", " ϕ = zeros(n_upperbound+1)\n", "\n", " ϕ[1] = 1.0\n", " ϕ[2] = κ₂*V/κ₁ * ϕ[1];\n", "\n", " for n in 1:n_upperbound-1\n", " ϕ[n+2] = (κ₁*n*ϕ[n+1] + κ₂*V*ϕ[n+1] - κ₂*V*ϕ[n])/(κ₁*(n+1));\n", " end\n", " return ϕ./(sum(ϕ))\n", "end" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "ϕ = compute_ϕ(n_upperbound, 0.1, 1000.0);" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot!(collect(0:n_upperbound), ϕ,\n", " linewidth = 3,\n", " colour = :red,\n", " label = \"master equation\")" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.6.2", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.2" } }, "nbformat": 4, "nbformat_minor": 4 }