31:T6b3,


    x = [32,37,42,47,52,57,62,67,72,77,82,87,92]   

    y = [749,1525,1947,2201,2380,2537,2671,2758,2803,2943,3007,2979,2992]   
 
    model = function(x,p){return x.map(function(x_i){return ( 

       p[0]*Math.exp(-0.5*(Math.pow((x_i-p[1])/p[2],2))) + p[3]    
 
    )})}  

    p_init = [-60000,-90,40,3100];  

    Opt = { 

     maxIter: 500,            // Optional term

     bounds: [                // No bounds: [-Infinity, Infinity],

          [-70000, -50000],  // Bounds for p[0] 

          [-90, 100],           // Bounds for p[1] 

          [1, 40],                    // Bounds for p[2] 

          [3000, 3100]         // Bounds for p[3] 

      ], 

    } 

   }; 

    fitting = fminsearch(model,p_init,x,y,Opt)

32:T964,// Modified fminsearch function with bounds: Please note that this implementation is not optimized. // Adapted from: https://github.com/jonasalmeida/fminsearch // See associated notebook: https://observablehq.com/@christophe-yamahata/nonlinear-regression-with-bounds-in-javascript fminsearch = function(fun, Parm0, x, y, Opt) { if (!Opt) { Opt = {} }; if (!Opt.maxIter) { Opt.maxIter = 1000 }; if (!Opt.step) { // initial step is 1/100 of initial value (remember not to use zero in Parm0) Opt.step = Parm0.map(function(p) { return p / 100 }); Opt.step = Opt.step.map(function(si) { if (si == 0) { return 1 } else { return si } }); // convert null steps into 1's }; if (typeof (Opt.display) == 'undefined') { Opt.display = true }; if (!Opt.objFun) { Opt.objFun = function(y, yp) { return y.map(function(yi, i) { return Math.pow((yi - yp[i]), 2) }).reduce(function(a, b) { return a + b }) } } //SSD if (!Opt.bounds) { Opt.bounds = Parm0.map(function() { return [-Infinity, Infinity] }) }; // Default no bounds var cloneVector = function(V) { return V.map(function(v) { return v }) }; var ya, y0, yb, fP0, fP1; var P0 = cloneVector(Parm0), P1 = cloneVector(Parm0); var n = P0.length; // var step = Opt.step; let step = cloneVector(Opt.step); // Ensure step is reinitialized var funParm = function(P) { return Opt.objFun(y, fun(x, P)) } //function (of Parameters) to minimize // Apply bounds var applyBounds = function(P, bounds) { return P.map(function(p, i) { return Math.max(bounds[i][0], Math.min(bounds[i][1], p)); }); } // Silly multi-univariate screening for (var i = 0; i < Opt.maxIter; i++) { for (var j = 0; j < n; j++) { // take a step for each parameter P1 = cloneVector(P0); P1[j] += step[j]; P1 = applyBounds(P1, Opt.bounds); // Apply bounds after updating the parameter if (funParm(P1) < funParm(P0)) { // if parm value going in the right direction step[j] = 1.2 * step[j]; // then go a little faster P0 = cloneVector(P1); } else { step[j] = -(0.5 * step[j]); // otherwise reverse and go slower } } if (Opt.display) { if (i > (Opt.maxIter - 10)) { console.log(i + 1, funParm(P0), P0) } } } return P0 }33:T6ea,python_script = ` """ This Python script performs nonlinear regression with bounded parameters to fit a model to data points. The script uses least squares optimization from the SciPy library. For further details and a similar implementation in JavaScript, visit: https://observablehq.com/@christophe-yamahata/nonlinear-regression-with-bounds-in-javascript Author: C. Yamahata Date: May 2024 """ import numpy as np from scipy.optimize import least_squares def model(x, params): """ Gaussian model with an offset: A * exp(-0.5 * ((x - μ) / σ)^2) + B """ amplitude, mean, stddev, offset = params return amplitude * np.exp(-0.5 * np.power((x - mean) / stddev, 2)) + offset def objective_function(params, x, y): """ Objective function for least squares optimization. Returns the residuals between the model and the observed data. """ return model(x, params) - y # Given data points x_data = np.array([32, 37, 42, 47, 52, 57, 62, 67, 72, 77, 82, 87, 92]) y_data = np.array([749, 1525, 1947, 2201, 2380, 2537, 2671, 2758, 2803, 2943, 3007, 2979, 2992]) # Initial guess for the parameters initial_params = np.array([-60000, -90, 40, 3100]) # Define lower and upper bounds for each parameter lower_bounds = [-70000, -90, 1, 3000] upper_bounds = [-50000, 100, 40, 3100] bounds = (lower_bounds, upper_bounds) # Optimize parameters using least squares with bounds result = least_squares(objective_function, initial_params, bounds=bounds, args=(x_data, y_data)) # Print best-fit parameters param_names = ["Amplitude, A (p_0)", "Mean, μ (p_1)", "stddev, σ (p_2)", "offset, B (p_3)"] print("\\nBest-fit parameters:") for name, value in zip(param_names, result.x): print(f"{name:<16} : {value:.6f}") `34:[{"id":"88924e508f20d73a94646045d223a6dbc244aa64fef2498100f2a37bc97bb4f8ed4c03b438888a18b70fbacf46756028c3cc0989365603a88fc02a18e6aee53f","href":"https://static.observableusercontent.com/files/88924e508f20d73a94646045d223a6dbc244aa64fef2498100f2a37bc97bb4f8ed4c03b438888a18b70fbacf46756028c3cc0989365603a88fc02a18e6aee53f","name":"Python_output.png","type":"image/png","lastModified":1716805161290.474,"size":6075}] 2f:{"title":"Nonlinear regression with bounds in JavaScript","cells":[{"id":0,"value":"# Nonlinear regression with bounds in JavaScript","pinned":false,"hidden":false,"mode":"md"},{"id":922,"value":"

