У меня есть следующий код R, который я пытаюсь преобразовать в MATLAB.(Нет, я не хочу запускать код R в MATLAB, как показано здесь ).
Код R здесь:
# model parameters
dt <- 0.001
t <- seq(dt,0.3,dt)
n=700*1000
D = 1
d = 0.5
# model
ft <- n*d/sqrt(2*D*t^3)*dnorm(d/sqrt(2*D*t),0,1)
fmids <- n*d/sqrt(2*D*(t+dt/2)^3)*dnorm(d/sqrt(2*D*(t+dt/2)),0,1)
plot(t,ft*dt,type="l",lwd=1.5,lty=2)
# simulation
#
# simulation by drawing from uniform distribution
# and converting to time by using quantile function of normal distribution
ps <- runif(n,0,1)
ts <- 2*pnorm(-d/sqrt(2*D*t))
sumn <- sapply(ts, FUN = function(tb) sum(ps < tb))
lines(t[-length(sumn)],sumn[-1]-sumn[-length(sumn)],col=4)
И код MATLAB, который я сделал до сих пор:
% # model
ft = (n*d)./sqrt(2*D.*t.^3).*normpdf(d./sqrt(2*D.*t),0,1);
fmids = (n*d)./sqrt(2*D*((t+dt)./2).^3).*normpdf(d./sqrt(2*D.*((t+dt)./2)),0,1);
figure;plot(t,ft.*dt);
% # simulation
% #
% # simulation by drawing from uniform distribution
% # and converting to time by using quantile function of normal distribution
ps = rand(1,n);
ts = 2*normcdf(-d./sqrt(2*D*t));
Итак, вот где я застрял.Я не понимаю, что делает функция sumn = sapply(ts, FUN = function(tb) sum(ps < tb))
и откуда появился параметр 'tb
'.Он также не определен в данном R-коде.
Может кто-нибудь сказать мне, что эквивалент этого кода функции R в MATLAB?
[РЕДАКТИРОВАТЬ 1: ОБНОВЛЕНИЕ]
Итак, на основеВ комментариях @Croote я получил следующий код для функции, определенной в sapply()
sumidx = bsxfun(@lt,ps,ts');
summat = sumidx.*repmat(ps,300,1);
sumn = sum(summat,2);
sumnfin = sumn(2:end)-sumn(1:end-1);
plot(t(1:length(sumn)-1),sumnfin)
Однако я не получаю желаемых результатов.Кривые должны перекрываться друг с другом: синяя кривая правильна, поэтому оранжевый нужно перекрывать с синей кривой .
Что мне здесь не хватает?R pnorm()
эквивалентно normcdf()
MATLAB, как я сделал здесь?
[РЕДАКТИРОВАТЬ 2: НАЙТИ ЖУКУ!]
Итак, поигравшись, я обнаружил, что все, что мне нужно было сделать, это получить количество вхождений tb < pb
.Линия summat = sumidx.*repmat(ps,300,1)
не должна быть там.После удаления этой строки и сохранения sumn = sum(sumidx,2);
я получаю желаемый результат.