34:null 35:["$","title",null,{"children":"Finite Mixtures of Gaussian Linear Regressions | Dylan Green | Observable"}] 3b:T932, ### An Example Let's work an example that demonstrates how this is useful. Let's assume that ${tex`Y`} is drawn from a distribution such that ${tex.block` y_i = \boldsymbol{x_i} \boldsymbol{\theta} + \epsilon_i`} and for this example let's assume that ${tex`\boldsymbol{x_i}`} is one dimensional (${tex`\boldsymbol{x_i} = x_i`}), and that we're fitting a line: ${tex.block` y_i = \theta_1 x_i + \theta_2 + \epsilon_i`} Evidently ${tex`\theta_1`} is the slope and ${tex`\theta_2`} will be the intercept of the fitted model, and ${tex`\epsilon_i`} represents some noise around the line (so not all points lie directly on the line we fit, you can think of this as the scatter away from the fitted line). Finally let's assume that we have ${tex`N`} points, and that the noise is drawn from a zero-mean Gaussian ${tex.block` \epsilon \sim \mathcal{N}(0, \sigma I_N)`} which implies that the data is drawn from its own Gaussian distribution ${tex.block` Y \sim \mathcal{N}(\boldsymbol{X}\boldsymbol{\theta}, \sigma I_N)`} Here ${tex`\boldsymbol{X}`} is the "augmented" data matrix (the original data vector with a column/row of ones appended). Assuming the noise is Gaussian isn't strictly necessary, but then your data distribution will no longer be a Gaussian. The likelihood is just the marginal probability ${tex.block` L = P(Y | \boldsymbol{X}, \boldsymbol{\theta}) = \prod_{i=1}^N \frac{1}{\sigma\sqrt{2\pi}}\exp{\left(\frac{-(y_i - (\theta_1 x_i + \theta_2))^2}{2\sigma^2}\right)}`} Fitting the parameters (in this case ${tex`\theta_1`} and ${tex`\theta_2`}) requires finding those which maximize the likelihood (a Maximum Likelihood Estimate). Skipping the tedious math (which is left as an exercise for the reader) you can find that the best fit vector ${tex`\boldsymbol{\theta}_{MLE}`} is given by the same result as solving an ordinary least squares solution ${tex.block` \boldsymbol{\theta}_{MLE} = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T \boldsymbol{y}`} An advantage to doing this through this method, however, is that you can also find a Maximum Likelihood Estimate of ${tex` \sigma_{MLE}`}, which can provide an estimate of your noise. Again the math is left as an exercise for the reader, but the answer is ${tex.block` \sigma_{MLE} = \sum_{i=1}^N (y_i - \boldsymbol{x_i}\boldsymbol{\theta})^2 / N`} 3c:T71a,## Finite Mixtures of Gaussian Linear Regressions ### (or, how to fit *two* lines with Gaussians) Something very exciting that we can do with a Gaussian Linear Regression fit is that you can fit a *mixture* of lines. This is very closely related to a Gaussian Mixture Model (GMM), since the latter is simply a very specific case of the much more general "Finite Mixture Model". Since our line fitting produces a likelihood, we can produce a compound (mixture) likelihood in the same way we would with any other kind of probabilistic likelihood. If we assume that the dataset is comprised of ${tex`K`} possible lines, then our likelihood would be ${tex.block`L = P(Y | \boldsymbol{X}, \boldsymbol{\Theta}) = \sum_{j=1}^K \pi_j \left(\prod_{i=1}^N \frac{1}{\sigma_j\sqrt{2\pi}}\exp{\left(\frac{-(y_i - (\theta_{1,j} x_i + \theta_{2,j}))^2}{2\sigma_j^2}\right)}\right)`} with the *mixture weights* ${tex`\pi_i`} defining the amount of each line that contributes to a given point. For this extremely rudimentary example, the Maximum Likelihood Estimate for this (in this example) will be ${tex`\pi_j = 1/K`}. Note that here we now have ${tex`2K + K = 3K`} parameters to solve for (${tex`K`} from the mixing weights, and ${tex`2K`} from the two individual parameters ${tex`\theta_{1}`}, ${tex`\theta_{2}`} for each line). Due to the summation, though, there's no closed form solution for the maximum likelihood estimation for the parameters ${tex`\boldsymbol{\Theta}`}², so this equation is typically solved using a different method. ² Note that for notational purposes ${tex`\boldsymbol{\Theta}`} refers to the collective parameters of the entire mixture model, and ${tex`\boldsymbol{\theta_k}`} will refer to the parameters of an individual model within the mixture. 3d:T6ba, Interestingly, this likelihood can be solved for in two different (but equally viable) ways using the Expectation Maximizatino algorithm. The first would be to apply the above likelihood directly to the EM algorithm. This is the method used in this notebook, with details provided in ["Finite Mixtures of Generalized Linear Regression Models", Grün & Leish (2007)](https://epub.ub.uni-muenchen.de/2087/1/tr013.pdf). In short: **E-Step:** Calculate the membership probabilities for iteration ${tex`j`} (the probability that point ${tex`i`} is drawn from line ${tex`k`}) as ${tex.block`\hat{p}_{ik}^{(j)} = \frac{\pi_k^{(j-1)} P(Y_i | \boldsymbol{X}, \boldsymbol{\Theta}_k^{(j-1)})}{\sum_{l=1}^{K} \pi_l^{(j-1)} P(Y_i | \boldsymbol{X}, \boldsymbol{\Theta}_l^{(j-1)})}`} **M-Step:** Using the membership probabilities, update the estimate of the parameters. The mixture weights are easy: ${tex.block` \pi_k^{(j)} \leftarrow \frac{1}{N} \sum_{i=1}^N \hat{p}_{ik}^{(j)}`} The functional parameters are slightly more challenging, but I will leave it as an exercise for the reader to prove that the MLE parameters are the solution to a weighted linear system such that ${tex.block` \boldsymbol{\theta_k}^{(j)} \leftarrow (\boldsymbol{X}^T\boldsymbol{W_k}^{(j)}\boldsymbol{X})^{-1}\boldsymbol{X}^T \boldsymbol{W_k}^{(j)}\boldsymbol{y}`} where the weights matrix ${tex`\boldsymbol{W_k}^{(j)}`} is a diagonal matrix whose entries are given by the ${tex`k`}-th column of ${tex`\hat{p}_{ik}^{(j)}`}. And finally the estimate for the noise is somewhat simple, also left as an exercise for the reader: ${tex.block` \sigma_{k, MLE}^{(j)} \leftarrow \sum_{i=1}^N \hat{p}_{ik} (y_i - \boldsymbol{x_i}\boldsymbol{\theta}^{(j)})^2`} 3e:T40e,{ let x_lines = [0, 10] let big_X = augment_x(tf.tensor(x_lines)) let y_fit_1 = tf.dot(big_X, two_line_fit[1]).dataSync() let y_fit_2 = tf.dot(big_X, two_line_fit[2]).dataSync() var c1 = "red" var c2 = "blue" if (two_line_fit[1].dataSync()[0] > two_line_fit[2].dataSync()[0]){ c1 = "blue" c2 = "red" } return Plot.plot({ marks: [ Plot.dot(x_1, {x: x_1, y: y_1}), Plot.dot(x_2, {x: x_2, y: y_2}), Plot.ruleY([0]), Plot.ruleY([42], {opacity:0}), Plot.ruleX([0]), Plot.line(x_lines, {x: x_lines, y: y_fit_1, stroke: c1}), Plot.line(x_lines, {x: x_lines, y: y_fit_2, stroke: c2}), Plot.areaY(x_lines, {x: x_lines, y1: y_fit_1.map(ele => ele + two_line_fit[3].dataSync()[0]), y2: y_fit_1.map(ele => ele - + two_line_fit[3].dataSync()[0]), opacity:0.2, fill:c1}), Plot.areaY(x_lines, {x: x_lines, y1: y_fit_2.map(ele => ele + two_line_fit[4].dataSync()[0]), y2: y_fit_2.map(ele => ele - + two_line_fit[4].dataSync()[0]), opacity:0.2, fill:c2}) ] }) }3f:Tb5a,// This is the cell that's going to fit the double regression two_line_history = { return tf.tidy(() => { var history = [] // Returning the history so we can plot the steps individually var x_cat = tf.concat([x_1, x_2]); var y_cat = tf.concat([y_1, y_2]); var big_X = augment_x(x_cat); // Randomly assign class probabilities (we are going to // use the true value of 2 regressions here) // We start with an initial guess of the probabilities because it's // easier to pick a reasonable starting point that trying to pick // a good initial estimate of the regressions. // Speaking of, I pick a "good" starting point" by essentially // assigning the probability that it matches the "correct" line to be // between 0.45 and 0.75 var p_k_1 = tf.randomUniform([x_1.length, 1]).mul(0.3).add(0.45); var p_k_2 = tf.randomUniform([x_2.length, 1]).mul(0.3).add(0.45); var p_k = tf.concat([p_k_1, tf.sub(1, p_k_2)], 0); p_k = tf.concat([p_k, tf.sub(1, p_k)], 1); // Finds the best fitting parameters using the p_ks as weights in // a weighted least squares fit function M_step(p_k){ let w_1 = p_k.gather([0], 1).flatten(); let w_2 = p_k.gather([1], 1).flatten(); let theta_1 = linear_fit(big_X, w_1, y_cat); let theta_2 = linear_fit(big_X, w_2, y_cat); let pi = tf.sum(p_k, 0).div(x_cat.shape[0]); // Calculating the sigma for the first line var y_pred = tf.dot(big_X, theta_1) var sum_w = tf.sum(w_1) let sigma_1 = tf.sqrt(tf.sum(tf.squaredDifference(y_pred, y_cat).mul(w_1)).div(sum_w)) // Calculating the sigma for the second line y_pred = tf.dot(big_X, theta_2) sum_w = tf.sum(w_2) let sigma_2 = tf.sqrt(tf.sum(tf.squaredDifference(y_pred, y_cat).mul(w_2)).div(sum_w)) return [pi, theta_1, theta_2, sigma_1, sigma_2] } // Finds the class probabilities using the parameters to generate two // lines, which then we calculate the p(y|X) likelihood from the mixture // weights pi function E_step(pi, theta_1, theta_2, sigma_1, sigma_2) { // Normal likelihood for the first line var y_pred = tf.dot(big_X, theta_1) let p_1 = normal(y_cat, y_pred, sigma_1).reshape([-1, 1]) // Normal likelihood for the second line y_pred = tf.dot(big_X, theta_2) let p_2 = normal(y_cat, y_pred, sigma_2).reshape([-1, 1]) // Concat and then normalize to sum to 1 along the rows let p_k = tf.concat([p_1, p_2], 1) return p_k.div(p_k.sum(1).reshape([-1, 1])) } // Run this E/M for idk 5 iterations var fit_pars = M_step(p_k) history.push(fit_pars) for (let i=0; i<3; i++){ p_k = E_step(fit_pars[0], fit_pars[1], fit_pars[2], fit_pars[3], fit_pars[4]) fit_pars = M_step(p_k) history.push(fit_pars) } return history }); }40:[] 39:{"id":"b442a54c9622cfef","slug":"finite-mixtures-of-gaussian-linear-regressions","title":"Finite Mixtures of Gaussian Linear Regressions","thumbnail":"37dbb8d107b8b4008fb954e5d65f6f24899dff31f0e7982d5c4e0f0b211300e7","access_level":"public","create_time":"$D2023-05-17T06:37:48.329Z","update_time":"$D2023-05-23T20:56:54.256Z","view_time":null,"first_public_time":"$D2023-05-17T06:37:48.329Z","has_unpublished_edits":null,"fork":false,"stars":0,"starred":false,"owner":{"type":"team","id":"7247f6820a7fa0d8","login":"dylanagreen","name":"Dylan Green","avatar_url":"https://avatars-next.observableusercontent.com/fc84910f0bd7293f70f25cdf3b35c9e385100abd98e1784446d23a0578697cb4","tier":"starter_2024","bio":"I've tried to write a chess engine for 10 years. Currently a graduate student in Physics.","following":false},"version":877,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"cdd8b4cee95ab7fe","login":"dylanagreen","name":"Dylan Green","avatar_url":"https://avatars-next.observableusercontent.com/fc84910f0bd7293f70f25cdf3b35c9e385100abd98e1784446d23a0578697cb4","tier":null,"bio":"I've tried to write a chess engine for 10 years. Currently a graduate student in Physics.","following":false}],"members":[{"type":"individual","id":"cdd8b4cee95ab7fe","login":"dylanagreen","name":"Dylan Green","avatar_url":"https://avatars-next.observableusercontent.com/fc84910f0bd7293f70f25cdf3b35c9e385100abd98e1784446d23a0578697cb4","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":1,"collection":null,"type":"notebook","body":{"titling":"auto","title":"Finite Mixtures of Gaussian Linear Regressions","cells":[{"id":603,"value":"# Finite Mixtures of Gaussian Linear Regressions\n¹Note that linear in means linear in the *parameters*, not linear in the *variables*. \n\n## Preliminaries\n### (or, how to fit a line using a Gaussian)\n When fitting a Generalized Linear Model (GLM), the central idea is that your response, ${tex`Y`}, is drawn from some distribution that is proportional to (some function of) ${tex`x`}. A GLM can use any number of exponential family distributions, but for the purposes of explanation let us use the Gaussian (or Normal) distribution, because it is one of the most common and because it helps illustrate this point.\n\n You can use this probabilistic interpretation to fit a line with a Gaussian distribution\n\n ${tex.block`P(y|x) \\propto \\exp{\\left(\\frac{-(y - g(x))^2}{2 \\sigma^2}\\right)}`}\n\n where ${tex `g(x)`} is the function you're trying to fit. \n","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":9,"value":"$3b","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":614,"value":"\nPlay with the sliders below, and see how the best fit line changes! The green represents one ${tex`\\sigma_{MLE}`} away from the line, and you can see in the table below how good of an estimate our best fit line is (for all parameters, ${tex`\\sigma_{MLE}`} included).","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":22,"value":"viewof sigma = Inputs.range([0, 3], {label: tex`\\sigma=`, step: 1e-3, value: 0.5})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":29,"value":"viewof slope = Inputs.range([0, 3], {label: tex`m=`, step: 1e-3})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":38,"value":"viewof intercept = Inputs.range([0, 10], {label: tex`b=`, step: 1e-3})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":75,"value":"viewof n_points = Inputs.range([5, 30], {label: tex`N=`, step: 1, value: 10})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":42,"value":"{\n let y_pred = await fit_par[0].data()\n return Plot.plot({\n marks: [\n Plot.dot(x, {x: x, y: y}),\n Plot.ruleY([0]),\n Plot.ruleY([45], {opacity:0}),\n Plot.ruleX([0]),\n Plot.line(x, {x: x, y: y_pred, stroke: \"green\"}),\n Plot.areaY(x, {x: x, y1: y_pred.map(ele => ele + sigma_pred), y2: y_pred.map(ele => ele - sigma_pred), opacity:0.2, fill:\"green\"})\n ]\n })\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":453,"value":" | ${tex`b`} | ${tex`m`} | ${tex`\\sigma`}\n-------- | -------- | -------- | --------\nTrue | ${intercept} | ${slope} | ${sigma} | \nPredicted | ${theta_fit[0].toFixed(3)} | ${theta_fit[1].toFixed(3)}| ${sigma_pred.toFixed(3)}","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":301,"value":"$3c","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":647,"value":"$3d","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":653,"value":"\nAn alternative method, which I won't go into detail here but which is detailed in [\"Many regression algorithms, one unified model: A review\", Stulp & Sigaud (2015)](https://www.sciencedirect.com/science/article/pii/S0893608015001185) , is to instead fit a mixture of gaussians (a GMM) to the joint ${tex`P(Y, \\boldsymbol{X})`} density, and then convert that to a marginal density using some simple integrational formulas. See page 68 of the aforementioned paper for more details. \n","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":656,"value":"Given all of that, we can generate a dataset drawn from two different lines, and then simultaneously fit 2 lines to the dataset to try and recover *both* of the constituent lines *simultaneously*. There are some sliders below to adjust the constituent dataset, and an additional plot to step through each of the iterations of the EM algorithm. **Note: I always terminate the EM algorithm after 3 iterations, and I also pick a semi clever starting point (because I know which point belongs to which line), since this is just for demonstration purposes and not real scientific value.**","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":713,"value":"**Line 1**","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":660,"value":"viewof sigma_1 = Inputs.range([0, 3], {label: tex`\\sigma_1=`, step: 1e-3, value: 1.25})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":661,"value":"viewof slope_1 = Inputs.range([0, 3], {label: tex`m_1=`, step: 1e-1})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":662,"value":"viewof intercept_1 = Inputs.range([0, 10], {label: tex`b_1=`, step: 1e-3})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":715,"value":"**Line 2**","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":701,"value":"viewof sigma_2 = Inputs.range([0, 3], {label: tex`\\sigma_2=`, step: 1e-3, value: 1.25})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":702,"value":"viewof slope_2 = Inputs.range([-3, 0], {label: tex`m_2=`, step: 1e-1})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":703,"value":"viewof intercept_2 = Inputs.range([15, 25], {label: tex`b_2=`, step: 1e-2})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":165,"value":"$3e","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":774,"value":"viewof two_line_idx = Inputs.range([1, two_line_history.length], {label: md`**Iteration Number:**`, step: 1, value: two_line_history.length})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":745,"value":"| ${tex`m_1`} | ${tex`b_1`} | ${tex`\\sigma_1`} | ${tex`m_2`} | ${tex`b_2`} | ${tex`\\sigma_2`}\n-------- | -------- | -------- | --------\nTrue | ${slope_1} | ${intercept_1} | ${sigma_1} | ${slope_2} | ${intercept_2} | ${sigma_2} | \nPredicted | ${fits_ordered[0][1].toFixed(3)} | ${fits_ordered[0][0].toFixed(3)}| ${fits_ordered[2][0].toFixed(3)}| ${fits_ordered[1][1].toFixed(3)} | ${fits_ordered[1][0].toFixed(3)}| ${fits_ordered[3][0].toFixed(3)}","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":847,"value":"---","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":469,"value":"# Behind the Scenes\n## Generic\n","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":34,"value":"x = {\n // We randomly slightly off set the x values so that \n // they're not uniformly distributed\n let xoffset = d3.randomNormal.source(d3.randomLcg(0.100821))(0, 0.1)\n return linspace(0.1, 10, n_points).map(ele => ele + xoffset(1))\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":358,"value":"// x should be a tensor, mu can be either\nfunction normal(x, mu, sigma) {\n let exponential = tf.exp(tf.mul(-0.5, tf.squaredDifference(x, mu).div(sigma)))\n let pre = tf.pow(tf.mul(sigma, Math.sqrt(2 * Math.PI)), -1)\n return tf.mul(pre, exponential)\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":116,"value":"function linear_fit(X, w, y) {\n let A_tf = tf.dot(tf.transpose(X), tf.mul(w.reshape([-1, 1]), X));\n let qr = tf.linalg.qr(A_tf); // This might end up being slow for large matrices...\n\n let y_tf = tf.dot(tf.transpose(X), tf.mul(w, y));\n return backsolve(qr[1], tf.dot(tf.transpose(qr[0]), y_tf));\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":141,"value":"// In lieu of an actual linspace option this mapping will have to do.\nfunction linspace(start, stop, num) {\n return d3.range(num).map(ele => start + (ele * (stop - start)) / num)\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":129,"value":"function augment_x(x) {\n // Tensorflow functions are almost identical to python versions\n // Reshape x to have an additional dimension, the dimension we concatenate along\n // Note that this is now in the \"correct\" form, i.e. 1s on column, x on column\n return tf.concat([tf.ones([x.shape[0], 1]), x.reshape([-1, 1])], 1);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":92,"value":"import {backsolve} from \"@dylanagreen/tf-linalg\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":94,"value":"tf = require('@tensorflow/tfjs')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":512,"value":"## 1-Line","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":36,"value":"y = x.map(ele => slope * ele + intercept + noise(1))","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":16,"value":"noise = d3.randomNormal.source(d3.randomLcg(0.91701))(0, sigma)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":325,"value":"theta_fit = fit_par[1].data()","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":98,"value":"sigma_pred = {\n var y_pred = fit_par[0].dataSync();\n return Math.sqrt(d3.sum(y_pred.map((ele, i) => (ele - y[i])**2)) / n_points);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":90,"value":"fit_par = {\n // tf.tidy deallocates any variables or space we use that we don't need again when we're done\n return tf.tidy(() => {\n let X_tf = augment_x(tf.tensor(x))\n let theta_tf = linear_fit(X_tf, tf.onesLike(x), y);\n return [tf.dot(X_tf, theta_tf), theta_tf];\n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":471,"value":"## 2-Line","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":749,"value":"fits_ordered = {\n let fit_1 = two_line_fit[1].dataSync()\n let fit_2 = two_line_fit[2].dataSync()\n\n let s_1 = two_line_fit[3].dataSync()\n let s_2 = two_line_fit[4].dataSync()\n\n if (fit_1[1] > fit_2[1]) {\n return [fit_1, fit_2, s_1, s_2]\n }\n else {\n return [fit_2, fit_1, s_1, s_2]\n }\n\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":139,"value":"x_1 = {\n let xoffset = d3.randomNormal.source(d3.randomLcg(0.100821))(0, 0.25)\n return linspace(0.1, 10, 30).map(ele => ele + xoffset(4))\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":155,"value":"x_2 = {\n let xoffset = d3.randomNormal.source(d3.randomLcg(0.062796))(0, 0.25)\n return linspace(0.1, 10, 30).map(ele => ele + xoffset(1))\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":446,"value":"noise_2 = d3.randomNormal.source(d3.randomLcg(0.91701))(0, sigma_2)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":694,"value":"noise_1 = d3.randomNormal.source(d3.randomLcg(0.91701))(0, sigma_1)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":177,"value":"y_1 = x_1.map(ele => slope_1 * ele + intercept_1 + noise_1(1))","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":180,"value":"y_2 = x_2.map(ele => slope_2 * ele + intercept_2 + noise_2(2))","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":768,"value":"two_line_fit = two_line_history[two_line_idx-1]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":195,"value":"$3f","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":823,"value":"","pinned":true,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"}],"files":"$W40","stdlib":"1"}} 3a:"$39:body"