数学建模常用的分析法及其MATLAB实现

[toc]

介绍

本篇分别有主成分分析,层次分析,聚类分析,因子分析的简单介绍以及在MATLAB(2018a)上使用的简单介绍 ( 附源码 )

主成分分析

主成分分析-百度百科
简单讲:主成分分析就是取出几个较少的变量尽可能多地反映原来变量的信息
(只看主要对结果有影的几个变量)
主成分分析适用于数学建模中筛选出少量变量对结果造成主要影响并用函数关系表征的情况

%% 主成分分析 (降维)
clc, clear

load example_1.txt  %数据导入(数据要求:前几列为自变量,最后一列为因变量)
data = example_1;

[m, n] = size(data); 
num = 3;   % 选取的主成分的个数

mu = mean(data); 
sigma = std(data);  %标准差
%z-score标准化方法适用于属性A的最大值和最小值未知的情况,或有超出取值范围的离群数据的情况
%标准化的新数据=(原数据-均值)/标准差
std_data = zscore(data);
b = std_data(: , 1:end-1);     % 四个变量x1, x2, x3, x4
r = cov(b);                                % 变量的协方差矩阵

% 运用协方差矩阵进行PCA
[PC, latent, explained] = pcacov(r);  %返回主成分(PC)、协方差矩阵X的特征值 (latent)和每个特征向量表征在观测量总方差中所占的百分数(explained)
% 新的主成分z1 = PC(1,1)*x1 + PC(2,1)*x2 + PC(3,1)*x3 + PC(4,1)*x4  , z2 = ...
f = repmat(sign(sum(PC)), size(PC, 1), 1);            %sum(PC)表示对矩阵PC的列求和
PC = PC .* f;


%1.普通的最小二乘法回归
regress_args_b = [ones(m, 1), b] \ std_data(:, end);   %标准化数据的回归方程系数
bzh = mu ./ sigma;
% 原始数据的常数项
ch10 = mu(end) - bzh(1:end-1) * regress_args_b(2:end) * sigma(end);
fr_1 = regress_args_b(2:end); fr_1 = fr_1';
% 原始数据的自变量的系数
ch1 = fr_1 ./ sigma(1:end-1) * sigma(end);
% 此时模型为 y = ch10 + ch1[1]*x1 + ch1[2] * x2 + ch1[3] * x3 + ch1[4] * x4
% 计算均方误差
check1 = sqrt(sum( (data(:, 1:end-1) * ch1' + ch10 - data(:, end)) .^2 ) / (m - n))


%2.主成分回归模型
pca_val = b * PC(:, 1:num);
%主成分数据的回归方程系数
regress_args_pca = [ones(m, 1), pca_val] \ std_data(:, end);
beta = PC(:, 1:num) * regress_args_pca(2:num+1);   %标准化数据的回归方程系数
% 原始数据的常数项
ch20 = mu(end) - bzh(1:end-1) * beta * sigma(end);
fr_2 = beta';
% 原始数据的自变量的系数
ch2 = fr_2 ./ sigma(1:end-1) * sigma(end);
% 此时模型为 y = ch20 + ch2[1]*x1 + ch2[2] * x2 + ch2[3] * x3 + ch2[4] * x4
% 计算均方误差
check2 = sqrt(sum( (data(:, 1:end-1) * ch2' + ch20 - data(:, end)) .^2 ) / (m - num - 1))

层次分析

层次分析-百度百科
简单讲:就是自己将两两影响因素作重要性对比得出判断矩阵,然后得出各个影响因素的权重
层次分析适用于数学建模多个变量对结果不同影响程度并用函数关系表征的情况

%层次分析法(AHP)
disp('请输入判断矩阵A(n阶)');
A = input('A=');
[n,n] = size(A);
x = ones(n,100);
y = ones(n,100);
m = zeros(1,100);
m(1) = max(x(:,1));
y(:,1) = x(:,1);
x(:,2) = A*y(:,1);
m(2) = max(x(:,2));
y(:,2) = x(:,2)/m(2);
p=0.0001; i=2; k=abs(m(2)-m(1));
while k>p
    i=i+1;
    x(:,i) = A*y(:,i-1);
    m(i) = max(x(:,i));
    y(:,i) = x(:,i)/m(i);
    k=abs(m(i)-m(i-1));
end
a = sum(y(:,i));
w = y(:,i)/a;
t = m(i);
disp(w);
%一致性检验
CI = (t-n)/(n-1);
RI = [0 0 0.52 0.89 1.12 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];
CR = CI/RI(n);
if CR<0.10
    disp('此矩阵一致性可以接受!');
    disp('CI=');disp(CI);
    disp('CR=');disp(CR);
end

聚类分析

聚类分析-百度百科
简单讲:就是多个影响因素按不同联系程度分类

%聚类分析
X;   %X为导入数据的名称(类型:数值矩阵)
X=X';    %将矩阵X转置,按情况看是否需要此语句
Y=pdist(X);     %X为要聚类的数组
SF=squareform(Y);
Z=linkage(Y,'average');
dendrogram(Z,0);
T=cluster(Z,'maxclust',3);

因子分析

因子分析-百度百科
可看作是主成分分析的推广,涉及的计算与主成分分析很类似

%因子分析
clc,clear
load ssgs.txt   %把原始数据保存在纯文本文件ssgs.txt中
n=size(ssgs,1);
x=ssgs(:,[1:4]); y=ssgs(:,5); %分别提出自变量x1...x4和因变量x的值
x=zscore(x); %数据标准化
r=corrcoef(x)  %求相关系数矩阵
[vec1,val,con1]=pcacov(r)  %进行主成分分析的相关计算
f1=repmat(sign(sum(vec1)),size(vec1,1),1);
vec2=vec1.*f1;     %特征向量正负号转换
f2=repmat(sqrt(val)',size(vec2,1),1); 
a=vec2.*f2   %求初等载荷矩阵
num=input('请选择主因子的个数:');  %交互式选择主因子的个数
am=a(:,[1:num]);  %提出num个主因子的载荷矩阵
[bm,t]=rotatefactors(am,'method', 'varimax') %am旋转变换,bm为旋转后的载荷阵
bt=[bm,a(:,[num+1:end])];  %旋转后全部因子的载荷矩阵,前两个旋转,后面不旋转
con2=sum(bt.^2)       %计算因子贡献
check=[con1,con2'/sum(con2)*100]%该语句是领会旋转意义,con1是未旋转前的贡献率
rate=con2(1:num)/sum(con2) %计算因子贡献率
coef=inv(r)*bm          %计算得分函数的系数
score=x*coef           %计算各个因子的得分
weight=rate/sum(rate)  %计算得分的权重
Tscore=score*weight'   %对各因子的得分进行加权求和,即求各企业综合得分
[STscore,ind]=sort(Tscore,'descend')      %对企业进行排序
display=[score(ind,:)';STscore';ind'] %显示排序结果
[ccoef,p]=corrcoef([Tscore,y])    %计算F与资产负债的相关系数
[d,dt,e,et,stats]=regress(Tscore,[ones(n,1),y]);%计算F与资产负债的方程
d,stats  %显示回归系数,和相关统计量的值

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,265评论 6 490
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,078评论 2 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 156,852评论 0 347
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,408评论 1 283
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,445评论 5 384
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,772评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,921评论 3 406
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,688评论 0 266
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,130评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,467评论 2 325
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,617评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,276评论 4 329
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,882评论 3 312
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,740评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,967评论 1 265
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,315评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,486评论 2 348