博文
随机抽样一致性算法matlab示例代码
|
2017年05月05日 21:05:05
前一篇博文曾提到随机抽样一致性算法, 这里给出一段利用这种方法计算扩散系数的matlab示例代码. 代码中同时测试了matlab的全局优化算法.
| msd.m | |
|---|---|
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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87 | functionMSD
clc;clearall;closeall;globaltmsdNincRepsIzro%%文件名称及Ransac设定
file='msd.xvg';
Npnt=5;%随机选择点数
Iter=1000;%最大循环数
Reps=1.5;%最大偏差
Rinc=0.5;%使用点数的比例
Izro=0;%拟合直线是否过零点,0:过;非零:不过%%获取文件头行数
fid=fopen(file);
txt=fgetl(fid);
idx=length(strfind(txt,'#'))+length(strfind(txt,'@'));
Line=0;whileidx>0Line=Line+1;txt=fgetl(fid);idx=length(strfind(txt,'#'))+length(strfind(txt,'@'));end%%读入数据,准备作图
[t,msd]=textread(file,'%f%f','headerlines',Line);
h=plot(0,[0;0],'EraseMode','background');
set(h(1),'XData',t,'YData',msd);
axis([0max(t)0max(msd)]);%%RANSAC
Ndat=length(msd);
Ninc=round(Rinc*Ndat);
Nbst=0;Pbst=[];
Pcst=zeros(Npnt,1);ifIzro~=0;Pcst=ones(Npnt,1);end
warning('off','MATLAB:rankDeficientMatrix');fori=1:Iteridx=randperm(Ninc,Npnt);%随机选点t0=[t(idx),Pcst];par=(t0msd(idx))';%或使用polyfit(x,y,n),通用但速度慢dst=abs(par(1)*t(1:Ninc)+par(2)-msd(1:Ninc));%点间距离%dst=abs(par(1)*t(1:Ninc)+par(2)-msd(1:Ninc))/sqrt(P(1)*P(1)+1);%垂直距离num=length(find(dst<Reps));ifmod(i,100)==0title(sprintf('Iter:%dNbst:%d',i,Nbst));drawnow;endifnum>NbstNbst=num;Pbst=par;set(h(2),'XData',t,'YData',Pbst(1)*t+Pbst(2));drawnow;endend%%全局搜索方法
opt=optimset('Algorithm','sqp');%interior-point/sqp/active-set
fmin=createOptimProblem('fmincon','objective',@ObjFun,...'x0',[0,0],'lb',[0,0],'ub',[inf,inf],'options',opt);
gs=GlobalSearch('Display','iter');
[Pgs,Ngs]=run(gs,fmin);ifIzro==0;Pgs(2)=0;end%%最小二乘方法
t0=[t(1:Ninc),zeros(Ninc,1)];ifIzro~=0;t0=[t(1:Ninc),ones(Ninc,1)];end
Plsq=(t0msd(1:Ninc))';%%作图
plot(t,msd,'k.',...t,polyval(Plsq,t),'r',...t,polyval(Pbst,t),'g',...t,polyval(Pgs,t),'b','linewidth',2)
legend('MSD','LSQ','RANSAC','GlobalSearch');
axis([0max(t)0max(msd)]);end%%functionObjFun=ObjFun(P)globaltmsdNincRepsIzroifIzro==0;P(2)=0;enddst=abs(P(1)*t(1:Ninc)+P(2)-msd(1:Ninc));%点间距离%dst=abs(P(1)*t(1:Ninc)+P(2)-msd(1:Ninc))/sqrt(P(1)*P(1)+1);%垂直距离ObjFun=length(find(dst>Reps));end |
◆本文地址: http://jerkwin.github.io/2017/05/05/随机抽样一致性算法matlab示例代码/, 转载请注明◆
https://blog.sciencenet.cn/blog-548663-1053261.html
上一篇：扩散模式的分类以及扩散系数的计算
下一篇：xvg曲线中心移动平均平滑脚本
