备战数学建模31-数据插值与曲线拟合3

Posted 2022-06-27 nuist__NJUPT

目录

一、数据插值经典案例

二、曲线拟合经典案例

1-多项式拟合

2-非线性拟合

---

一、数据插值经典案例

下面是一号池随着1，3，5...15周的各种指标的情况，对该池子中的指标做出插值，并将插值后绘制的所有的图像显示在一幅图上。

数据较多，如果选择每一行分别插值比较麻烦，可以循环插值，然后subplot函数绘制多个图，当然，对于一维插值，我们首选三次样条插值。

具体的MATLAB代码如下：

clear; clc
load('Z.mat')
x = Z(1,:) ; %得到一个行向量，周数
[n,m] = size(Z) ; %得到原始数据的行数和列数
%每幅图y轴的标签
ylab='周数','轮虫','溶氧','COD','水温','PH值','盐度','透明度','总碱度','氯离子','透明度','生物量';
P = zeros(11,15) ;

for i = 2 : n %从第2行开始插值
    y = Z(i,:) ;
    new_x = 1:15 ;
    p1 = interp1(x,y,new_x,'spline') ; %三次样条插值
    subplot(4,3,i-1) ;
    plot(x,y,'ro',new_x,p1,'-');
    axis([0 15,-inf,inf]) ;
    xlabel('周数') ;
    ylabel(ylab(i)) ;
    P(i-1,:) = p1 ;
end
legend('原始数据','三次样条插值数据','Location','SouthEast'); %图例
disp(P) ;

    

绘制的图形如下所示：

二、曲线拟合经典案例

1-多项式拟合

我们看一下简单例子，下面数据进行一次多项式拟合，当然可以直接拟合求解，也就是求解一元线性回归方程，也可以最小二乘法求解，我们除了绘制拟合曲线，也计算了拟合优度。

MATLAB代码如下：
方法1，直接一次多项式拟合

clear; clc
load('nihe.mat');
%plot(x,y,'o') ;
xlabel('x的值') ;
ylabel('y的值') ;
n = size(x,1) ; %样本点个数
x = x' ;
y = y' ;
xp = 2.7 : 0.1 : 6.9 ;
p = polyfit(x,y,1) ;
yp = polyval(p,xp);
plot(x,y,'o',xp,yp) ;
legend('样本数据','拟合函数','location','SouthEast') ;

方法2，最小二乘法参数估计，然后根据表达式绘图，此处计算了拟合优度。

clear; clc
load('nihe.mat');
plot(x,y,'o') ;
xlabel('x的值') ;
ylabel('y的值') ;
n = size(x,1) ; %样本点个数
%最小二乘法进行参数估计
k = (n*sum(x.*y)-sum(x)*sum(y))/(n*sum(x.*x)-sum(x)*sum(x)) ;
b = (sum(x.*x)*sum(y)-sum(x)*sum(x.*y))/(n*sum(x.*x)-sum(x)*sum(x)) ;
hold on ;
grid on ;
f = @(x) k*x + b ; %匿名函数
fplot(f,[1,7.5]) ;
legend('样本数据','拟合函数','location','SouthEast') ;
y_hat = k * x + b ;
SSR = sum((y_hat-mean(y)).^2) ; %回归平方和
SSE = sum((y_hat-y).^2) ; %误差平方和
SST = sum((y-mean(y)).^2) ; %总体平方和
disp('拟合优度值如下：') ;
R_2= SSR / SST ;
disp(R_2) ;

关于拟合优度的解释如下：拟合优度越接近于1，误差平方和越小，说明拟合效果越好。 

绘制的图形如下：

2-非线性拟合

我们看一下羡慕的人口预测，可以直接使用非线性拟合函数，也可以使用拟合工具箱cftool

MATLAB代码如下，此处直接进行非线性拟合，需要用到匿名函数。

clear; clc
x = 1790:10:2000;
y = [3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.0,92.0,106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4,281.4];
plot(x,y,'o')
y1 = @(b,t) b(1) ./ (1 + (b(1) / 3.9 - 1)*exp(-b(2)*(t-1790))); %匿名函数
b0 = [100 0.1] ; %初值的选取很重要
a = nlinfit(x,y,y1,b0)  ; %非线性拟合后的系数
xp = 1790 : 10 : 2030 ;
yp = y1(a,xp) ; %将系数和数值代入表达式求值
hold on ;
plot(xp, yp, '-') ;
legend('原始散点','拟合曲线') ;

绘制图形如下：