Документация

Предсказание населения США

Этот пример показывает, что экстраполирование данных с помощью многочленов даже скромного градуса опасно и ненадежно.

Этот пример является более старым, чем MATLAB®. Это запустилось как упражнение в Компьютерных Методах для Математических вычислений, Форсайтом, Малкольмом и Молером, опубликованным Prentice Hall в 1 977.

Теперь, MATLAB делает намного легче отличаться параметры и видеть результаты, но базовые математические принципы неизменны.

Создайте и постройте график двух векторов с данными переписи США от 1 900 до 2000.

% Time interval
t = (1900:10:2000)';

% Population
p = [75.995 91.972 105.711 123.203 131.669 ...
   150.697 179.323 203.212 226.505 249.633 281.422]';

% Plot
plot(t,p,'bo');
axis([1900 2020 0 400]);
title('Population of the U.S. 1900-2000');
ylabel('Millions');

Каково ваше предположение для населения в году 2010?

p

p = 11×1

   75.9950
   91.9720
  105.7110
  123.2030
  131.6690
  150.6970
  179.3230
  203.2120
  226.5050
  249.6330
      ⋮

Соответствуйте данным многочленом в t и используйте его, чтобы экстраполировать население в t = 2010. Получите коэффициенты в многочлене путем решения линейной системы уравнений, включающих 11 11 матрицу Вандермонда, с элементами как степени масштабированного времени, A(i,j) = s(i)^(n-j).

n = length(t);
s = (t-1950)/50;
A = zeros(n);
A(:,end) = 1;
for j = n-1:-1:1
   A(:,j) = s .* A(:,j+1);
end

Получите коэффициенты c для многочлена градуса d, который соответствует данным p путем решения линейной системы уравнений, включающих последние столбцы d+1 матрицы Вандермонда:

A:, n-d:n) *c ~ = p

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

c = A(:,n-3:n)\p

c = 4×1

    1.2629
   23.7261
  100.3659
  155.9043

Теперь оцените многочлен в каждый год от 1 900 до 2010 и постройте график результатов.

v = (1900:2020)';
x = (v-1950)/50;
w = (2010-1950)/50;
y = polyval(c,x);
z = polyval(c,w);

hold on
plot(v,y,'k-');
plot(2010,z,'ks');
text(2010,z+15,num2str(z));
hold off

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

c = A(:,n-4:n)\p;
y = polyval(c,x);
z = polyval(c,w);

hold on
plot(v,y,'k-');
plot(2010,z,'ks');
text(2010,z-15,num2str(z));
hold off

Когда градус увеличивается, экстраполяция становится еще более ошибочной.

cla
plot(t,p,'bo')
hold on
axis([1900 2020 0 400])
colors = hsv(8);
labels = {'data'};
for d = 1:8
   [Q,R] = qr(A(:,n-d:n));
   R = R(1:d+1,:);
   Q = Q(:,1:d+1);
   c = R\(Q'*p);    % Same as c = A(:,n-d:n)\p;
   y = polyval(c,x);
   z = polyval(c,11);
   plot(v,y,'color',colors(d,:));
   labels{end+1} = ['degree = ' int2str(d)];
end
legend(labels, 'Location', 'NorthWest')
hold off

Смотрите также

полисоответствие