Matlab案例分析

2021-07-18 数学建模 matlab应用 Comments

[目录]

找出10万以内的所有亲和数
编写猜数字游戏
钻井问题
BP网络预测
M/M/1排队问题
利用相关系数修复碎纸条
多元统计分析
主成分回归估计

一、找出10万以内的所有亲和数

1.1 题目介绍

亲和数是这样一对正整数a和b，使得a的所有真因子的和等于b，而b的所有真因子的和等于a。亲和数问题最早由毕达哥拉斯学派发现和研究的。他们在研究数字的规律的时候发现有以下的性点的两个数：

$1+2+4+71+142=220 \\ 1+2+4+5+10+11+20+22+44+55+110=284$

就是220的真因子是$1、2、4、5、10、11、20、22、44、55、110$，它们的和是$284；284$的真因子是$1、2、4、71、142$,其和恰好是$220$。这是最早发现的一对亲和数，也是最小的一对亲和数。

考虑到1是每个整数的因子，把除去整数本身之外的所有因子叫做这个数的“真因子”。如果两个整数，其中每一个数的真因子的和都恰好等于另一个数，那么这两个数，就构成一对“亲和数”。

1.2 程序

%寻找100000以内的亲和数对
clc
E=[];
for j=1:100000;   
    T=tuo1(j);
    M=tuo1(tuo1(j));
    if j==tuo1(tuo1(j)) && M~=T
        C=[M,T];
        E=[E;C];
    end
end
F=E(1:2:end,:)

function T= tuo1(a)%函数1，寻找数的所用因数
    A=[];
    for i=1:100000
        if (mod(a,i)==0)
            A=[A,i];
        end
    end
    T=sum(A)-a;
end

1.3 结果

二、编写猜数字游戏

2.1 题目介绍

编写程序实现猜数字游戏游戏说明：

计算机随机生成一个4位数（4个数字互不相同）然后玩家输入一个4位数，计算机提示（？A？B）:

A 的含义是数字位置正确的有多少个?

B 的意思是仅数字正确的有多少个?

如计算机随机产生的为1234，玩家输入1325，则返回1A2B （数字1位置正确，而2和3仅数字正确，位置不对），程序在4A或玩家猜测10次时结束。

2.2 程序

clear
clc
format 
disp('猜数字开始了，现在已经产生了一个四位随机数')
X=randperm(9);
a1=X(1);
b1=X(2);
c1=X(3);
d1=X(4);
t=10;

while i~=4 && t>0
    M=input('请输入一个四位数：')

    a2=fix(M/1000);
    b2=fix((M-1000*a2)/100);
    c2=fix(M/10)-100*a2-10*b2;
    d2=M-1000*a2-100*b2-10*c2;


    i=0;%数字位置全对
    j=0;%仅数字对

    if a1==a2 
        i=i+1;
    end
    if b1==b2 
        i=i+1;
    end
    if c1==c2 
        i=i+1;
    end
    if d1==d2 
        i=i+1;
    end

    if a1==b2 ||a1==c2 ||a1==d2
        j=j+1;
    end
    if b1==a2 ||b1==c2 ||b1==d2
        j=j+1;
    end
    if c1==a2 ||c1==b2 ||c1==d2
        j=j+1;
    end
    if d1==a2 ||d1==b2 ||d1==c2
        j=j+1;
    end
    t=t-1;
    disp([num2str(i),'A',num2str(j),'B'])
end

2.3 结果

三、钻井问题

3.1 题目介绍

试编写程序计算经过怎样的坐标平移旋转，可以使得这些点横纵坐标小数部分都小于0.05 的数目达到最大。

十五口井的坐标如下：

$(0.5,2)|(1.41,3.5)|(3,1.5)|(3.37,3.51)|(3.4,5.5)|(4.72,2)|(4.72,6.24)|(5.34,4.10)$

$(7.57,2.01)|(8.38,4.50)|(8.98,3.41)|(9.5,0.8)|(9.57,1.2)|(9.58,1.58)|(9.58,8.05)$

3.2 程序

clear
clc
format 
t=-1;
G=[];

for t=1:100000
    t=t+1
    T=tuo1(t);
    G=([G;T]);
end
G;
H=find(G(:,4)>=6)
J=G(H(1:end),:)


function T=tuo1(t)
b=-1;
N=[];
E=[];
%montecarlo随机取点
% i\j  :蒙特卡洛取得随机点，代替坐标原点
i=10*rand(1);
j=9*rand(1);

