3D кривая

Чтобы подогнать кривую к набору точек, мы можем использовать обычный метод наименьших квадратов. Существует страница решения от MathWorks, описывающая процесс.

В качестве примера начнем с некоторых случайных данных:

% some 3d points
data = mvnrnd([0 0 0], [1 -0.5 0.8; -0.5 1.1 0; 0.8 0 1], 50);

Как показал @BasSwinckels, создав нужный матрица проектирования, вы можете использовать mldivide или pinv на решить переопределенную систему, выраженную как Ax=b:

% best-fit plane
C = [data(:,1) data(:,2) ones(size(data,1),1)] \ data(:,3);    % coefficients

% evaluate it on a regular grid covering the domain of the data
[xx,yy] = meshgrid(-3:.5:3, -3:.5:3);
zz = C(1)*xx + C(2)*yy + C(3);

% or expressed using matrix/vector product
%zz = reshape([xx(:) yy(:) ones(numel(xx),1)] * C, size(xx));

Далее визуализируем результат:

% plot points and surface
figure('Renderer','opengl')
line(data(:,1), data(:,2), data(:,3), 'LineStyle','none', ...
    'Marker','.', 'MarkerSize',25, 'Color','r')
surface(xx, yy, zz, ...
    'FaceColor','interp', 'EdgeColor','b', 'FaceAlpha',0.2)
grid on; axis tight equal;
view(9,9);
xlabel x; ylabel y; zlabel z;
colormap(cool(64))

Как уже упоминалось, мы можем получить полиномиальную аппроксимацию более высокого порядка, добавив больше членов в матрицу независимых переменных (A в Ax=b).

Скажем, мы хотим подобрать квадратичную модель с постоянными, линейными, взаимодействующими и квадратичными членами (1, x, y, xy, x^2, y^2). Мы можем сделать это вручную:

% best-fit quadratic curve
C = [ones(50,1) data(:,1:2) prod(data(:,1:2),2) data(:,1:2).^2] \ data(:,3);
zz = [ones(numel(xx),1) xx(:) yy(:) xx(:).*yy(:) xx(:).^2 yy(:).^2] * C;
zz = reshape(zz, size(xx));

В наборе инструментов статистики также есть вспомогательная функция x2fx, помогающая в построении схемы. матрица для пары типовых заказов:

C = x2fx(data(:,1:2), 'quadratic') \ data(:,3);
zz = x2fx([xx(:) yy(:)], 'quadratic') * C;
zz = reshape(zz, size(xx));

Наконец, есть отличная функция polyfitn на файловом обмене Джона Д'Эррико, которая позволяет указать все виды полиномиальных порядков и условий:

model = polyfitn(data(:,1:2), data(:,3), 2);
zz = polyvaln(model, [xx(:) yy(:)]);
zz = reshape(zz, size(xx));

На файловом обмене могут быть некоторые лучшие функции, но один из способов сделать это вручную:

x = a(:); %make column vectors
y = b(:);
z = c(:);

%first order fit
M = [ones(size(x)), x, y];
k1 = M\z; 
%least square solution of z = M * k1, so z = k1(1) + k1(2) * x + k1(3) * y

Точно так же вы можете выполнить подгонку второго порядка:

%second order fit
M = [ones(size(x)), x, y, x.^2, x.*y, y.^2];
k2 = M\z;

у которого, кажется, есть числовые проблемы для ограниченного набора данных, который вы предоставили. Введите help mldivide для более подробной информации.

Чтобы сделать интерполяцию по некоторой регулярной сетке, вы можете сделать так:

ngrid = 20;
[A,B] = meshgrid(linspace(min(a), max(a), ngrid), ...
                 linspace(min(b), max(b), ngrid));
M = [ones(numel(A),1), A(:), B(:), A(:).^2, A(:).*B(:), B(:).^2];
C2_fit = reshape(M * k2, size(A)); % = k2(1) + k2(2)*A + k2(3)*B + k2(4)*A.^2 + ...

%plot to compare fit with original data
surfl(A,B,C2_fit);shading flat;colormap gray
hold on
plot3(a,b,c, '.r')

Подгонку 3-го порядка можно выполнить с помощью формулы, приведенной TryHard ниже, но формулы быстро становятся утомительными, когда порядок увеличивается. Лучше напишите функцию, которая может построить M по данным x, y и order, если вам нужно сделать это более одного раза.

Это больше похоже на философский вопрос, чем на конкретную реализацию, особенно на бит — «как найти формулу, которая наилучшим образом соответствует набору данных?» По моему опыту, это выбор, который вы должны сделать в зависимости от того, чего вы пытаетесь достичь.

Что для вас определяет «лучший»? Для задачи подбора данных вы можете продолжать добавлять все больше и больше полиномиальных коэффициентов и улучшать значение R ^ 2... но в конечном итоге данные будут "переоцениваться". Недостатком полиномов высокого порядка является поведение за пределами выборочных данных, которые вы использовали для соответствия вашей поверхности отклика - они могут быстро уйти в какое-то дикое направление, которое может не подходить для того, что вы пытаетесь смоделировать. .

Есть ли у вас представление о физическом поведении системы/данных, которые вы подгоняете? Это можно использовать в качестве основы для того, какой набор уравнений использовать для создания математической модели. Я бы порекомендовал использовать самую экономичную (простую) модель, с которой вы можете сойти с рук.