Matlab-数据建模

Posted on | In , | Visitors:

This is my blog.
这一部分貌似和概率论与数理统计很有关系
唔～八大分布还记得…
怎么最近看的感觉我都认识
但是怎么看解析还是费劲呢？！！

数据建模

云模型

云数字特征

  1. 期望：\(E_x\)
  2. 熵： \(E_n\)不确定度
  3. 超熵： \(H_e\) 即熵的熵

触发机制与逆向触发机制

触发机制：

  1. 生成以En为期望，以\({H_e}^2\)为方差的正态随机数\(E_n’\)
  2. 生成以Ex为期望，以\({E_n’}^2\)为方差的正态随机数\(x\)
  3. 计算确定度\(u = e^{-\frac{(x-E_x)^2}{2{E_n’}^2}}\)
  4. 重复1 2 3，产生足够多的云滴

而通常我们是通过云滴来得到期望值的。
逆向触发机制如下：

  1. 计算样本均值\(\overline{x}\)和方差\(S^2\)
  2. \(E_x = \overline{x}\)
  3. \(E_n = \sqrt{\frac{\pi}{2}}*\frac{1}{n}\sum_{1}^n\left|x-E_x\right|\)
  4. \(H_e = \sqrt{S^2-{E_n}^2}\)

Matlab

could_transform.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
function [x,y,Ex,En,He] = could_transform(y_spor,n)
 
Ex = mean(y_spor);
En = mean(abs(y_spor - Ex)) .* sqrt(pi ./ 2);
He = sqrt(var(y_spor) - En.^2);
 
for q = 1:n
    Enn = randn(1).* He + En;
    x(q) = randn(1).* Enn + Ex;
    y(q) = exp(-(x(q) - Ex) .^ 2 ./ (2 .* Enn .^2));
end
x;
y;
end

test.m

1
2
3
4
5
6
7
8
9
10
N = 1500;
Y = [9.5 10.3 10.1 8.1
    10.3 9.7 10.4 10.1
    10.6 8.6 9.2 10.0]';