% X1/Y1：相对于原坐标井的位置  
X1=[0.5 1.41 3.00 3.37 3.4 4.72 4.72 5.43 7.57 8.38 8.98 9.50 9.57 9.58 9.58];
Y1=[2 3.5 1.5 3.51 5.5 2 6.24 4.1 2.01 4.5 3.41 0.8 1.2 1.58 8.05];

% X2/Y2:X1/Y1减去i/j后的相对位置
X2=X1-i;
Y2=Y1-j;

% a  :各点到相对原点的距离
% A  ：各点相对相对原点的角度
a=sqrt(X2.^2+Y2.^2);
A=atand(Y2./X2);

% B  :产生0-90的角度偏移
for B=0:90
    b=b+1;
    
    n=0;
    for m=1:15
        
        X3=a.*cosd(A-B);
        Y3=a.*sind(A-B);
        [X3;Y3];
        X=abs(X3-fix(X3));
        Y=abs(Y3-fix(Y3));

        if X(m)<0.05 && Y(m)<0.05
            n=n+1;
        end
        N=[N,n];
    end
    N=max(N);
    e=[i,j,B,N];
    E=[E;e];   
end
T=E;
end

3.3 结果

四、利用BP神经网络预测上海市2009年6月的入境旅游人数

4.1 题目介绍

说明：以2008年1月到2009年4 月的数据为输入，2008年2 月到2009年5月的数据为输出，训练模型，最终依靠训练出的模型预测2009年6月的旅游人数。

时间	人数	时间	人数
2008-1	517640	2008-10	627675
2008-2	528780	2008-11	599271
2008-3	576064	2008-12	419753
2008-4	575336	2009-1	414747
2008-5	556293	2009-2	468685
2008-6	503906	2009-3	510617
2008-7	510391	2009-4	549291
2008-8	483831	2009-5	490005
2008-9	504473

4.2 程序

clear
clc
format bank
% 读取、处理数据数据
A=[517640 528780 576064 575336 556293 503906 510391 483831 504473 627675 599271 419753 414747 468685 510617 549291 490005];
p1=[A(1:4);A(2:5);A(3:6);A(4:7);A(5:8);A(6:9);A(7:10);A(8:11);A(9:12);A(10:13);A(11:14);A(12:15);A(13:16)]'
t1=[A(5:17)]
[P,minp,maxp,T,mint,maxt]=premnmx(p1,t1);

%迭代学习
net=newff(minmax(P),[12,5,1],{'tansig','tansig','purelin'},'traingdm')
net.trainParam.epochs=1000000;
net.trainParam.goal=0.00001
net.trainParam.lr=0.01;
[net,tr]=train(net,P,T)

%反向预测
y=sim(net,P)
Y=postmnmx(y,mint,maxt)
t=1:13;
plot(t,t1,'r*',t,Y,'bo');
title('入境人数实际值与预测值对比')
legend('真实值','预测值')
B=[t;t1;Y]'

%展开预测2009年6月的数据
c=A(14:17);
c1=premnmx(c);
c2=c1';
C1=sim(net,c2);
C=postmnmx(C1,mint,maxt);
disp(['经过多次迭代，预测2009年6月份入境旅游人数为：',num2str(C)])

4.3 结果

五、 M/M/1 问题

5.1 题目介绍

随机模拟问题：考虑一个简单的排队系统$M/M/1$，具体理论可参看《算法大全》中排队论模型一章，设顾客到达服务窗口的规律服从参数为 $\lambda$ 的泊松分布（等价于两个顾客之间的间隔时间服从均值为 $\frac{1}{\lambda}$ 的指数分布，又称参数为 $\lambda$ 的负指数分布，密度函数为：

$f(x)= \begin{cases} \lambda e^{-\lambda x} & \text{x>0} \\ 0 & \text{else}\end{cases}$

其中，其数学期望$E(X)=\frac{1}{\lambda}$

（注：一般来说，到达服务窗口的顾客流是一个泊松流，即在$(0,t)$时间间隔内，到达服务窗口的顾客人数服从参数为 $\lambda t$ 的泊松分布，我们知道，该泊松分布的数学期望为 $\lambda t$ ；该泊松流等价于任意相邻两个顾客到达的间隔时间服从参数为 $\lambda$ 的负指数分布。比如，若 $\lambda=5$ ，表示单位时间内平均（数学期望）有 $5$ 个顾客到达窗口，等价于任意相邻两个顾客到达的间隔时间为$\frac{1}{5}$ ）。

每个顾客在窗口接受服务的时间服从参数为 $\mu$ 的负指数分布，试在以下两种情况下作随机模拟：

