2013-08-31 1 views
12

У меня есть дискретная регулярная сетка из a,b точек и их соответствующих значений c, и я интерполирую ее дальше, чтобы получить плавную кривую. Теперь из данных интерполяции я также хочу создать полиномиальное уравнение для подгонки кривой. Как разместить 3D-график в полиноме?3D curvefitting

Я пытаюсь сделать это в MATLAB. Я использовал панель инструментов для поверхностного монтажа в MATLAB (r2010a) для трехмерных данных. Но как найти формулу, которая наилучшим образом соответствует набору данных в MATLAB/MAPLE или любом другом программном обеспечении. Любой совет? Также наиболее полезными были бы некоторые реальные примеры кода для просмотра, файлы PDF, в Интернете и т. Д.

Это всего лишь небольшая часть моих данных.

a = [ 0.001 .. 0.011]; 

b = [1, .. 10]; 

c = [ -.304860225, .. .379710865]; 

Заранее благодарен.

ответ

15

Чтобы соответствовать кривой на множестве точек, мы можем использовать регрессию ordinary least-squares. В MathWorks описывается процесс solution page.

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

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

Как @BasSwinckels показал, построив нужный design matrix, вы можете использовать mldivide или pinv к solve the overdetermined system выражается в 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)) 

1st_order_polynomial


Как упоминались выше, мы можем получить полиномиальную подгонку высшего порядка, добавляя больше терминов к независимой матрице переменных (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 в Statistics Toolbox, который помогает в создании дизайна матрицы для нескольких модельных заказов:

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

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

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

2nd_order_polynomial

+0

Как я могу выполнить такую ​​же операцию в python ..!? Бит-помощь будет оценена ... @Amro – diffracteD

+2

@diffracteD: Я перевел код на Python: https://gist.github.com/amroamroamro/1db8d69b4b65e8bc66a6 – Amro

+0

Спасибо за вашу помощь ... я обязательно попробую. ! – diffracteD

3

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

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') 

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

+1

+1 - третий порядок, соответствующий «M = [единицы (размер (x)), x, y, x.^2, x. * Y, y.^2, x.^3, x .^2. * y, x. * Y.^2, y.^3] 'делает довольно хорошую работу ... –

+0

Извините за поздний ответ! Спасибо, что ответили Бас Суинкельс и старались. Хорошее решение. Но есть ли способ найти наилучшую формулу формулы/уравнения через код или с помощью любого инструмента в MATLAB? Как найти формулу, которая лучше всего соответствует указанному набору данных? – Syeda

+0

Я тоже получаю это предупреждение. Предупреждение: Ранг дефицит, ранг = 5, tol = 9.9961e-013. Не могли бы вы помочь мне понять, что это значит? Спасибо. – Syeda

2

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

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

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