This notebook is a fork of a previous notebook on nonlinear regression in JavaScript.

\n\n

\nThe aim here is to implement a similar nonlinear regression algorithm while imposing bounding limits on the parameters.\n

","pinned":false,"hidden":false,"mode":"html"},{"id":665,"value":"
\n

Nonlinear regression with or without bounds

The modified fminsearch.js JavaScript function used in this notebook is based on the function originally implemented by @jonasalmeida. The original function can be found on GitHub.\n

","pinned":false,"hidden":false,"mode":"html"},{"id":845,"value":"
\n

Example Usage with Bounds

","pinned":false,"hidden":false,"mode":"html"},{"id":833,"value":"$31","pinned":false,"hidden":false,"mode":"html"},{"id":792,"value":"

\n First, the x and y data are stored in two arrays:\n

","pinned":false,"hidden":false,"mode":"html"},{"id":12,"value":"x = [32,37,42,47,52,57,62,67,72,77,82,87,92]","pinned":true,"hidden":false,"mode":"ojs"},{"id":14,"value":"y = [749,1525,1947,2201,2380,2537,2671,2758,2803,2943,3007,2979,2992]","pinned":true,"hidden":false,"mode":"ojs"},{"id":1368,"value":"

Next, we define the model. Here, we assume the function ${tex`y = f(x)`}:

","pinned":false,"hidden":false,"mode":"html"},{"id":1360,"value":"tex.block`\n f(x) = a_0 e^{-\\frac{1}{2} (\\frac{x-a_1}{a_2})^2 } + a_3\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1373,"value":"

where ${tex`a_0`}, ${tex`a_1`}, ${tex`a_2`}, ${tex`a_3`} are the fitting parameters to be determined.

","pinned":false,"hidden":false,"mode":"html"},{"id":886,"value":"

From the analytical model given above, we write:

","pinned":false,"hidden":false,"mode":"html"},{"id":18,"value":"model = function(x,p){return x.map(function(x_i){return (\n p[0]*Math.exp(-0.5*(Math.pow((x_i-p[1])/p[2],2))) + p[3]\n)})}","pinned":true,"hidden":false,"mode":"ojs"},{"id":786,"value":"

Then, the initial guess parameters are loaded into the array p_init:

","pinned":false,"hidden":false,"mode":"html"},{"id":28,"value":"p_init = [-60000,-90,40,3100];","pinned":true,"hidden":false,"mode":"ojs"},{"id":1167,"value":"

Next, we define the bounds for these parameters:

","pinned":false,"hidden":false,"mode":"html"},{"id":1154,"value":"// Define bounds for the parameters\nOpt = {\n return {\n maxIter: 500, // Optional\n bounds: [\n [-70000, -50000], // Bounds for p[0] => [-Infinity, Infinity], // No bounds\n [-90, 100], // Bounds for p[1]\n [1, 40], // Bounds for p[2]\n [3000, 3100] // Bounds for p[3]\n ],\n }\n};","pinned":true,"hidden":false,"mode":"ojs"},{"id":790,"value":"

\n Finally, the nonlinear regression with bounding limits is completed using fminsearch as follows:\n fitting_with_bounds = fminsearch(model,p_init,x,y,Opt)\n

","pinned":false,"hidden":false,"mode":"html"},{"id":1525,"value":"The resulting parameters with bounds are:","pinned":false,"hidden":false,"mode":"html"},{"id":22,"value":"fitting_with_bounds = fminsearch(model,p_init,x,y,Opt)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1304,"value":"

\n This can be compared with the solution obtained when the function is called without the optional bounds:\n fitting_without_bounds = fminsearch(model,p_init,x,y)\n

","pinned":false,"hidden":false,"mode":"html"},{"id":1532,"value":"Resulting in:","pinned":false,"hidden":false,"mode":"html"},{"id":1300,"value":"fitting_without_bounds = fminsearch(model,p_init,x,y)","pinned":false,"hidden":false,"mode":"ojs"},{"id":901,"value":"
\n

Graphic illustration

\n In the following graph, you can compare the final results \n
(__ with bounds / _ _ without bounds) to the initial estimate.\n

","pinned":false,"hidden":false,"mode":"html"},{"id":61,"value":"Plot.plot({\n inset: 0,\n grid: true,\n height: 450,\n width: 800, \n marginLeft: 60,\n marginTop: 40, \n marginBottom: 50,\n x: {\n label: \"X → \",\n tickPadding: 10, \n labelOffset: 40,\n\n },\n y: {\n label: \"↑ Y\",\n tickPadding: 10,\n labelOffset: 20, \n },\n style: {\n fontSize: \"14px\",\n \"font-weight\": \"bold\",\n },\n marks: [\n Plot.ruleX([x_min,x_max]),\n Plot.ruleY([y_min,y_max]),\n \n Plot.line(initial_curve, {x: \"x\", y: \"y\", curve: \"catmull-rom\", stroke: color_initial, strokeOpacity: 0.5}), \n \n Plot.line(fitting_with_bounds_curve, {x: \"x\", y: \"y\", curve: \"catmull-rom\", stroke: color_stroke1}),\n\n Plot.line(fitting_without_bounds_curve, {x: \"x\", y: \"y\", curve: \"catmull-rom\", stroke: color_stroke2, strokeDasharray: \"6,10\"}), \n \n Plot.dot(data_curve, {x: \"x\", y: \"y\", curve: \"catmull-rom\", fill: \"yellow\", stroke: color_dot, r: 3.6, strokeOpacity: 0.7}), \n\n\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1493,"value":"

Nonlinear regression with bounds in Python

\nThe following script provides a similar implementation in Python. The script uses least squares optimization from the SciPy library.\n

","pinned":false,"hidden":false,"mode":"html"},{"id":1495,"value":"Inputs.button(\"Copy Python script to clipboard\", {value: null, reduce: () => navigator.clipboard.writeText(python_script)})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1542,"value":"
\n

\n
Results Generated by the Python Script.\n

\n
","pinned":false,"hidden":false,"mode":"html"},{"id":1499,"value":"html`