1）$\lambda=5$ 人/分钟，$\mu=6$ 人/分钟；

2） $\lambda=\mu=5$人/分钟。

（Matlab中指数分布随机数生成函数为$exprnd(mu)$，其中$mu$为数学期望值，可以模拟很长一段时间排队系统的工作状况，比如取 $T=5000$分钟，进而得到系统的以下参数的模拟值：

1）队长：指排队系统中等待服务的顾客人数的平均值。

2）逗留时间：指一个顾客在排队系统中的平均停留时间。

3）忙期：指从顾客到达空闲服务机构起，到服务机构再次为空闲止这段时间的平均长度。

请将模拟值与理论值比较一下。

5.2 程序

%% 计算模块
 clear

 tic
 N=0;
 Tim=0;
 J=1;
 t=5000;
 lamda=5;
 mu=6;

 for j=1:J
     [nn,tim]=tuo(j)
     Tim=Tim+tim;
     N=N+nn;
 end
 N=N/J;
 Tim=Tim/J;
 Jj=lamda/(mu*(mu-lamda))+1/mu;
 disp('当lamda=5,mu=6时：')
 disp(['每次模拟时间为5000min,经过',num2str(J),'次模拟，得到如下平均结果：'])
 disp(['平均服务的客户人数：Lq=',num2str(N)])
 disp(['客户平均逗留时间及平均逗留忙期：Ws & B=',num2str(Tim)])
 disp('******************分界线******************')
 disp('理论结果如下：')
 disp(['平均服务的客户人数：Lq`=',num2str((lamda)^2/(mu*(mu-lamda)))])
 disp(['客户的平均逗留时间及平均服务忙期：Ws` & B`=','num2str(Jj)'])
 disp('******************分界线******************')
 toc


 %% 函数模块
 function [nn,tim]=tuo(t)
     t=5000;
     lamda=5;
     mu=6;
     i=1;
     w=0;
     e0=0;c0=0;
     x(1)=exprnd(1/lamda);
     c(1)=c0+x(1);
     b(1)=c(1);
     while b(i)<=1
         y(i)=exprnd(1/lamda);
         e(i)=b(i)+y(i);
         wait(i)=b(i)+y(i);
         w=w+wait(i);
         i=i+1;
         x(i)=exprnd(1/lamda);
         c(i)=c(i-1)+x(i);
         b(i)=max(c(i),e(i-1));
     end
     n=i-1;
     nn=n/t;
     tim=w/n;
 end

5.3 结果

六、利用相关系数修复碎纸条

6.1 题目说明

将一张写满内容A4纸平均分成19份，随机打乱顺序后尝试拼接完整。

由于第I片的右侧和第I+1片的左侧存在相关性，故可以通过求取相关系数得出联系！

$\rho_{xy}=\frac{cov(X,Y)}{\sqrt{D(X)}\sqrt{D(Y)}}=\frac{E((X-EX)(Y-EY))}{\sqrt{D(X)\sqrt{D(Y)}}}=\frac{EXY-EX·EY}{\sqrt{D(X)}\sqrt{D(Y)}}$

$r=\frac{\sum_{i=1}^n(x_i-\bar x(y_i-\bar y))}{\sqrt{\sum_{i=1}^n(x_i-\bar x)^2}\sum_{i=1}^n(y_i-\bar y)^2}$

6.2 程序

clear
clc


%% 函数准备
a1=imread('000.bmp');
[m,n]=size(a1);
dirname = 'ImageChips';
files = dir(fullfile(dirname, '*.bmp'));
a=zeros(m,n,19);
pic=[];

%% 算第一张纸条
E=[1];
F=[1:19];
for j=1:19
    e=tuo(j);
    E=[E,e];
end
f1 = intersect(E,F);
f=setdiff(F,f1);
disp(['利用相关系数，得出',num2str(f),'张纸片左侧与其他纸片右侧关系最小'])

%% 调用函数，测算结果
c=0;
TT=[9];
j=9;
D=[];
for c=1:18
    T=tuo(j);
    TT=[TT,T];
    j=T;
    T=tuo(j);
    c=c+1;
end
disp('纸片按照现行顺序，排列如下所示：')
TT
for ii = 1:length(files)
filename = fullfile(dirname, files(TT(ii)).name);
a(:,:,TT(ii))=imread(filename);
pic=[pic,a(:,:,TT(ii))];
end
%double(pic);
figure
imshow(pic,[])

