數學建模多元分析實例
原文:http://blog.csdn.net/qq_34861102/article/details/77101690
聚類分析
相關命令
pdist
Y = pdist(X,'metric')
X
必選,'metric'
指定計算矩陣中對象的距離的方式pdist(X,'minkowski',p)
p
在閔氏距離計算矩陣中用到的指數值linkage
Z = linkage(Y,'method')
Y
上面生成的距離(必選),method
按照指定的算法生成聚類樹(可選)cluster
T = cluster(Z,'cutoff',c)
zsore
zsore(X)
dendrogram
H = dendrogram(Z,P)
P
為節點數,默認是30clusterdata
squareform
pdist
轉化為方陣實例
a = [1,0;1,1;3,2;4,3;2,5];
y = pdist(a,'cityblock');
yc = squareform(y);
z = linkage(y);
dendrogram(z);
T = cluster(z,'maxclust',3)
主成分分析(降低維度)
使用步驟
R
p
個主成分,計算綜合評價值使用實例:
%獲取數據
gj = load('data.txt');
%數據標准化
gi = zscore(gj);
%計算相關系數矩陣
r = corrcoef(gj);
%x的列為r的特征向量,即主成分的系數
%y為r的特征值
%z為各個主成分的貢獻率
[x,y,z] = pcacov(r);
%構造與x矩陣同維數的元素為+1或-1的矩陣
f = repmat(sign(sum(x)),size(x,1),1);
%使得每個特征向量的分量和為正
x = x.*f;
%選取主成分的個數
num = 3;
%計算各個主成分的得分
df = gj*x(:,1:num);
%計算綜合得分
tf = df*z(1:num)/100;
[stf,ind] = sort(tf,'descend');
stf = stf';
ind = ind';
因子分析
R
p
個主因子使用實例
%獲取數據
gj = load('data.txt');
%數據標准化
gi = zscore(gj);
%計算相關系數矩陣
r = corrcoef(gj);
%x的列為r的特征向量,即主成分的系數
%y為r的特征值
%z為各個主成分的貢獻率
[x,y,z] = pcacov(r);
%構造與x矩陣同維數的元素為+1或-1的矩陣
f = repmat(sign(sum(x)),size(x,1),1);
%使得每個特征向量的分量和為正
x = x.*f;
%初等荷載矩陣
f2 = repmat(sqrt(y)',size(x,1),1);
a = x.*f2;
%主因子個數
num = input('請選擇主因子的個數:');
am = a(:,[1,num]);
[b,t] = rotatefactors(am,'method','varimax');
bt = [b,a(:,[num+1:end])];
%計算共同度
degree = sum(b.^2,2);
%計算因子貢獻
contr = sum(bt.^2);
%計算因子貢獻率
rate = contr(1:num)/sum(contr);
%計算得分函數的系數
coef = inv(r)*b;
%計算得分
score = x*coef;
判別分析(判別個體所屬類別)
實例:
a = [9 7 8 8 9 8 7 4 3 6 2 1 6 8 2;8 6 7 5 9 9 5 4 6 3 4 2 4 1 4;7 6 8 5 3 7 6 4 6 3 5 2 5 3 5];
train = a(:,[1:12])';
sample = a(:,[13:end])';
group = [ones(7,1);2*ones(5,1)];
%馬氏距離分類
[x1,y1] = classify(sample,train,group,'mahalanobis');
%線性分類
[x2,y2] = classify(sample,train,group,'linear');
%二次分類
[x3,y3] = classify(sample,train,group,'quadratic');
典型相關分析
x = load x.txt;
y = load y.txt;
p = size(x,2);
q = size(y,2);
%標准化數據
x = zscore(x);
y = zscore(y);
%觀測數據的個數
n = size(x,1);
%a1,b1返回的是典型變量的系數。r返回的是典型相關系數
%u1,v1返回的是典型變量的值。status返回的是假設檢驗的一些統計量的值
[a1,b1,r,u1,v1,status] = cononcorr(x,y);
%修正每一列的正負號,使系數和為正
a = a1.*repmat(sign(sum(a1)),size(a1,1),1);
b = b1.*repmat(sign(sum(b1)),size(b1,1),1);
u = u1.*repmat(sign(sum(u1)),size(u1,1),1);
v = v1.*repmat(sign(sum(v1)),size(v1,1),1);
%計算x,u y,v x,v x,u的相關系數
x_u_r = x'*u/(n-1);
y_v_r = y'*v/(n-1);
x_v_r = x'*v/(n-1);
y_u_r = y'*u/(n-1);
%X組原始變量被u_i解釋的方差比例
%方差累積比例
ux = sum(x_u_r.^2)/p;
ux_cum = cumsum(ux);
vx = sum(x_v_r.^2)/p;
vx_cum = cumsum(vx);
uy = sum(y_u_r.^2)/q;
uy_cum = cumsum(uy);
vy = sum(y_v_r.^2)/q;
vy_cum = cumsum(vy);
%典型相關系數的平方
val = r.^2;
對應分析(Q型和R型的結合)
a = [543 342 453 609 261 360 243 183;245 785 630 597 311 233 108 69;300 200 489 740 365 324 327 228;401 396 395 693 350 309 263 143;147 117 410 726 366 447 329 420];
%行和
a_i_dot = sum(a,2);
%列和
a_dot_j = sum(a);
%數據和
T = sum(a_i_dot);
%對應矩陣
P = a/T;
%邊緣分布
r = sum(P,2);
c = sum(P);
%計算行輪廓分布陣
Row_prifile = a./repmat(sum(a,2),1,size(a,2));
%標准化數據陣
B = (P - r*c) ./ sqrt((r*c));
%對標准化后的數據陣作奇異值分解
[u s v] = svd(B,'econ');
%修改特征向量的符號矩陣
w1 = sign(repmat(sum(v),size(v,1),1));
%使得v中每一個行向量的分量和大於0
w2 = sign(repmat(sum(v),size(u,1),1));
%修改特征向量的正負號
vb = v.*w1;
ub = u.*w2;
%計算慣量
lamda = diag(s).^2;
%計算卡方統計量的分解
ksi2square = T*(lamda);
%計算總卡方統計量
T_ksi2square = sum(ksi2square);
%計算貢獻率
con_rate = lamda/sum(lamda);
%計算累積貢獻率
cum_rate = cumsum(con_rate);
%求加權特征向量
beta = diag(r.^(-1/2))*ub;
%求行輪廓坐標
G = beta*s;
%求加權特征向量
alpha = diag(c.^(-1/2))*vb;
%求列輪廓坐標
F = alpha*s;
%樣本點的個數
num1 = size(G,1);
%行坐標的取值范圍
rang = minmax(G(:,[1,2])');
%畫圖調整
delta = (rang(:,2)-rang(:,1))/(4*num1);
chrow = {'A','B','C','D','E'};
strcol = {'少男','少女','白領','工人','農民','士兵','主管','教授'};
plot(G(:,1),G(:,2),'*','Color','k','LineWidth',1.3);
text(G(:,1),G(:,2)-delta(2),chrow);
hold on;
plot(F(:,1),F(:,2),'H','Color','k','LineWidth',1.3);
text(F(:,1)-delta(1),F(:,2)+1.2*delta(2),strcol);
xlabel('dim1');
ylabel('dim2');
多維標度分析
D = [0 1 sqrt(3) 2 sqrt(3) 1 1; zeros(1,2) 1 sqrt(3) 2 sqrt(3) 1; zeros(1,3) 1 sqrt(3) 2 1; zeros(1,4) 1 sqrt(3) 1; zeros(1,5) 1 1;zeros(1,6) 1;zeros(1,7)];
d = D + D';
[y,eigvals] = cmdscale(d);
plot(y(:,1),y(:,2),'o','Color','k','LineWidth',1.3)
本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。