\n${python_script}\n

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":510,"value":"---\n## Nonlinear regression with bounds: JavaScript function","pinned":false,"hidden":false,"mode":"md"},{"id":5,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":25,"value":"---\n## Appendix","pinned":false,"hidden":false,"mode":"md"},{"id":1185,"value":"data_curve = {\n var data = [];\n for(var i = 0; i < x.length; i++) \n data.push({\"x\":x[i],\"y\":y[i]});\n return data;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":157,"value":"initial_curve = {\n var initial_curve = [];\n var x_ = [];\n var y_ = [];\n const n = 500;\n const dx = (x_max-x_min)/n;\n for(var i = 0; i < n; i++) \n {\n x_[i] = x_min + dx*i;\n y_[i] = p_init[0]*Math.exp(-0.5*(Math.pow((x_[i]-p_init[1])/p_init[2],2))) + p_init[3]\n\n initial_curve.push({\"x\":x_[i],\"y\":y_[i]});\n }\n\n return initial_curve;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":134,"value":"fitting_with_bounds_curve = {\n var fitting_curve = [];\n var x_ = [];\n var y_ = [];\n const n = 500;\n const dx = (x_max-x_min)/n;\n for(var i = 0; i < n; i++) \n {\n x_[i] = x_min + dx*i;\n y_[i] = fitting_with_bounds[0]*Math.exp(-0.5*(Math.pow((x_[i]-fitting_with_bounds[1])/fitting_with_bounds[2],2))) + fitting_with_bounds[3];\n\n fitting_curve.push({\"x\":x_[i],\"y\":y_[i]});\n }\n\n return fitting_curve;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1314,"value":"fitting_without_bounds_curve = {\n var fitting_curve = [];\n var x_ = [];\n var y_ = [];\n const n = 500;\n const dx = (x_max-x_min)/n;\n for(var i = 0; i < n; i++) \n {\n x_[i] = x_min + dx*i;\n y_[i] = fitting_without_bounds[0]*Math.exp(-0.5*(Math.pow((x_[i]-fitting_without_bounds[1])/fitting_without_bounds[2],2))) + fitting_without_bounds[3];\n\n fitting_curve.push({\"x\":x_[i],\"y\":y_[i]});\n }\n\n return fitting_curve;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":122,"value":"color_dot = \"darkblue\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":454,"value":"color_stroke1 = \"orange\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1312,"value":"color_stroke2 = \"purple\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":465,"value":"color_initial = \"blue\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":287,"value":"x_min = 30;","pinned":false,"hidden":false,"mode":"ojs"},{"id":291,"value":"x_max = 100;","pinned":false,"hidden":false,"mode":"ojs"},{"id":339,"value":"y_min = 0;","pinned":false,"hidden":false,"mode":"ojs"},{"id":341,"value":"y_max = 3500;","pinned":false,"hidden":false,"mode":"ojs"},{"id":569,"value":"// https://observablehq.com/@jashkenas/markdown-style-footnotes\nreference = (referencingCellName, number, text) => {\n return html`

↫ ${text}

`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1491,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}