中南大学实用测量数据处理实验

中南大学《实用测量数据处理》实验（自用）

实验一：总体最小二乘估计

为了确定直线方程$y={\beta }_0+{\beta }_{1}x$中的两个参数，获取了四对观测值$(1,2),(2,6),(4,4),(6,1)$，试根据最小二乘和总体最小二乘求出参数估值。

A = [1 1;1 2;1 4;1 6];
L = [2;6;4;1];

X0 = (A'*A)\A'*L;
X1 = X0;
while true

    V1 = (A*X1-L)'*(A*X1-L)/(1+X1'*X1);
    X2 = (A'*A)\A'*L+(A'*A)\X1*V1;
    
    if all(abs(X2-X1)<10e-4)
        break;
    else
        X1 = X2;
    end
end
% 可视化
figure
scatter(A(:,2),L)
hold on
x = 0:7;
plot(x,X0(1)+X0(2)*x)
plot(x,X1(1)+X1(2)*x)
legend("观测点","最小二乘","整体最小二乘")

实验二：线性回归分析

序号	坝段1温度/(℃)	坝段2温度/(℃)	水平位移/mm
1	6.81	6.57	15.40
2	9.88	9.53	13.20
…	…	…	…

导入数据

clear;close all;clc;

data=[6.81	6.57	15.40
    9.88	9.53	13.20
    9.52	7.79	15.26
    9.70	12.32	11.33
    6.52	9.88	13.32
    5.34	8.26	14.26
    9.50	7.55	14.16];

A = [ones(size(data,1),1),data(:,1:end-1)];
y = data(:,end)

计算回归结果

% 回归
beta = (A'*A)\A'*y
y_ = A*beta

显著性检验

m = size(A,2)-1;
n = size(A,1);
% 1、总体显著性检验
% 原假设：beta为0
avg = mean(y);

SSR = sum((y_ - avg).*(y_ - avg));  % 回归
SST = sum((y_ - y).*(y_ - y));  % 残差

F = (SSR/m)/(SST/(n-m-1));
if F>finv(0.950,m,n-m-1)
    disp("拒绝原假设，回归方程是显著的。")
end

% 2、回归系数显著性检验
% (1)使用t检验
C = diag(inv(A'*A));
sigma0 = sqrt(SST/(n-m-1));
for i=2:size(C,1)
    t = beta(i)/(sigma0*sqrt(C(i)));
    if abs(t)>tinv(0.975,n-m-1)
        fprintf("自变量X%d对y显著的。\n",i)
    else
        fprintf("自变量X%d对y不显著\n",i)
    end
end
% (2)使用f检验
for i=2:size(C,1)
    F = (beta(i)*beta(i)/C(i))/(SST/(n-m-1));
    if F>finv(0.950,1,n-m-1)
        fprintf("自变量X%d对y显著的。\n",i)
    else
        fprintf("自变量X%d对y不显著\n",i)
    end
end

实验三：非线性回归分析

导入数据

clc;clear;close all;
% 导入数据
y = load("data3.txt");

非线性回归

A = ones(size(y,1),4);
t = (1:size(y,1))';

% 对θ设置初始值
theta1 = [max(y),max(y),1,0]'; 

% 开始迭代
while true
    
    A(:,2) = -exp(-theta1(3)*(t.^theta1(4)));
    A(:,3) = -theta1(2)*A(:,2).*(t.^theta1(4));
    A(:,4) = A(:,3).*theta1(3).*log(t); 

    L = y - (f(t,theta1));

    % 平差,防止奇异使用岭估计
    dtheta = (A'*A + diag(ones(1,4)*0.001))\A'*L; 
    
    theta2 = theta1 + dtheta;
    
    % 迭代退出条件
    if all(abs(theta2-theta1)<10e-4)
        break;
    else
        theta1 = theta2;
    end
end

假设检验

% 假设检验
y_ = f(t,theta2);  % 在最后面定义了f函数
e = y - y_;
mu = mean(e);  % 均值
% sigma = std(e);  % 标准差
n = size(e,1);
sigma = sqrt(sum((e-mu).*(e-mu))/(n-1));

% 统计每个区间的频数
L = 10;
[m,~] = hist(e)  % m用于存储频数，注意hist只能统计10个区间。

l = min(e);r = max(e);d = (r-l)/L; % 间隔
x = l:d:r;
q = (normcdf(x(2:end),mu,sigma) - normcdf(x(1:end-1),mu,sigma))*n  % 理论频数

% 1\计算卡方统计量结果
kafang = sum((m-q).*(m-q)./q)
if abs(kafang)>chi2inv(0.975,L-3)
    disp("残差不服从正态分布")
else
    disp("残差服从正态分布")
end

% 2\残差均值是否为零的检验
U = mu/(sigma/sqrt(n))

if abs(U)<norminv(0.975,0,1)
    disp("残差均值服从为0的检验")
else
    disp("不服从")
end

function y=f(t,theta)
    y = theta(1)-theta(2)*exp(-theta(3)*t.^theta(4));
end

实验四：平滑与滤波

clc;clear;close all;
% 读取数据
y = load("data4.dat");

% 可视化数据
tiledlayout(5,1);
x = 1:size(y,1);
nexttile
plot(x,y)

% 三点滑动平均
g1 = (y(1:end-2) + y(2:end-1) + y(3:end))/3;
nexttile
plot(x(2:end-1),g1)

% 五点滑动平均（线性）
g2 = (y(1:end-4) + y(2:end-3) + y(3:end-2) + y(4:end-1) + y(5:end))/5;
nexttile
plot(x(3:end-2),g2)


% 五点滑动平均（二次曲线）
g3 = (-3*y(1:end-4) + 12*y(2:end-3) + 17*y(3:end-2) + 12*y(4:end-1) - 3*y(5:end))/35;
nexttile
plot(x(3:end-2),g3)


% 高斯权平滑
N = 5;a = 5;
st = (N-1)/2+1;ed = size(y,1)-(N-1)/2;
g4 = zeros(ed-st+1,1);
w = weight(N,a);  % 计算权重
for i=st:ed
    g4(i-st+1) = w'*y(i-2:i+2);
end
nexttile
plot(x(st:ed),g4)

function w = weight(N,sigma)
    n = (N-1)/2+1;
    w = zeros(N,1);
    for i=1:N
        w(i) = guassin(n-i,sigma);
    end
    w = w/sum(w);
end

function g = guassin(x,sigma)
    g = 1/sqrt(2*pi*sigma*sigma)*exp(-x*x/(2*sigma*sigma));
end