2b:null 2c:["$","title",null,{"children":"Multivariate Normal Distribution | Joe | Observable"}] 31:T424,function multivariateNormal(mean, covArray) { const n = mean.length; const cov = math.matrix(covArray); return { // Probability Density Function pdf: (x) => { const c = 1 / (math.sqrt(2 * math.PI) ** n * math.sqrt(math.det(cov))); return ( c * math.exp( -(1 / 2) * math.multiply( math.subtract(math.matrix(x), math.matrix(mean)), math.inv(cov), math.subtract(math.matrix(x), math.matrix(mean)) ) ) ); }, // Differential entropy entropy: 0.5 * math.log(math.det(cov)) + 0.5 * n * (1 + math.log(2 * math.PI)), // Generate n samples using Cholesky Decomposition sample: (n_samples) => Array(n_samples) .fill() .map((_) => { const L = cholesky2(covArray); // const L = choleskyDecomposition(cov); const z = boxMuller(n); return math .add(math.matrix(mean), math.multiply(L, math.matrix(z))) .toArray(); }) }; }32: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(); }33:[] 2f:{"id":"6298f34a299294a7","slug":"multivariate-normal-distribution2","title":"Multivariate Normal Distribution","thumbnail":"49e3b8e7b5c25c76853ac31b4ed41e895e2facb91573182f4e68e995d0ccb899","access_level":"public","create_time":"$D2024-07-05T17:23:44.214Z","update_time":"$D2024-09-11T14:59:56.354Z","view_time":null,"first_public_time":"$D2024-07-05T17:23:44.214Z","fork":{"id":"e2b0cf7f5c18f3e8","slug":"multivariate-normal-distribution","title":"Multivariate Normal Distribution","owner":{"id":"06475ae5027cde36","login":"sw1227","name":"sw1227"}},"stars":0,"starred":false,"owner":{"type":"team","id":"d4bde6b3ee650269","login":"josephsmann","name":"Joe","avatar_url":"https://avatars-next.observableusercontent.com/8089eca3d0d1bfc18ba5f7492ea1574d05db45eebbc8e0f3af1c72105b236d4c","tier":"starter_2024"},"version":937,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"fc88ff3942f7785e","login":"josephsmann","name":"Joe","avatar_url":"https://avatars-next.observableusercontent.com/8089eca3d0d1bfc18ba5f7492ea1574d05db45eebbc8e0f3af1c72105b236d4c","tier":null}],"members":[{"type":"team","id":"d4bde6b3ee650269","login":"josephsmann","name":"Joe","avatar_url":"https://avatars-next.observableusercontent.com/8089eca3d0d1bfc18ba5f7492ea1574d05db45eebbc8e0f3af1c72105b236d4c","tier":"starter_2024","role":"editor"},{"type":"individual","id":"fc88ff3942f7785e","login":"josephsmann","name":"Joe","avatar_url":"https://avatars-next.observableusercontent.com/8089eca3d0d1bfc18ba5f7492ea1574d05db45eebbc8e0f3af1c72105b236d4c","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":0,"collection":null,"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":927,"value":"import { cholesky2 } from \"@josephsmann/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":"$31","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":326,"value":"// Usage example\n{\n const norm = multivariateNormal(\n [0, 0],\n [\n [1, 0],\n [0, 2]\n ]\n );\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":"$32","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(100) //.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":"$W33","stdlib":"1","titling":"auto"}} 30:"$2f:body"