%% 排序函数
function T=tuo(j)
    a1=imread('000.bmp');
    [m,n]=size(a1);
    %imshow(a1);%显示图像
    %image(a1);

    dirname = 'ImageChips';
    files = dir(fullfile(dirname, '*.bmp'));
    a=zeros(m,n,19);
    pic=[];
    for ii = 1:length(files)
    filename = fullfile(dirname, files(ii).name);
    a(:,:,ii)=imread(filename);
    pic=[pic,a(:,:,ii)];
    end
    %double(pic);
    %figure
    %imshow(pic,[])
    a=double(a);
    b=a./255;
    r=[];
    m=1;
    RRr=[];
    for i=1:19
        x=b(:,72,j);
        y=b(:,1,i);
        r1=sum((x-mean(x)).*(y-mean(y)));
        r2=sqrt(sum((x-mean(x)).^2))*sqrt(sum((y-mean(y)).^2));
        r=[r;m,(r1/r2)];
        %rr=corrcoef(x,y);%或利用库函数来求样本相关系数
        m=m+1;
        R=max(r(:,2));
        Rj=r(r(:,2)==R);%已知向量某一列特殊值求向量位置
    end

        RRr=[RRr;Rj];
        T=RRr;
end

6.3 结果

6.4 程序分析

6.4.1 图片及其灰度

%% 函数准备
a1=imread('000.bmp');
[m,n]=size(a1);
%imshow(a1);%显示图像
%image(a1);
dirname = 'ImageChips';
files = dir(fullfile(dirname, '*.bmp'));
a=zeros(m,n,19);
pic=[];

imread函数，读取图片灰度，数组。
$[m,n]$是一个取值范围在0~255之间的二维数组，记录图片灰度。
imshow展示图片灰度数组。
image将灰度向量转变为图片
dir - 列出文件夹内容:(列出“dirname”这个文件夹内所有后缀为.bmp的文件)

6.4.2 对文件夹下的文件重命名

%% 对文件夹下的文件重命名
files = dir('F:\Matlab\*.mat');
len=length(files);
for i=1:len
    oldname=files(i).name;
    newname=strcat(num2str(i), 'new.mat');
    command = ['rename' 32 oldname 32 newname];
    status = dos(command);
    if status == 0
        disp([oldname, ' 已被重命名为 ', newname])
    else
        disp([oldname, ' 重命名失败!'])
    end
end

七、多元统计分析法

7.1 题目分析

下表是多个省市全年经济发展指标，试对其发展进行评价!

版幅限制我对列表进行了转置。

省份	安徽	河南	湖北	湖南	江西	山西
生产总值	4812.70	8815.09	6320.40	5612.26	3495.94	3000.00
第一产业	932.40	1647.48	1020.09	1155.85	711.70	253.40
第二产业	2169.80	4515.35	2994.67	2214.41	1595.74	1810.10
第三产业	1710.50	2652.26	2305.72	2242.00	1188.50	978.90
人均总值	7448.60	9071.80	10505.80	8379.40	8161.30	9024.10
粮食产量	2743.00	4260.00	2100.12	2810.26	1803.40	1062.00
工业增加	1736.00	3862.18	1664.73	1781.14	1110.70	1568.50
固定投资	1914.20	3099.00	2356.38	1981.29	1820.00	1458.50
消费总额	1503.10	2808.17	2667.48	2069.84	1059.90	884.80
进出口	72.10	66.10	67.72	54.38	35.30	53.80
预算收入	274.40	428.80	303.69	348.80	205.80	255.20
居民收入	7511.40	7704.90	8023.00	8617.48	7560.00	7902.90
农民收入	2499.30	2553.15	2897.00	2838.00	2953.00	2589.60

7.2 解答过程

Spss求解
- 1、对其中6(m)个变量13(z)个观测值,在分析-描述统计-描述中进行标准化处理，形成$[m\times z]$的矩阵A；
- 2、对标准化处理后的矩阵在分析-降维-因子中进行因子降维，用主成分分析法提取出特征值大于特定值的n个重要成分，形成$[z\times n]$的矩阵B；
- 3、将矩阵$A\times B$，即$[m,z]\times[z,n]=[m,n]$，的到成分系数与标准化参数的乘积矩阵SS;
- 4、对第二步中的选取的n个成分系数的权重重新分配，得到$[n,1]$的权重矩阵L；
- 5、矩阵$L\times S$，即得到评价值。

6-14
6-15

Matlab求解

A=xlsread('D:\File\文档\MATLAB\暑期训练\多元统计分析\各省GDP占比权重排名\a.xlsx','Sheet1','B1:N7')%去取数据
[B,B_mean,B_std]=zscore(A)%对A进行标准化处理