for i = 1:size(Y,1)
    subplot(size(Y,1)/2,2,i)    
    [x,y,Ex,En,He] = could_transform(Y(i,:),N);
    plot(x,y,'r.'); 
    axis([8,12,0,1];
end

Logistic回归

模型

\(ln\frac{\pi}{1-\pi} = \beta_0+\beta_1x_1+\cdot\cdot\cdot+\beta_kx_k\(
在Matlab中我们使用regress函数
感觉像是0-1分布，我瞎说的，别理我

Matlab

已知：
\(p\(表示评价结果 \(q\(表示某概率
\(P = 0,q\leq0.5\(
\(P = 1,q>0.5\(
为了方便做回归运算，取区间的中值作为\(q\(的值，即
\(P = 0,q=0.25\(
\(P = 1,q=0.75\(

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
X0 = []; %回归数据
XE = []; %验证与预测数据
Y0 = []; %回归数据P值
 
% 数据转化和参数回归
n = size(Y0,1);
for i = 1:n
    if Y0(i) == 0
        Y1(i,1) = 0.25;
    else
        Y1(i,1) = 0.75;
    end
end
 
X1 = ones(size(X0,1),1); %构建常数项系数
X = [X1,X0];
Y = log(Y1 ./ (1 - Y1));
b = regress(Y,X);
 
%模型的验证
for i = 1:size(XE,1)
    q0 = exp(b(1) + b(2) * XE(i,1) + b(3) * XE(i,2) + b(4) * XE(i,3))/(1+...
    exp(b(1) + b(2) * XE(i,1) + b(3) * XE(i,2) + b(4) * XE(i,3)));
    if q0 <= 0.5
        P(i) = 0;
    else
        P(i) = 1;
    end
end
 
disp(['回归系数'num2str(b')]);
disp(['评价结果'num2str(P)]);

主成分分析(PCA)

对指标体系进行降维处理
看到某O奖论文里用这个了～

步骤

  1. 对原始数据进行标准化处理
    \(x_{ij}^\ast=\frac{x_{ij}-\overline{x_j}}{\sqrt{Var(x_j)}}\) \(i=1,2,\cdot\cdot\cdot,n;j=1,2,\cdot\cdot\cdot,p)\)
    其中，\(\overline{x_j}=\frac{1}{n}\sum_{i=1}^nx_{ij}\)
    \(Var(x_j)=\frac{1}{n-1}\sum_{i=1}^n(x_{ij}-\overline{x_j})^2 (j=1,2,\cdot\cdot\cdot,p)\)
  2. 计算样本相关系数矩阵 corrcoef()
    \(r_{ij}=cov(x_i,x_j)\)
  3. 计算相关系数矩阵\(R\(的特征值 eig()
  4. 计算贡献率，选择主成分
    \(贡献率=\frac{\lambda_i}{\sum_{i=1}^p\lambda_i}\\0
  5. 计算主成分得分
    \(F_{ij}=\sum_{i=1}^pa_{jk}x_{ik} ,i=1,2,\cdot\cdot\cdot,n;j=1,2,\cdot\cdot\cdot,p\)
  6. 依据主成分得分数据，进一步对问题进行后续的分析和建模

Matlab

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
%数据导入
A = [];
%	标准化处理
a = size(A,1);
b = size(A,2);
for i = 1:b
	SA(:,i) = (A(:,i) - mean(A(:,i))) / std(A(:,i));
end
%计算相关系数矩阵
CM = corrcoef(SA);
[V,D] = eig(CM);
for j = 1:b
    DS(j,1) = D(b + 1 - j,b + 1 - j);
end
for i = 1:b
	DS(i,2) = DS(i,1) / sum(DS:,1));%贡献率
	DS(i,3) = sum(DS(1:i,1)) / sum(DS(:,1));%累积贡献率
end
%选择主成分
T = 0.9; %主成分信息保留率
for k = 1:b
	if DS(k,3) >= T
		Com_num = k;
		break;
	end
end
%提取主成分对应的特征向量
for j = 1:Com_num
	PV(:,j) = V(:,b + 1 - j);
end
%计算主成分得分
new_score = SA * PV;
for i = 1:a
	total_score(i,1) = sum(new_score(i,:));
	total_score(i,2) = i;
end
result_report = [new_score,total_score];
result_report = sortrows(result_report,-4);%将总分按降序排序

支持向量机(SVM)

一个小故事

也可以看看视频



Matlab

看不懂书上的代码，从网上扒了一个
线性的

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 function [ classification ] = SVM_L( train,test )
%  进行SVM线性可分的二分类处理
%  1、首先需要一组训练数据train，并且已知训练数据的类别属性，在这里，属性只有两类，并用1，2来表示。
%  2、通过svmtrain（只能处理2分类问题）函数，来进行分类器的训练
%  3、通过svmclassify函数，根据训练后获得的模型svm_struct，来对测试数据test进行分类
train=[0 0;2 4;3 3;3 4;4 2;4 4;4 3;5 3;6 2;7 1;2 9;3 8;4 6;4 7;5 6;5 8;6 6;7 4;8 4;10 10];                                              %训练数据点
group=[1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2]'; %训练数据已知分类情况
                                                                            %与train顺序对应
test=[3 2;4 8;6 5;7 6;2 5;5 2];                                                    %测试数据
%训练分类模型
svmModel = svmtrain(train,group,'kernel_function','linear','showplot',true);
%分类测试
classification=svmclassify(svmModel,test,'Showplot',true);
end

非线性的

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function [ classfication ] = SVM_NL( train,tast )
%SVM对线性不可分的数据进行处理
%在选择核函数时，尝试用linear以外的rbf,quadratic,polynomial等，观察获得的分类情况
%训练数据
train=[5 5;6 4;5 6;5 4;4 5;8 5;8 8;4 5;5 7;7 8;1 2;1 4;4 2;5 1.5;7 3;10 4;4 9;2 8;8 9;8 10]; 
%训练数据分类情况                                                                                      
group=[1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2];                          
%测试数据
test=[6 6;5.5 5.5;7 6;12 14;7 11;2 2;9 9;8 2;2 6;5 10;4 7;7 4];                 
%训练分类模型
svmModel = svmtrain(train,group,'kernel_function','rbf','showplot',true);
%分类
classification=svmclassify(svmModel,test,'Showplot',true);
end

有情景的

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
function [ classfication ] = SVM2_2( train,test )
%使用matlab自带的关于花的数据进行二分类实验（150*4），其中，每一行代表一朵花，
%共有150行（朵)，每一朵包含4个属性值（特征），即4列。且每1-50，51-100，101-150行的数据为同一类，分别为setosa青风藤类，versicolor云芝类，virginica锦葵类
%实验中为了使用svmtrain（只处理二分类问题）因此，将数据分为两类，51-100为一类，1-50和101-150共为一类
%实验先选用2个特征值，再选用全部四个特征值来进行训练模型，最后比较特征数不同的情况下分类精度的情况。
load fisheriris   %下载数据包含:meas（150*4花特征数据）
%和species（150*1 花的类属性数据）
meas=meas(:,1:2);  %选取出数据前100行，前2列
train=[(meas(51:90,:));(meas(101:140,:))];%选取数据中每类对应的前40个作为训练数据
test=[(meas(91:100,:));(meas(141:150,:))];%选取数据中每类对应的后10个作为测试数据
group=[(species(51:90));(species(101:140))];%选取类别标识前40个数据作为训练数据
%使用训练数据，进行SVM模型训练
svmModel = svmtrain(train,group,'kernel_function','rbf','showplot',true);
%使用测试数据，测试分类效果
classfication = svmclassify(svmModel,test,'showplot',true);
%正确的分类情况为groupTest，实验测试获得的分类情况为classfication
groupTest=[(species(91:100));(species(141:150))]; 
%计算分类精度
count=0;
for i=(1:20)
   if strcmp(classfication(i),groupTest(i))
      count=count+1;
   end
end
fprintf('分类精度为：%f\n' ,count/20);
end

K-均值（k-Means）

划分聚类的方法，就是变成k个堆

Matlab

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
%	这里假定每个样本具有两个特征，共有20个样本
x = [0 0;1 0;0 1;1 1;2 1;1 2];%我懒得打完了QAQ
z = zeros(2,2);
z1 = zeros(2,2);
z = x(1:2,1:2);
%寻找聚类中心
while 1
	count = zeros(2,1);
	allsum = zeros(2,2);
	for i = 1:20 % 对每一个样本i，计算到2个聚类中心的距离
		temp1 = sqrt((z(1,1) - x(i,1)) .^2 + (z(1,2) - x(i,2)) .^2);
		temp2 = sqrt((z(2,1) - x(i,1)) .^2 + (z(2,2) - x(i,2)) .^2);
		if(temp1<temp2)
			count(1) = count(1) + 1;
			allsum(1,1) = allsum(1,1) + x(i,1);
			allsum(1,2) = allsum(1,2) + x(i,2);
		else
			count(2) = count(2) + 1;
			allsum(2,1) = allsum(2,1) + x(i,1);
			allsum(2,2) = allsum(2,2) + x(i,2);
		end
	end
	z1(1,1) = allsum(1,1)/count(1);
	z1(1,2) = allsum(1,2)/count(1);
	z1(2,1) = allsum(2,1)/count(2);
	z1(2,2) = allsum(2,2)/count(2);
	if(z == z1)
		break;
	else
		z = z1;
	end
end
% 结果显示
disp(z1);%输出聚类中心
plot(x(:,1),x(:,2),'b*');
hold on
plot(z1(:,1),z1(:,2),'ro');
title('K均值法分类图');
xlabel('特征x1');
ylabel('特征x2');

额。。。网上的貌似直接有函数么。。。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
% 使用高斯分布（正态分布）
% 随机生成3个中心以及标准差
s = rng(5,'v5normal');
mu = round((rand(3,2)-0.5)*19)+1;
sigma = round(rand(3,2)*40)/10+1;
X = [mvnrnd(mu(1,:),sigma(1,:),200); ...
     mvnrnd(mu(2,:),sigma(2,:),300); ...
     mvnrnd(mu(3,:),sigma(3,:),400)];
%mvnrnd(mu,sigma,n);产生多维正态随机数，mu为期望向量，sigma为协方差矩阵，n为规模。
% 作图
P1 = figure;clf;
scatter(X(:,1),X(:,2),10,'ro');
title('研究样本散点分布图')
% 距离用传统欧式距离，分成两类
[cidx,cmeans,sumd,D] = kmeans(X,2,'dist','sqEuclidean');
%X N*P的数据矩阵
%cidx N*1的向量,存储的是每个点的聚类标号
%cmeans K*P的矩阵,存储的是K个聚类质心位置
%sumd 1*K的和向量,存储的是类间所有点与该类质心点距离之和
%D N*K的矩阵，存储的是每个点与所有质心的距离;
P2 = figure;
clf;%用来清除图形的命令
[silh,h] = silhouette(X,cidx,'sqeuclidean');

朴素贝叶斯判别模型

求后验概率时用到了Bayes公式

步骤

  1. 求解先验概率
  2. 计算条件概率
  3. 计算后验概率

Matlab

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
%导入数据
training_sample = [];
class = [];
test = [];
test_class = [];
%执行程序，注意不可有一元素完全相同的，这样就不能分类了，可以把这一列（行）删除
ObjBayes = NaiveBayes.fit(training_sample,class);
pred = ObjBayes.predict(test);
strcmpi(pred,test_class);
%计算分类准确率
A = sum(strcmpi(pred,test_class));
B = length(strcmpi(pred,test_class));
acc_tate = A/B
%分类错误的编号
find(strcmpi(pred,test_class) == 0)

应用和小废话

2012CUMCM A题
大数据的处理。。。
感觉这题很注重模型的使用，挑对了，不就对了吗？
感觉和美赛不一样
博博说A题有很多物理知识
我们应该不会选的
我是不是白白浪费了一天
唔～不开心
明天继续下一个专题喽

后记

最近想让自己静下来
看《瓦尔登湖》
做瑜伽
但是貌似书看的杂杂乱乱的
🧘‍♀️都没做完就去折腾💻了
唔～这样很不好的呀
说好的，
期末完以后要多动动，不能再总是坐着的
下雪了，而且操场还远，
计划泡汤
说好的，要翻译文章，
书不能借，现在手头只有《瓦尔登湖》是英文版的
这中文看的都累
唔～打算回家买书
顺便买📓，嘻嘻嘻
话说，学校怎么只有吉品轩开了呢
真心不好吃
腊八腊八，冻掉下巴？！！
之前倒是没听过
不过正好下雪了
算吃了腊八粥吧～
冷的…还不好吃

继续加油喽，
修身齐家治国平天下

转载请注明出处，谢谢。

愿 我是你的小太阳



买糖果去喽