2b:null 2c:["$","title",null,{"children":"Multivariate Normal Distribution | sw1227 | Observable"}] 31:T496,covMatChart = { const svg = d3.create("svg") .attr("width", width) .attr("height", 200); const size = 100; const xOffset = 50; const maxAbs = Math.max(...math.flatten(A).toArray().map(Math.abs)) const color = d => d3.interpolateRdBu((d+maxAbs)/(2*maxAbs)); svg.append("text") .attr("x", xOffset-10).attr("y", size) .attr("text-anchor", "end") .attr("dominant-baseline", "central") .text("A = "); const g = svg.selectAll(".row") .data(A.toArray()) .enter().append("g") .attr("transform", (d, i) => `translate(${xOffset}, ${i*size})`) .selectAll(".cell") .data(d => d) .enter(); g.append("rect") .attr("x", (d, i) => i*size) .attr("y", 0) .attr("width", size).attr("height", size) .attr("fill", color); g.append("rect") .attr("x", (d, i) => i*size + size/4) .attr("y", size/4) .attr("width", size/2).attr("height", size/2) .attr("fill", "#fff") .attr("opacity", 0.5); g.append("text") .attr("x", (d, i) => i*size + size/2).attr("y", size/2) .attr("text-anchor", "middle") .attr("dominant-baseline", "central") .text(d => d); return svg.node(); }32:[] 2f:{"id":"e2b0cf7f5c18f3e8","slug":"multivariate-normal-distribution","title":"Multivariate Normal Distribution","thumbnail":"ebbd7f51fba6d5d3aa0ee52dc9ca238660c8954e2ffd600dea8d464398abb25f","access_level":"public","create_time":"$D2019-08-05T11:21:19.384Z","update_time":"$D2020-11-28T06:31:51.488Z","view_time":null,"first_public_time":"$D2019-08-25T04:50:48.916Z","fork":false,"stars":20,"starred":false,"owner":{"type":"team","id":"06475ae5027cde36","login":"sw1227","name":"sw1227","avatar_url":"https://avatars-next.observableusercontent.com/3e09ff2c49317e618eedc394f3b33dbb5c952acf0d9654683424e681e7857f8c","tier":"starter_2024"},"version":925,"publish_version":925,"role":null,"authored":false,"authors":[{"type":"individual","id":"4a648b0bd17e2c08","login":"sw1227","name":"sw1227","avatar_url":"https://avatars-next.observableusercontent.com/3e09ff2c49317e618eedc394f3b33dbb5c952acf0d9654683424e681e7857f8c","tier":null}],"members":[{"type":"individual","id":"4a648b0bd17e2c08","login":"sw1227","name":"sw1227","avatar_url":"https://avatars-next.observableusercontent.com/3e09ff2c49317e618eedc394f3b33dbb5c952acf0d9654683424e681e7857f8c","tier":null,"role":"editor"}],"comments":[{"id":"b86f41aded53be97","reply_to_id":null,"content":"Thank you for your work. Just correcting typos: Covariande -> Covariance and definete -> definite","create_time":"$D2024-07-08T03:31:10.640Z","update_time":null,"resolved":false,"user":{"type":"individual","id":"fc88ff3942f7785e","login":"josephsmann","name":"Joe","avatar_url":"https://avatars-next.observableusercontent.com/8089eca3d0d1bfc18ba5f7492ea1574d05db45eebbc8e0f3af1c72105b236d4c","tier":null}},{"id":"7ad361df776c27c6","reply_to_id":null,"content":"shouldn't \nmath.multiply(cov, math.matrix(z)) be \nmath.multiply(L, math.matrix(z)) ?","create_time":"$D2024-07-08T03:41:57.272Z","update_time":"$D2024-07-08T03:44:23.544Z","resolved":false,"user":{"type":"individual","id":"fc88ff3942f7785e","login":"josephsmann","name":"Joe","avatar_url":"https://avatars-next.observableusercontent.com/8089eca3d0d1bfc18ba5f7492ea1574d05db45eebbc8e0f3af1c72105b236d4c","tier":null}}],"custom_thumbnail":true,"trashed":false,"subscribed":null,"forks":5,"collection":{"id":"c640659921d068e9","slug":"statistics","title":"Statistics","owner":{"id":"06475ae5027cde36","login":"sw1227","name":"sw1227"}},"type":"notebook","body":{"title":"Multivariate Normal Distribution","cells":[{"id":0,"value":"md`# Multivariate Normal Distribution\nImplementation and visualization of **multivariate** normal distribution.\n\n## 1. Definitions & properties\nProbability Density Fucntion (PDF) of ${tex`n`}-variate normal distribution is:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":225,"value":"tex.block`\nf(\\bm{x}) = \\frac{1}{(\\sqrt{2 \\pi})^n \\sqrt{|\\Sigma|}} \\exp \n\\bigl( -\\frac{1}{2} (\\bm{x} - \\bm{\\mu})^T \\bm{\\Sigma}^{-1} (\\bm{x} - \\bm{\\mu}) \\bigr)\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":255,"value":"md`\nwhere ${tex`\\bm{\\mu}`} is ${tex`n`}-dimensional mean vector and ${tex`|\\bm{\\Sigma}|`} is the determinant of ${tex`n \\times n`} covariance matrix ${tex`\\bm{\\Sigma}`}. Covariande matrix ${tex`\\bm{\\Sigma}`} is positive definete and symmetric.\n\nThe differential entropy of this distiribution is:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":206,"value":"tex.block`\nH(\\bm{x}) = \\frac{1}{2} \\ln | \\boldsymbol{\\Sigma}| + \\frac{n}{2}(1+\\ln(2\\pi))\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":368,"value":"md`## 2. Sampling from normal distribution\n\n### (1) Sampling by Box-Muller transform (1D)\nTwo independent random variables (${tex`Z_0, Z_1`}) with a standard normal distribution can be obtained by the following equations:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":388,"value":"tex.block`\n \\left\\{\n \\begin{array}{l}\n Z_0 = R \\operatorname{cos}(\\theta) = \\sqrt{-2\\operatorname{ln}U_1} \\operatorname{cos}(2 \\pi U_2) \\\\\n Z_1 = R \\operatorname{sin}(\\theta) = \\sqrt{-2\\operatorname{ln}U_1} \\operatorname{sin}(2 \\pi U_2)\n \\end{array}\n \\right.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":391,"value":"md`\nwhere ${tex`U_1, U_2`} are independent samples from the uniform distribution on ${tex`(0, 1)`}.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":541,"value":"md`Box-Muller transform in implemented in [my previous notebook: Box-Muller transform](https://observablehq.com/@sw1227/box-muller-transform).`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":546,"value":"import {boxMuller} from \"@sw1227/box-muller-transform\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":382,"value":"md`### (2) Sampling from multivariate normal distribution\nBelow is the widely used method [from Wikipedia](https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution).\n\nFirst, find a matrix ${tex`L`} such that ${tex`\\Sigma = LL^T`}. Since covariance matrix ${tex`\\Sigma`} is positive definite and symmetric, we can use [Cholesky decomposition](https://observablehq.com/@sw1227/cholesky-decomposition) of ${tex`\\Sigma`} to obtain a lower triangular matrix ${tex`L`} that meets the equation.\n\nLet ${tex`\\bm{z} = (z_1, z_2, ..., z_n)^T`} be samples from ${tex`n`} independent standard normal variates (which can be generated by the Box-Muller transform as mentioned above). Samples ${tex`\\bm{x}`} from multivariate normal distribution ${tex`N( \\bm{\\mu}, \\bm{\\Sigma})`} can be calculate by:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":410,"value":"tex.block`\n\\bm{x} = \\bm{\\mu} + \\bm{L} \\bm{z}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":417,"value":"md`\nWe can check ${tex`\\bm{x} \\sim N( \\bm{\\mu}, \\bm{\\Sigma})`} by 2 steps. Note that ${tex`\\bm{z} \\sim N(\\bm{0}, \\bm{I})`}.\n\n- ${tex`E(\\bm{x}) = E(\\bm{\\mu}) + E(\\bm{Lz}) = \\bm{\\mu} + \\bm{L} E(\\bm{z}) = \\bm{\\mu}`} \n- ${tex`E[(\\bm{x} - E(\\bm{x}))(\\bm{x} - E(\\bm{x}))^T] = E(\\bm{Lzz^TL^T}) = \\bm{L}E(\\bm{zz^T}) \\bm{L^T} = \\bm{L I L^T} = \\bm{LL^T} = \\bm{\\Sigma}`} (Covariance matrix).\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":554,"value":"md`Cholesky decomposition is implemented in [my previous notebook: Cholesky Decomposition](https://observablehq.com/@sw1227/cholesky-decomposition).`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":557,"value":"import {choleskyDecomposition} from \"@sw1227/cholesky-decomposition\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":289,"value":"md`## 3. Implementation\nImplemement a function to generate a normal distribution object that can:\n- Calculate probability density at any point\n- Calculate the entropy of the distribution\n- Generate samples from the distribution\n\n[math.js](https://mathjs.org/docs/datatypes/matrices.html) is used for matrices operation.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":194,"value":"function multivariateNormal(mean, covArray) {\n const n = mean.length;\n const cov = math.matrix(covArray);\n return {\n // Probability Density Function\n pdf: x => {\n const c = 1 / (math.sqrt(2*math.PI)**n * math.sqrt(math.det(cov)));\n return c * math.exp(\n -(1/2) * math.multiply(\n math.subtract(math.matrix(x), math.matrix(mean)),\n math.inv(cov),\n math.subtract(math.matrix(x), math.matrix(mean))\n )\n );\n },\n // Differential entropy\n entropy: 0.5*math.log(math.det(cov)) + 0.5*n*(1 + math.log(2*math.PI)),\n // Generate n samples using Cholesky Decomposition\n sample: n_samples => Array(n_samples).fill().map(_ => {\n const L = choleskyDecomposition(cov);\n const z = boxMuller(n);\n return math.add(\n math.matrix(mean),\n math.multiply(cov, math.matrix(z))\n ).toArray();\n }),\n };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":326,"value":"// Usage example\n{\n const norm = multivariateNormal([0, 0], [[1, 0], [0, 2]]);\n return [norm.entropy, norm.pdf([0, 0]), norm.sample(3)]\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":594,"value":"md`## 4. Visualization\nIn this section, we will visualize 2D normal distribution. \nFirst, let's create a slider to interactively change the value of covariance matrix ${tex`\\bm{A}`}.\nAll the visualizations below changes reactively as the slider value changes.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":733,"value":"tex.block`\n\\bm{A} = \n\\begin{pmatrix}\n{\\sigma_1}^2 & \\sigma_{12} \\\\\n\\sigma_{12} & {\\sigma_2}^2\n\\end{pmatrix}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":745,"value":"md`\nSince scaling all the element in ${tex`\\bm{A}`} does not change the shape of distribution, let's fix ${tex`{\\sigma_1}^2 = 1`}.\nNow we just need two sliders for ${tex`\\sigma_{12}`} and ${tex`{\\sigma_2}^2`} (because ${tex`\\bm{A}`} is symmetric).\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":683,"value":"md`Slider for ${tex`\\sigma_{12}`}: `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":604,"value":"viewof cov12 = {\n const form = html`\n`;\n form.i.oninput = () => form.o.value = `${form.value = form.i.valueAsNumber}`;\n form.i.oninput();\n return form;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":689,"value":"md`Slder for ${tex`{\\sigma_2}^2`}: `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":608,"value":"viewof v2 = {\n const form = html`\n`;\n form.i.oninput = () => form.o.value = `${form.value = form.i.valueAsNumber}`;\n form.i.oninput();\n return form;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":624,"value":"$31","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":915,"value":"md`Here, covariance matrix ${tex`A`} and normal distribution object is instantiated based on the slider value.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":912,"value":"A = math.matrix([[1, cov12], [cov12, v2]])","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":823,"value":"// mean: [0, 0]\nnormal = multivariateNormal([0, 0], A)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":831,"value":"md`### 4.1. Scatter plot and density contours of samples\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":580,"value":"samples2D = normal.sample(1000).map(d => ({x: d[0], y: d[1]}))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":578,"value":"import {chart as densityChart} with {samples2D as data} from \"@sw1227/fork-density-contours\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":587,"value":"densityChart","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":841,"value":"md`### 4.2. Probability distribution function\n#### (1) 2D Visualization\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":845,"value":"// Values of PDF in x: (-5, 5), y: (-5, 5)\ndata = d3.range(-5, 5, 0.1).map(y => (\n d3.range(-5, 5, 0.1).map(x => (\n normal.pdf([x, y])\n ))\n))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":59,"value":"imshow(data.reverse(), Math.min(width/100, 4), d3.interpolateViridis)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":116,"value":"md`#### (2) 3D Visualization`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":97,"value":"{\n const renderer = new THREE.WebGLRenderer({antialias: true});\n const controls = new THREE.OrbitControls(camera, renderer.domElement);\n invalidation.then(() => (controls.dispose(), renderer.dispose()));\n renderer.setSize(width, height);\n renderer.setPixelRatio(devicePixelRatio);\n controls.addEventListener(\"change\", () => renderer.render(scene, camera));\n renderer.render(scene, camera);\n return renderer.domElement;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":91,"value":"wireframe = (size, shape) => {\n // Geometry\n const wireframeGeometry = new THREE.PlaneGeometry(\n size.width, size.height, // width, height\n shape[0], shape[1] // Segments\n );\n // Material\n const wireframeMaterial = new THREE.MeshBasicMaterial({\n color: 0x2260ff,\n wireframe: true,\n transparent: true,\n depthWrite: false,\n side: THREE.DoubleSide,\n blending: THREE.AdditiveBlending\n });\n\n // Mesh\n const mesh = new THREE.Mesh(wireframeGeometry, wireframeMaterial);\n\n const heightArray = [].concat.apply([], data);\n // Update\n mesh.geometry.vertices.forEach((v, i) => {\n v.setZ(heightArray[i]);\n });\n mesh.geometry.verticesNeedUpdate = true;\n mesh.geometry.computeFaceNormals(); \n return mesh;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":87,"value":"scene = {\n const scene = new THREE.Scene();\n scene.background = new THREE.Color(0x333333);\n scene.add(new THREE.AxisHelper(10));\n scene.add(wireframe({width: 1, height: 1}, [99, 99]));\n return scene;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":84,"value":"camera = {\n const fov = 45;\n const aspect = width / height;\n const near = 0.1;\n const far = 1000;\n const camera = new THREE.PerspectiveCamera(fov, aspect, near, far);\n camera.position.set(0.3, -0.7, 0.2);\n camera.up.set(0, 0, 1);\n camera.lookAt(new THREE.Vector3(0, 0, 0));\n return camera;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":78,"value":"height = 500","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":61,"value":"import {imshow} from \"@sw1227/reusable-2d-array-image-function\";","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":74,"value":"THREE = {\n const THREE = window.THREE = await require(\"three@0.99.0/build/three.min.js\");\n await require(\"three@0.99.0/examples/js/controls/OrbitControls.js\").catch(() => {});\n return THREE;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":55,"value":"d3 = require(\"d3@5\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":4,"value":"math = require(\"mathjs@7\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W32","stdlib":"1","titling":"auto"}} 30:"$2f:body"