数据预处理——拟合matlab实现1

ゞ浴缸里的玫瑰 2022-12-23 13:30 471阅读 0赞

笔记整理来自《清风数学建模系列课程》

复习：[数据预处理——插值算法matlab实现][matlab]

**目录**

一、插值和拟合的区别

二、拟合——找到y和x之间的拟合曲线

1. 导入数据到matlab

2. 散点图分析

3. 确定拟合曲线的函数关系式

（1）最小二乘法的几何解释

（2）求解最小二乘法

（3）将散点图和拟合函数图像画在一起

(4) 评价拟合的好坏

（5）拟合优度的matlab求解

--------------------

# 一、插值和拟合的区别 #

与插值问题不同，在拟合问题中不需要曲线一定经过给定的点。拟问题的目标是寻求一个函数（曲线），使得该曲线在某种准则下与所有的数据点最为接近，即曲线拟合的最好（最小化损失函数）。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70][]

插值算法中，得到的多项式f(x)要**经过所有样本点**。但是如果样本点太多，那么这个多项式次数过高，会造成龙格现象。  
尽管我们可以选择分段的方法避免这种现象，但是更多时候我们更倾向于**得到一个确定的曲线**，尽管这条曲线**不能经过每一个样本点，但只要保证误差足够小**即可，这就是拟合的思想。(拟合的结果是得到一个确定的曲线)

# 二、拟合——找到y和x之间的拟合曲线 #

数据：

## 1. 导入数据到matlab ##

matlab中如何导入数据：

现在数据窗口新建两个变量x和y,把变量的数据从excel中粘贴进去。然后右击鼠标选择“save”，将数据保存为.mat格式即可，这里命名为data1.mat.

![20201125145501635.png][]

易错的细节：load data1.mat时导入出错是因为当前文件夹里面找不到data1.mat，所以需要把当前文件夹切换到有data1.mat的文件夹中。

![20201125150006914.png][]

## 2. 散点图分析 ##

clear;clc
    load  data1
    plot(x,y,'o')
    % 给x和y轴加上标签
    xlabel('x的值')
    ylabel('y的值')

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 1][]

通过散点图可以发现，数据变化的趋势有点像一次函数，所以我们可以定义y=f(x)的函数图像为y=kx+b,接下来只需要算出k和b的值即可。

## 3. 确定拟合曲线的函数关系式 ##

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 2][]

如何让样本点和拟合曲线最为接近呢？——最小二乘法。

### **（1）最小二乘法的几何解释** ###

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 3][]

![20201125151011547.png][]argmin表示使得|真实值与拟合值之差|的和最小的k和b的值。

第一种定义有绝对值，不容易求导，因此计算比较复杂。所以我们往往使用第二种定义，这也正是最小二乘的思想。最小二乘法的“二”其实指的就是平方的意思。  
为什么不用四次方？  
（1）避免极端数据对拟合曲线的影响。比如有一个异常值，如果用它的值减去拟合值再求四次方，会对结果有很大的影响。  
（2）最小二乘法得到的结果和MLE极大似然估计一致。  
不用奇数次方的原因：误差会正负相抵。

### （2）求解最小二乘法 ###

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 4][]

matlab求解：

![20201125153535654.png][]

size()函数的用法：

size（）：获取矩阵的行数和列数  
1、s=size(A),当只有一个输出参数时，返回一个行向量，该行向量的第一个元素是矩阵的行数，第二个元素是矩阵的列数。  
2、\[r,c\]=size(A),当有两个输出参数时，size函数将矩阵的行数返回到第一个输出变量r，将矩阵的列数返回到第二个输出变量c。  
3、size(A,n)如果在size函数的输入参数中再添加一项n，并用1或2为n赋值，则 size将返回矩阵的行数或列数。其中r=size(A,1)该语句返回的时矩阵A的行数， c=size(A,2) 该语句返回的时矩阵A的列数。

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

matlab计算得到：k = 2.0948，b =-1.0548

### （3）将散点图和拟合函数图像画在一起 ###

在指导拟合函数y=kx+b的函数关系式为y=2.0948x-1.0548,那么怎么画函数图像呢？

**方法1：我们可以自己生成x的值，再根据y=2.0948x-1.0548的值画图**

