matlab 趋势面分析 记录

1. 利用matlab进行趋势面分析并绘图

%% 方法1 仅二次趋势面模型

clc;clear;

num = xlsread("G:\Data\灾损曲线\趋势面分析-数据.xlsx",'B21:N37');

volume = num(:,1);

magnitude = num(:,3);

pop = num(:,8);

gdp = num(:,13);

% 绘制三维散点图

scatter3(volume, magnitude, gdp,'filled')

hold on;

% 估计模型，得到估计参数

v2 = volume.^2;

m2 = magnitude.^2;

vm = volume.*magnitude;

X = [ones(length(volume),1) volume magnitude v2 vm m2];

[b,bint,r,rint,stats] = regress(gdp,X,95);

% 生成格点经纬度，计算预测值

xfit = min(volume):0.1:max(volume);

yfit = min(magnitude):0.1:max(magnitude);

[XFIT,YFIT]= meshgrid (xfit,yfit);

ZFIT = b(1) + b(2) * XFIT + b(3) * YFIT + b(4) * XFIT.^2 + b(5) * XFIT .* YFIT + b(6) * YFIT.^2;

% 绘制拟合曲面

mesh(XFIT,YFIT,ZFIT);

clc;clear;
num = xlsread("G:\Data\灾损曲线\趋势面分析-数据.xlsx",'B21:N37');
volume = num(:,1);
magnitude = num(:,3);
pop = num(:,8);
gdp = num(:,13);

% 绘制三维散点图
scatter3(volume, magnitude, gdp,'filled')
hold on;
% 估计模型，得到估计参数
v2 = volume.^2;
m2 = magnitude.^2;
vm = volume.*magnitude;
X = [ones(length(volume),1) volume magnitude v2 vm m2];
[b,bint,r,rint,stats] = regress(gdp,X,95);
% 生成格点经纬度，计算预测值
xfit = min(volume):0.1:max(volume);
yfit = min(magnitude):0.1:max(magnitude);
[XFIT,YFIT]= meshgrid (xfit,yfit);
ZFIT = b(1) + b(2) * XFIT + b(3) * YFIT  + b(4) * XFIT.^2 + b(5) * XFIT .* YFIT + b(6) * YFIT.^2;
% 绘制拟合曲面
mesh(XFIT,YFIT,ZFIT);

%% 方法2  一次、二次及三次趋势面模型

% 文件信息读入

clc;clear;

num = xlsread("G:\Data\灾损曲线\趋势面分析-数据.xlsx",'B21:N37');

volume = num(:,1);

magnitude = num(:,3);

pop = num(:,8);

gdp = num(:,13);

opoX=num(:,1);

opoY=num(:,3);

oPH=num(:,13);

% oOM=info(:,4);

% oTN=info(:,5);

%{

% 2S法计算异常值

mPH=mean(oPH);

sPH=std(oPH);

num2=find(oPH>(mPH+2*sPH)|oPH<(mPH-2*sPH));

% 异常值剔除

PH=oPH;

for i=1:length(num2)

n=num2(i,1);

PH(n,:)=[0];

end

PH(all(PH==0,2),:)=[];

poX=opoX;

for i=1:length(num2)

n=num2(i,1);

poX(n,:)=[0];

end

poX(all(poX==0,2),:)=[];

poY=opoY;

for i=1:length(num2)

n=num2(i,1);

poY(n,:)=[0];

end

poY(all(poY==0,2),:)=[];

%}