[coeff,score,latent] = pca(B);
latent%特征值
l=latent(latent>1) %选取特征值大于一的
L=l/sum(l)%权重矩阵
score %成分系数与标准化参数的乘积矩阵 
S=score(:,1:3)
E=S(:,1)*L(1)+S(:,2)*L(2)+S(:,3)*L(3)
G=["安徽";"河南";"湖北";"湖南";"江西";"山西"]

disp('各省评价得分如下所示')
H=["省份","得分";G,E]

6.3 结果

6-12

6-13

八、主成分回归估计

8.1 题目分析

主成分估计（principal component estimate）是Massy在1965年提出的，它是回归系数参数的一种线性有偏估计（biased estimate），同其它有偏估计，如岭估计（ridge estimate）等一样，是为了克服最小二乘（LS）估计在设计阵病态（即存在多重共线性）时表现出的不稳定性而提出的。主成分估计采用的方法是将原来的回归自变量变换到另一组变量，即主成分，选择其中一部分重要的主成分作为新的自变量（此时丢弃了一部分，影响不大的自变量，这实际达到了降维的目的），然后用最小二乘法对选取主成分后的模型参数进行估计，最后再变换回原来的模型求出参数的估计。

求解水泥的影响因子与热量值之间的线性关系：

影响因子1	影响因子2	影响因子3	影响因子4	热量值
7	26	6	60	78.5
1	29	15	52	74.3
11	56	8	20	104.3
11	31	8	47	87.6
7	52	6	33	95.9
11	55	9	22	109.2
3	71	17	6	102.7
1	31	22	44	72.5
2	54	18	22	93.1
21	47	4	26	115.9
1	40	23	34	83.8
11	66	9	12	113.3
10	68	8	12	109.4

8.2 程序解答

clear
clc

A=xlsread('D:\File\文档\MATLAB\暑期训练\多元统计分析\对水泥性质影响因素的线性回归\水泥.xlsx','Sheet1','A2:D14');
B=xlsread('D:\File\文档\MATLAB\暑期训练\多元统计分析\对水泥性质影响因素的线性回归\水泥.xlsx','Sheet2','A2:A14');
[C,C_mean,C_std]=zscore(A)%对A标准化处理：C为标准化处理后的值；C_mean为均值；C_std为标准方差
[G,S,L]=pca(C);
L=L(1:3)%取特征值的前三个

%求特征向量
d1=size(C,1);
d=(C'*C)/(d1-1);
D=cov(d);
e=pcacov(D)';%特征向量
E=e(1:3,:)'
F=C*E %F=S

Q=E'./C_std %

syms X1 X2 X3 X4 F1 F2 F3
F1=Q(1,1)*(X1-C_mean(1))+Q(1,2)*(X2-C_mean(2))+Q(1,3)*(X3-C_mean(3))+Q(1,4)*(X4-C_mean(4));
F2=Q(2,1)*(X1-C_mean(1))+Q(2,2)*(X2-C_mean(2))+Q(2,3)*(X3-C_mean(3))+Q(2,4)*(X4-C_mean(4));
F3=Q(3,1)*(X1-C_mean(1))+Q(3,2)*(X2-C_mean(2))+Q(3,3)*(X3-C_mean(3))+Q(3,4)*(X4-C_mean(4));

y=95.423+9.883*F1+0.125*F2+4.555*F3;

Y=vpa(y,4)

8.3 结果

6-16

6-17

本文链接： http://tuoxiao.github.io/2021/07/18/Mathematical_modeling/matlab6/

TuTo

天地与我并生，万物与我为一。

Matlab案例分析

[目录]

一、 找出10万以内的所有亲和数

1.1 题目介绍

1.2 程序

1.3 结果

二、编写猜数字游戏

2.1 题目介绍

2.2 程序

2.3 结果

三、钻井问题

3.1 题目介绍

3.2 程序

3.3 结果

四、利用BP神经网络预测上海市2009年6月的入境旅游人数

4.1 题目介绍

4.2 程序

4.3 结果

五、 M/M/1 问题

5.1 题目介绍

5.2 程序

5.3 结果

六、利用相关系数修复碎纸条

6.1 题目说明

6.2 程序

6.3 结果

6.4 程序分析

6.4.1 图片及其灰度

6.4.2 对文件夹下的文件重命名

七、多元统计分析法

7.1 题目分析

7.2 解答过程

6.3 结果

八、主成分回归估计

8.1 题目分析

8.2 程序解答

8.3 结果

TuTo

一、找出10万以内的所有亲和数