clear;clc
    load  data1
    plot(x,y,'o')
    % 给x和y轴加上标签
    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 % 显示网格线
    
    % % 画出y=kx+b的函数图像 plot(x,y)
    % % 传统的画法：模拟生成x和y的序列，比如要画出[0,5]上的图形
    xx = 2.5: 0.1 :7  % 间隔设置的越小画出来的图形越准确
    yy = k * xx + b  % k和b都是已知值
    plot(xx,yy,'-')

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 5][]

clear;clc
    load  data1
    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))
    grid on % 显示网格线
    % % 画出y=kx+b的函数图像 plot(x,y)
    xx = 2.5: 0.1 :7  % 间隔设置的越小画出来的图形越准确
    yy = k * xx + b  % k和b都是已知值
    plot(x,y,'o',xx,yy,'r-')
    legend('真实值','拟合值（拟合函数）')
    % 给x和y轴加上标签
    xlabel('x的值')
    ylabel('y的值')

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 6][]

**方法2： 使用匿名函数**

% 匿名函数的基本用法。

handle = @(arglist) anonymous_function

% 其中handle为调用匿名函数时使用的名字。  
% arglist为匿名函数的输入参数（自变量），可以是一个，也可以是多个，用逗号分隔。  
% anonymous\_function为匿名函数的表达式。

举例：如要表示z=x^2+y^2

z=@(x,y) x^2+y^2; 
    z(1,2)

% % ans = 5  
% fplot函数可用于画出匿名一元函数的图形。

fplot(f,xinterval)

% fplot(f,xinterval) 将匿名函数f在指定区间xinterval绘图。xinterval = \[xmin xmax\] 表示定义域的范围

clear;clc
    load  data1
    plot(x,y,'o')
    % 给x和y轴加上标签
    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,[2.5,7]);
    legend('样本数据','拟合函数','location','SouthEast')

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 7][]

### (4) 评价拟合的好坏 ###

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 8][]

证明SST = SSE + SSR

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 9][]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 10][]

> 这里的“线性”指的是“对参数为线性”。y=a+bx^2有两个参数a和b,并没有出现a\*b，a^2等参数非线性的情况，所以这个函数为线性函数。

如何判断线性于参数的函数？

> 在函数中，参数仅以一次方出现，且不能乘以或除以其他任何的参数，并不能出现参数的复合函数形式。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 11][]

神经网络可以用拟合优度。

### （5）拟合优度的matlab求解 ###

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 12][]

y_hat = k*x+b; % y的拟合值
    SSR = sum((y_hat-mean(y)).^2)  % 回归平方和
    SSE = sum((y_hat-y).^2) % 误差平方和
    SST = sum((y-mean(y)).^2) % 总体平方和
    SST-SSE-SSR   % 5.6843e-14  =   5.6843*10^-14   matlab浮点数计算的一个误差
    R_2 = SSR / SST

计算得到R^2=0.9635

[matlab]: https://liangyan.blog.csdn.net/article/details/109992809
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70]: /images/20221120/a75198cd436b4dbf9919b1eba2374cdb.png
[20201125145501635.png]: /images/20221120/d30b1e69749f4db6a365f2bca32c603c.png
[20201125150006914.png]: /images/20221120/482a9effb11b4f02840d9e5919718978.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 1]: /images/20221120/acf64dcb22914f718877d661ce0a3886.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 2]: /images/20221120/2183f089055f4833b0a8e1946ab4bb46.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 3]: /images/20221120/7f234994d28640099b31e42a3586dabd.png
[20201125151011547.png]: /images/20221120/b82d0a6943fd4899a3e0a0b61f4ecb4f.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 4]: /images/20221120/1bd2468061454ff69f667a6c6f972801.png
[20201125153535654.png]: /images/20221120/a07295fc59d74009bcb06c8f475d3014.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 5]: /images/20221120/d146d4d5c0fd4ead830a1ddae40256a9.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 6]: /images/20221120/e944d74a95894b1a9a319f04103b087b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 7]: /images/20221120/4a60f8b85f894fa8bc73d6bb2984ee85.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 8]: https://img-blog.csdnimg.cn/2020112515571362.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc=,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 9]: https://img-blog.csdnimg.cn/20201125160037698.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc=,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 10]: https://img-blog.csdnimg.cn/20201125161117971.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc=,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 11]: https://img-blog.csdnimg.cn/20201125161421861.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc=,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc_size_16_color_FFFFFF_t_70 12]: https://img-blog.csdnimg.cn/20201125161820387.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3lhbnlhbndlbm1lbmc=,size_16,color_FFFFFF,t_70