% 最小二乘法擬合程序說明

% 這是一個集交互式和命令式于一體的一個程序。

% 首先程序提示你輸入數據，如果你輸入的數據正確，則程序繼續進行。

% 如果你輸入的x,y的大小不一致，程序會提示你輸入數據有誤，然后程序返回再次讓你輸入數據。

%待你輸入數據正確后，程序會給你一段提示信息，

    % 通過下面的交互式圖形，你可以事先估計一下你要擬合的多項式的階數，方便下面的計算

    % polytool()是交互式函數,在圖形上方[Degree]框中輸入階數,右擊左下角的[Export]輸出圖形

    % 回車打開polytool交互式界面

% 待你回車后，交互式圖形界面打開，你可以在圖形上方[Degree]框中輸入階數觀察多項式擬合的圖形,

% 選擇置信區間最小的多項式階數，以便在下面使用。

% 待你選擇好所要擬合的多項式階數后，程序提示你回車繼續進行擬合。

% 待你回車后，程序提示你輸入擬合的階數，

% 待你輸入擬合的階數后，回車后會顯示擬合多項式相關的數據。

% 最后程序提示你輸入所需擬合的數據點，輸入所需擬合的數據點，就會得到相應的結果。

% 如果沒有所需要擬合的數據點，直接回車跳過。

% 程序結束。
% 最小二乘法多項式擬合

% (1)-----選擇擬合多項式擬合的階數
disp('請以向量的形式輸入x,y.')
x=input('x=');
y=input('y=');
nx = length(x);
ny = length(y);
n  = length(x); 
if nx == ny
x1 = x(1); xn = x(n);
% n個數據可以擬合(n-1)階多項式,高階多項式多次求導,數值特性變差

disp('通過下面的交互式圖形，你可以事先估計一下你要擬合的多項式的階數，方便下面的計算.')

disp('polytool()是交互式函數,在圖形上方[Degree]框中輸入階數,右擊左下角的[Export]輸出圖形')

disp('回車打開polytool交互式界面')

pause;

polytool(x,y,1)

% 觀察多項式擬合的圖形,選擇置信區間最小的多項式階數

disp('回車繼續進行擬合')

pause;

% (2)-----計算多項式的各項系數和擬合值
m=input('    輸入多項式擬合的階數   m = ');
[p,S]=polyfit(x,y,m);
disp '      輸出多項式的各項系數'
fprintf (1,'         a = %3.16f \n',p)
disp '      輸出多項式的有關信息 S'
disp (S)
[yh,delta]=polyconf(p,x,S);
disp '         觀測數據     擬合數據'
disp '       x         y       yh'
for i = 1 : n
   xy = [x(i) y(i) yh(i)];
   disp (xy)
end
% (3)-----繪制觀測數據離散點圖和多項式曲線

plot(x,y,'r.')

title('\bf 實驗數據離散點圖  /  多項式曲線 \it y = a0+a1x+a2x^2+a3x^3+...')
grid
hold on; 
xi=[x1:0.1:xn];
yi=polyval(p,xi);
plot(xi,yi,'k-')
% (4)-----擬合效果和精度檢驗
Q=sum((y-yh).^2);
SGM = sqrt(Q / (n - 2));
RR = sum((yh-mean(y)).^2)/sum((y-mean(y)).^2);
fprintf (1,'      剩余平方和            Q = %3.6f \n',Q)
fprintf ('\n')
fprintf (1,'        標準誤差        Sigma = %3.6f \n',SGM)
fprintf ('\n')
fprintf (1,'        相關指數           RR = %3.6f \n',RR)
fprintf ('\n')
disp('請輸入你所需要擬合的數據點，若沒有請按回車鍵結束程序.')
fprintf ('\n')
x0=input('    輸入插值點             x0 = ');
y0=polyval(p,x0);
fprintf (1,'    輸出插值點擬合函數值   y0 = %3.4f \n',y0)
else
disp('輸入的數據有誤，請重新運行程序并輸入正確的數據。')
clear
zxecf
end