Кубическая регрессия (линия наилучшего соответствия) в JavaScript

Вот настоящий рабочий код для решения этой кубической задачи с использованием uncmin неограниченного минимизатора библиотеки numeric.js в качестве задача наименьших квадратов (jsbin здесь):

var data_x = [2,5,7,12,20,32,50];
var data_y = [5,10,15,20,25,30,35];

var cubic = function(params,x) {
  return params[0] * x*x*x +
    params[1] * x*x +
    params[2] * x +
    params[3];
};

var objective = function(params) {
  var total = 0.0;
  for(var i=0; i < data_x.length; ++i) {
    var resultThisDatum = cubic(params, data_x[i]);
    var delta = resultThisDatum - data_y[i];
    total += (delta*delta);
  }
  return total;
};

var initial = [1,1,1,1];
var minimiser = numeric.uncmin(objective,initial);

console.log("initial:");
for(var j=0; j<initial.length; ++j) {
  console.log(initial[j]);  
}

console.log("minimiser:");
for(var j=0; j<minimiser.solution.length; ++j) {
  console.log(minimiser.solution[j]);
}

Я получаю результаты:

 0.0005750849851827991
-0.05886106462847641
 2.1839575020602164
 1.1276055079334206

Поясню: у нас есть функция 'cubic', которая оценивает общую кубическую функцию для набора параметров params и значения x. Эта функция обернута для создания целевой функции, которая принимает набор параметров и пропускает каждое значение x из нашего набора данных через целевую функцию и вычисляет сумму квадратов. Эта функция передается в uncmin из numeric.js с набором начальных значений; uncmin выполняет тяжелую работу и возвращает объект, свойство solution которого содержит оптимизированный набор параметров.

Чтобы сделать это без глобальных переменных (капризно!), вы можете иметь фабрику целевых функций следующим образом:

var makeObjective = function(targetFunc,xlist,ylist) {
  var objective = function(params) {
    var total = 0.0;
    for(var i=0; i < xlist.length; ++i) {
      var resultThisDatum = targetFunc(params, xlist[i]);
      var delta = resultThisDatum - ylist[i];
      total += (delta*delta);
    }
    return total;
  };
  return objective;
};

Которые вы можете использовать для создания целевых функций:

var objective = makeObjective(cubic, data_x, data_y); // then carry on as before

Знание того, как это сделать на практике, очень помогло бы многим людям, поэтому я рад, что это всплыло.

Изменить: разъяснение по cubic

var cubic = function(params,x) {
  return params[0] * x*x*x +
    params[1] * x*x +
    params[2] * x +
    params[3];
};

Cubic определяется как функция, которая принимает массив параметров params и значение x. Учитывая params, мы можем определить функцию f(x). Для куба это f(x) = a x^3 + b x^2 + c x + d, поэтому есть 4 параметра (от [0] до [3]), и с учетом этих 4 значений параметров у нас есть одна функция f(x) с 1 входом x.

Код структурирован так, чтобы вы могли заменить cubic другой функцией той же структуры; это может быть linear с двумя параметрами:

var linear = function(params, x) {
    return params[0]*x + params[1];
};

Остальная часть кода будет смотреть на длину params, чтобы узнать, сколько параметров нужно изменить.

Обратите внимание, что весь этот фрагмент кода пытается найти набор значений параметров, которые создают кривую, которая наилучшим образом соответствует всем данным; если вы хотите найти соответствие для последних 4 точек некоторых данных, вы должны передать только эти значения в data_x и data_y.

Я бы сформулировал это как задачу наименьших квадратов. Пусть M будет матрицей n×4, сформированной следующим образом:

x_1^3  x_1^2  x_1  1
x_2^3  x_2^2  x_2  1
  ⋮       ⋮      ⋮
x_n^3  x_n^2  x_n  1

Затем вычислите матрицу 4 × 4 A=M^T⋅M и вектор-столбец 4 × 1 b=M^T⋅y и решить линейную систему уравнений Aξ=б. Результирующий вектор ξ будет содержать ваши коэффициенты от a до d.

Приведенное выше описание позволяет легко понять, что происходит, математически. Однако для реализации, особенно для очень больших n, описанный выше подход может оказаться неприемлемым. В таких случаях вы можете построить A и b напрямую, без явного построения M. Например, A_1,2=sum(x_i^3 * x_i^2 for all i). Таким образом, вы можете перебрать все i и добавить соответствующие значения в соответствующие элементы матрицы и вектора.