随机抽样一致性算法matlab示例代码

类别:    标签: gmx matlab   阅读次数:   版权: (CC) BY-NC-SA

前一篇博文曾提到随机抽样一致性算法, 这里给出一段利用这种方法计算扩散系数的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
function MSD
clc; clear all; close all;
global t msd Ninc Reps Izro

%% 文件名称及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;
while idx>0
  Line= 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([0 max(t) 0 max(msd)]);

%% RANSAC
Ndat = length(msd);
Ninc = round(Rinc*Ndat);
Nbst = 0; Pbst=[];
Pcst = zeros(Npnt,1);
if Izro~=0; Pcst=ones(Npnt,1); end
warning('off', 'MATLAB:rankDeficientMatrix');
for i = 1:Iter
  idx = randperm(Ninc, Npnt); % 随机选点
  t0  = [t(idx),  Pcst];
  par = (t0\msd(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));
  if mod(i,100)==0
    title(sprintf('Iter: %d    Nbst: %d', i, Nbst));
    drawnow;
  end
  if num>Nbst
    Nbst = num; Pbst = par;
    set(h(2), 'XData',t, 'YData',Pbst(1)*t+Pbst(2));
    drawnow;
  end
end

%% 全局搜索方法
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);
if Izro==0; Pgs(2)=0; end

%% 最小二乘方法
t0  = [t(1:Ninc), zeros(Ninc,1)];
if Izro~=0; t0 = [t(1:Ninc), ones(Ninc,1)]; end
Plsq = (t0\msd(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([0 max(t) 0 max(msd)]);
end

%%
function ObjFun = ObjFun(P)
  global t msd Ninc Reps Izro
  if Izro==0; P(2)=0; end
  dst = 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

参考资料

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 扩散模式的分类以及扩散系数的计算
后一篇: xvg曲线中心移动平均平滑脚本

访问人次(2017年1月27日起): | 最后更新: 2017-12-13 23:49:42 UTC | 版权所有 © 2008 - 2017 Jerkwin