- 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 | 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
|
