2a:[] 28:{"title":"Multivariate normal distribution","cells":[{"id":0,"value":"# Multivariate normal distribution","pinned":false,"hidden":false,"mode":"md"},{"id":147,"value":"textdist = md`

${tex`\\boldsymbol{x}\\sim \\mathrm{N}(\\boldsymbol{\\mu},\\boldsymbol{\\Sigma})`}, for ${tex`\\boldsymbol{x} \\in \\mathbb{R}^p`}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":271,"value":"The density function is\n\n${tex`f(\\boldsymbol{x}) = \\vert 2\\pi \\boldsymbol{\\Sigma} \\vert^{-1/2} \\exp\\Big( -\\frac{1}{2}(\\boldsymbol{x}-\\boldsymbol{\\mu})^\\top\\boldsymbol{\\Sigma}^{-1}(\\boldsymbol{x} - \\boldsymbol{\\mu}) \\Big) `} ","pinned":false,"hidden":false,"mode":"md"},{"id":295,"value":"The interactive graph below lets you explore the bivariate case with ${tex`p=2`} variables.","pinned":false,"hidden":false,"mode":"md"},{"id":166,"value":"moments = tex`\n\\begin{aligned}\n\\mathrm{E}( \\boldsymbol{x}) &= \\mu = \n\\begin{pmatrix} \n \\mu_1 \\\\ \n \\mu_2 \n\\end{pmatrix}\n=\n\\begin{pmatrix} \n ${mu[0].toPrecision(2)} \\\\ \n ${mu[1].toPrecision(2)} \n\\end{pmatrix} \\\\[0.5em]\n\\mathrm{Cov}( \\boldsymbol{x}) &=\\Sigma = \n\\begin{pmatrix} \n \\sigma_1^2 & \\rho \\sigma_1 \\sigma_2\\\\ \n \\rho \\sigma_1 \\sigma_2 & \\sigma_2^2 \n\\end{pmatrix}\n=\n\\begin{pmatrix} \n ${(sigma1**2).toPrecision(3)} & ${(rho*sigma1*sigma2).toPrecision(3)} \\\\ \n ${(rho*sigma1*sigma2).toPrecision(3)} & ${(sigma1**2).toPrecision(3)} \n\\end{pmatrix}\n\\end{aligned}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":69,"value":"viewof param = Inputs.form([\n Inputs.range([-1, 1], {label: tex`\\mu_1`, step: 0.1}),\n Inputs.range([-1, 1], {label: tex`\\mu_2`, step: 0.1}),\n Inputs.range([0.1, 2], {label: tex`\\sigma_1`, step: 0.1}),\n Inputs.range([0.1, 2], {label: tex`\\sigma_2`, step: 0.1}),\n Inputs.range([-0.99, 0.99], {label: tex`\\rho`, step: 0.01})\n])","pinned":false,"hidden":false,"mode":"ojs"},{"id":243,"value":"plttrue = Plot.plot({\n width: 500, // Set the width\n height: 500, // Set the height equal to the width for a square aspect\n x: {label: \"x₁\", domain: [-5,5]},\n y: {label: \"x₂\", domain: [-5,5]},\n color: {\n scheme: \"Blues\",\n legend: true,\n type: \"sequential\",\n label: \"pdf\",\n domain: [0,maxpdf]\n },\n marks: [\n Plot.contour(mvnormal_on_grid, {filter: d => d.pdf > Number.EPSILON, x: \"x1\", y: \"x2\", fill: \"pdf\", stroke: \"black\", blur: 4, \n levels: [0.01, 0.02, 0.3]})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\n---\n## 🙈 Appendix 🙈\nHere follow code and details that you probably do not want to read.","pinned":false,"hidden":false,"mode":"md"},{"id":156,"value":"## Implementation of Multivariate normal pdf and sampler","pinned":false,"hidden":false,"mode":"md"},{"id":20,"value":"These are naive implementations of sampling and pdf evaluation of some Multivariate distributions (currently only Normal, but more will come). Gets the job done, at least for reasonably sized and not poorly conditioned problems.","pinned":false,"hidden":false,"mode":"md"},{"id":29,"value":"**Sampling from standard multivariate normal**, similar to `randn` in Matlab/Julia","pinned":false,"hidden":false,"mode":"md"},{"id":7,"value":"function randn(n){\n return d3.range(n).map(d => jstat.normal.sample(0,1))\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":33,"value":"**Sampling from general multivariate normal distribution** ${tex`x \\sim \\mathrm{N}(\\mu,\\Sigma)`}","pinned":false,"hidden":false,"mode":"md"},{"id":9,"value":"function mvnrand(mu, Sigma, n = 1){\n const p = math.matrix(mu).size()[0]\n if (n==1){\n return math.add(mu, math.multiply(jstat.cholesky(Sigma),randn(p))) \n }else{\n return d3.range(n).map(d => math.add(mu, math.multiply(jstat.cholesky(Sigma),randn(p))) )\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":39,"value":"**Evaluating (log)pdf of the multivariate normal distribution** ${tex`x \\sim \\mathrm{N}(\\mu,\\Sigma)`}","pinned":false,"hidden":false,"mode":"md"},{"id":25,"value":"function mvnpdf(x, mu, Sigma, log = false){\n const p = math.matrix(mu).size()[0];\n if (log){\n return -(p/2)*math.log(2*math.pi) -0.5*math.log(math.det(Sigma)) \n -0.5*math.multiply( math.multiply( math.transpose(math.subtract(x,mu)), math.inv(Sigma)), math.subtract(x,mu))\n }else{\n return (2*math.pi)**(-p/2)*math.det(Sigma)**(-0.5)*\n math.exp(-0.5*math.multiply( math.multiply( math.transpose(math.subtract(x,mu)), math.inv(Sigma)), math.subtract(x,mu)))\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":13,"value":"#### Imports","pinned":false,"hidden":false,"mode":"md"},{"id":16,"value":"jstat = require('jstat')","pinned":true,"hidden":false,"mode":"ojs"},{"id":18,"value":"math = require(\"mathjs\");","pinned":true,"hidden":false,"mode":"ojs"},{"id":88,"value":"#### Data","pinned":false,"hidden":false,"mode":"md"},{"id":109,"value":"mu = [param[0],param[1]]","pinned":false,"hidden":false,"mode":"ojs"},{"id":112,"value":"sigma1 = param[2]","pinned":false,"hidden":false,"mode":"ojs"},{"id":114,"value":"sigma2 = param[3]","pinned":false,"hidden":false,"mode":"ojs"},{"id":116,"value":"rho = param[4]","pinned":false,"hidden":false,"mode":"ojs"},{"id":90,"value":"maxpdf = mvnpdf([param[0],param[1]], [param[0],param[1]], [[param[2]**2,param[2]*param[3]*param[4]],[param[2]*param[3]*param[4],param[3]**2]])","pinned":false,"hidden":false,"mode":"ojs"},{"id":246,"value":"mvnormal_on_grid = {\n const ngridx1 = 50;\n const ngridx2 = 50;\n var griddata = [];\n for (let x1 of d3.range(-5, 5, 10/ngridx1)){\n for (let x2 of d3.range(-5, 5, 10/ngridx2)){\n griddata.push({x1: x1, x2: x2, pdf: mvnpdf([x1, x2], [mu[0],mu[1]], [[sigma1**2,rho*sigma1*sigma2],[rho*sigma1*sigma2,sigma2**2]]) })\n }\n }\n return griddata\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":203,"value":"#### Functions","pinned":false,"hidden":false,"mode":"md"},{"id":137,"value":"#### Style","pinned":false,"hidden":false,"mode":"md"},{"id":139,"value":"","pinned":false,"hidden":false,"mode":"md"},{"id":144,"value":"viewof md_colour = Inputs.color({label: \"Highlight color cells\", value: \"#F5F5DC\"})","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W2a","stdlib":"1","titling":"auto"}