Блог

Главная - Вопросы по программированию - Преобразование R-кода в код MATLAB: застрял в sapply ()

Преобразование R-кода в код MATLAB: застрял в sapply ()

#r #matlab #debugging

#r #matlab #отладка

Вопрос:

У меня есть следующий 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)

Однако я не получаю желаемых результатов. Кривые должны перекрываться друг с другом: синяя кривая правильная, поэтому оранжевая должна перекрываться с синей кривой.

введите описание изображения здесь

Чего мне здесь не хватает? pnorm() Эквивалентен ли R коду MATLAB normcdf() , как я сделал здесь?

[ПРАВКА 2: НАЙДЕНА ОШИБКА!]

Итак, повозившись, я обнаружил, что все, что мне нужно было сделать, это получить количество вхождений tb < pb . Строка summat = sumidx.*repmat(ps,300,1) не должна быть там. После удаления этой строки и сохранения sumn = sum(sumidx,2); я получаю желаемый результат.

введите описание изображения здесь

Комментарии:

1. sapply применяет функцию, определенную как function(tb){sum(ps < tb)} , к списку или вектору ts и возвращает список. Это похоже на forloop, применяющий функцию к каждому элементу в списке. Следовательно tb , существует только как переменный заполнитель в применяемой функции.

2. Таким образом, функция принимает одно значение tb и идентифицирует значения, ps которые меньше этого конкретного tb значения, и суммирует их. Затем повторите то же самое для остальных значений tb ? Или сначала он выполняет поиск всех tb значений и суммирует их в конце цикла?

Ответ №1:

Итак, основываясь на комментариях от @Croote и повозившись, я придумал следующий код для функции, определенной в sapply() 

 sumidx  = bsxfun(@lt,ps,ts');
sumn    = sum(sumidx,2);

И для графика я закодировал его как 

 sumnfin = sumn(2:end)-sumn(1:end-1);
plot(t(1:length(sumn)-1),sumnfin)

Наконец, я получаю желаемый результат

введите описание изображения здесь