Преобразуйте R-код в MATLAB

Question

r код для алгоритма с фиктивными данными: Я работаю, чтобы перевести это на MATLAB, но борясь с вычислениями, которые выполняются внутри цикла. Любая помощь будет оценена.

data <- c(-0.39, 0.12, 0.94, 1.67, 1.76, 2.44, 3.72, 
     4.28, 4.92, 5.53, 0.06, 0.48, 1.01, 1.68, 1.80, 
     3.25, 4.12, 4.60, 5.28, 6.22) 
pi <- 0.546 
sigmas1 <- 0.87 
sigmas2 <- 0.77 
mu1 <- numeric(0) 
mu2 <- numeric(0) 
r <- numeric(0) 
R1 <- matrix (0 ,20 ,100) 
mu1[1] <- 4.62 
mu2[1] <- 1.06 

for(j in 1:100){ 
    for (i in 1:20){ 
    r [i] <- pi * dnorm (data[i] , mu2[j], sigmas2^(1/2))/((1- pi)*dnorm(data[i], 
     mu1[j], sigmas1^(1/2))+ pi*dnorm(data[i], mu2[j], sigmas2^(1/2))) 
R1[i, j] <- r[i] 
} 
r 
mu1[j+1] <- sum((1-r)*data)/sum(1-r) 
mu2[j+1] <- sum(r*data)/sum(r) 
Muu1 <- mu1[j+1] 
Muu2 <- mu2[j+1] 
} 

Muu1 
Muu2 
x11() 
layout(matrix(c(1, 2))) 
plot(mu1, type="l", main="", xlab="EM Iteration for the Fictitious Data") 
plot(mu2, type="l", main="", xlab='EM Iteration for the Fictitious Data') 

Answer 1

MATLAB эквивалент dnorm функции R является normpdf. Аргументы те же, что и в R:

normpdf(X,mu,sigma) 

С этим можно легко адаптировать петлю. Поскольку функция normpdf позволяет векторам в качестве входов, вы можете сбросить внутренний цикл и вместо этого использовать векторизованный подход. Всегда имейте в виду, что * и / являются матричным умножением и делением в MATLAB. Чтобы получить элементарные операторы, используйте вместо этого .* и ./.

Обратите внимание, что в MATLAB лучше предусмотреть все переменные. Так как mu1 и mu2 идут от 1 до 100, но на каждом шаге вы устанавливаете значение mu[j+1], оно будет иметь размер 1x101. Для r и R1 размер ясный думаю.

Все вместе, это дало бы следующий код:

data = [-0.39, 0.12, 0.94, 1.67, 1.76, 2.44, 3.72,... 
    4.28, 4.92, 5.53, 0.06, 0.48, 1.01, 1.68, 1.80,... 
    3.25, 4.12, 4.60, 5.28, 6.22]; 

pi=0.546; 
sigmas1 = 0.87; 
sigmas2 = 0.77; 
mu1 = zeros(1,101); 
mu2 = zeros(1,101); 
r = zeros(1,20); 
R1 = zeros(20,100); 
mu1(1) = 4.62; 
mu2(1) = 1.06; 

for j=1:100 
    r= pi*normpdf(data,mu2(j),sigmas2^(1/2)) ./ ... 
     ((1-pi)*normpdf(data,mu1(j),sigmas1^(1/2)) + ... 
     pi*normpdf(data,mu2(j),sigmas2^(1/2))); 
    R1(:,j) = r; 
    mu1(j+1) = sum((1-r).*data)/sum(1-r); 
    mu2(j+1) = sum(r.*data)/sum(r); 
end 

figure; 
subplot(1,2,1); 
plot(mu1); 
subplot(1,2,2); 
plot(mu2); 

Если это не работает для вас, или у Вас есть какие-либо вопросы по коду, не стесняйтесь комментировать.