Matlab - Schlecht konditionierte Polynom

Question

Ich habe eine Frage mich zu fragen, lineare, quadratische und kubische fit für Bevölkerungsdaten zu tun, dann schätzen die Bevölkerung im Jahr 1915. Die lineare und quadratische fits arbeiten, aber die kubische scheint einen Fehler zu erhöhen mir zu sagen, das Polynom ist schlecht konditioniert. Die Grafik ist ziemlich nah an den Datenwerten und scheint eine gute Anpassung zu sein. Was kann ich tun, um dies zu beheben? Der Code ist:

clear; 
clc; 
close all; 

year = [1815,1845,1875,1905,1935,1965]; 
population = [8.3,19.7,44.4,83.2,127.1,190.9]; 

rlinear = polyfit(year,population,1); 
rquadratic = polyfit(year,population,2); 
rcubic = polyfit(year,population,3); 

newTime = linspace(1815,1965,100); 
vrlinear = polyval(rlinear,newTime); 
vrquadratic = polyval(rquadratic,newTime); 
vrcubic = polyval(rcubic,newTime); 

subplot(2,2,1) 
plot(year,population,'ob',newTime,vrlinear) 
xlabel('Year') 
ylabel('Population (millions)') 
title('Year vs. US population') 

subplot(2,2,2) 
plot(year,population,'ob',newTime,vrquadratic) 
xlabel('Year') 
ylabel('Population (millions)') 
title('Year vs. US population') 

subplot(2,2,3) 
plot(year,population,'ob',newTime,vrcubic) 
xlabel('Year') 
ylabel('Population (millions)') 
title('Year vs. US population') 

estimate = polyval(rquadratic,1915); 
fprintf('The estimated population in the year 1915 is %d million. \r',estimate) 

Answer 1

Im Anschluss an die Warnmeldung:

Warnung: Polynomial ist schlecht konditioniert. In Punkte mit unterschiedlichen X Werte, reduzieren den Grad des Polynoms, oder versuchen, Zentrieren und Skalierung wie in HELP POLYFIT beschrieben.

Zentrieren und Skalieren löst das Problem:

[rcubic2,~,mu] = polyfit(year,population,3); 
vrcubic2 = polyval(rcubic2,newTime,[],mu); 

subplot(2,2,3) 
plot(year,population,'ob',newTime,vrcubic1) 
hold on 
plot(newTime,vrcubic2,'--r') 
xlabel('Year') 
ylabel('Population (millions)') 
title('Year vs. US population') 

Die Ergebnisse zeigen, dass in diesem Fall die beiden kubischen passt praktisch gleich sind. Siehe polyfit Hilfe für weitere Details zu diesem Thema.

Answer 2

In neueren Versionen von MATLAB, die Warnung sagt auch, sollten Sie read the help. Insbesondere müssen Sie Ihre Daten neu skalieren, um numerische Probleme zu vermeiden. Glücklicherweise ist dies etwas, das polyfit und polyval für Sie tun kann.

Der Grund, warum Sie dies tun müssen, ist, dass die Anpassung eines Polynoms ziemlich schlecht konditioniert ist. Hinter den Kulissen baut es eine Matrix auf, die Ihre Jahre enthält, die zu jeder Potenz im Polynom erhoben werden. Sie haben also eine Matrix mit Einträgen wie 1815 zusammen mit 1815^3 = 6e9, das ist nicht ideal für die Konditionierung. Die Skalierung stellt sicher, dass Sie keine großen Zahlen in der Matrix erhalten und somit die Konditionsnummer verbessert wird.

Praktisch haben Sie polyfit und polyval ein bisschen anders zu nennen:

[p,s,m] = polyfit(year, population, order) 
[populationAtOtherYears] = polyval(p, otherYears, s, m) 

So Ihr Skript etwas aussehen würde:

clear; 
clc; 
close all; 

year  = [1815,1845,1875,1905,1935,1965]; 
population = [8.3,19.7,44.4,83.2,127.1,190.9]; 

[linear.poly, linear.sigma, linear.mu] = polyfit(year,population,1); 
[quadratic.poly, quadratic.sigma, quadratic.mu] = polyfit(year,population,2); 
[cubic.poly, cubic.sigma, cubic.mu] = polyfit(year,population,3); 

newTime = linspace(1815,1965,100); 
% next line is to make a function for easier calling 
evaluatePolynomial = @(fit, t) polyval(fit.poly, t, fit.sigma, fit.mu); 
linear.evaluated = evaluatePolynomial(linear, newTime); 
quadratic.evaluated = evaluatePolynomial(quadratic, newTime); 
cubic.evaluated  = evaluatePolynomial(cubic, newTime); 

subplot(2,2,1) 
plot(year,population,'ob',newTime,linear.evaluated) 
xlabel('Year') 
ylabel('Population (millions)') 
title('Year vs. US population') 

subplot(2,2,2) 
plot(year,population,'ob',newTime,quadratic.evaluated) 
xlabel('Year') 
ylabel('Population (millions)') 
title('Year vs. US population') 

subplot(2,2,3) 
plot(year,population,'ob',newTime,cubic.evaluated) 
xlabel('Year') 
ylabel('Population (millions)') 
title('Year vs. US population') 

estimate = evaluatePolynomial(quadratic,1915); 
fprintf('The estimated population in the year 1915 is %d million. \r',estimate)