用R语言优化并更改以下代码的变量名称set.seed(123) n <- 1000 mu1 <- c(0,4) mu2 <- c(-2,0) Sigma1 <- matrix(c(3,0,0,0.5),nr=2,nc=2) Sigma2 <- matrix(c(1,0,0,2),nr=2,nc=2) phi <- c(0.6,0.4) X <- matrix(0,nr=2,nc=n) for (i in 1:n) { if (runif(1)<=phi[1]) { X[,i] <- mvrnorm(1,mu=mu1,Sigma=Sigma1) }else{ X[,i] <- mvrnorm(1,mu=mu2,Sigma=Sigma2) } } ##initial guess for parameters mu10 <- runif(2) mu20 <- runif(2) Sigma10 <- diag(2) Sigma20 <- diag(2) phi0 <- runif(2) phi0 <- phi0/sum(phi0) ##EM algorithm k=2 prob <- matrix(rep(0,k*n),ncol = 2) weight <- matrix(rep(0,k*n),ncol = 2) phi <- phi0 mu <- matrix(c(mu10,mu20),nr=2) Sigma <- matrix(c(Sigma10,Sigma20),nr=2) #for loop,set up max iterations (200) for (step in 1:200) { for (j in 1:k) { for (i in 1:1000) { prob[i,j] <- dmvnorm(X[,i], mu[,j], Sigma[,(2*j-1):(2*j)]) weight[i,j] <- phi[j] * prob[i,j] } } row_sum <- rowSums(weight) prob <- weight/row_sum #prob denotes "wij" hear # note the parameters of the last iteration oldphi <- phi oldmu <- mu oldSigma <- Sigma # M-step:calculate the next theta by maximizing g(theta) for (j in 1:k) { sum1 <- sum(prob[, j]) sum2 <- X%*%prob[, j] phi[j] <- sum1/n mu[,j] <- sum2/sum1 sum3 <- matrix(c(0,0,0,0),nr=2) for (m in 1:n) { sum30 <- ((X[,m]-mu[,j])%*%t(X[,m]-mu[,j]))*prob[m,j] sum3 <- sum3+sum30 } Sigma[,(2*j-1):(2*j)] <- sum3/sum1 } # Set threshold: convergence is considered when the parameter obtained from the previous iteration has little change from the parameter obtained from the next iteration threshold <- 1e-5 if (sum(abs(phi - oldphi)) < threshold & sum(abs(mu - oldmu)) < threshold & sum(abs(Sigma - oldSigma)) < threshold) break #print the parameters in every iteration cat('step', step, 'phi', phi, 'mu', mu, 'Sigma', Sigma, '\n') }
时间: 2023-06-04 17:06:20 浏览: 771
以下是对该代码的变量名称和优化的建议:
set.seed(123) # 设置种子,保证每次运行结果相同
n <- 1000 # 样本数
mu1 <- c(0, 4) # 正态分布1的均值向量
mu2 <- c(-2, 0) # 正态分布2的均值向量
Sigma1 <- matrix(c(3, 0, 0, 0.5), nr = 2, nc = 2) # 正态分布1的协方差矩阵
Sigma2 <- matrix(c(1, 0, 0, 2), nr = 2, nc = 2) # 正态分布2的协方差矩阵
phi <- c(0.6, 0.4) # 混合比例
X <- matrix(0, nr = 2, nc = n) # 初始化数据矩阵
# 用循环生成样本数据
for (i in 1:n) {
if (runif(1) <= phi[1]) { # 以phi[1]为概率生成正态分布1的样本
X[, i] <- mvrnorm(1, mu = mu1, Sigma = Sigma1)
} else { # 以phi[2]为概率生成正态分布2的样本
X[, i] <- mvrnorm(1, mu = mu2, Sigma = Sigma2)
}
}
变量名称优化建议:
1. 使用有意义的变量名,避免使用单个字母的变量名。
2. 变量名最好全部使用小写字母。
3. 尽量缩短变量名,但要保证可读性。
例如:
set.seed(123) # 设置随机数种子
n <- 1000 # 样本数
mean1 <- c(0, 4) # 正态分布1的均值向量
mean2 <- c(-2, 0) # 正态分布2的均值向量
cov1 <- matrix(c(3, 0, 0, 0.5), nr = 2, nc = 2) # 正态分布1的协方差矩阵
cov2 <- matrix(c(1, 0, 0, 2), nr = 2, nc = 2) # 正态分布2的协方差矩阵
ratio <- c(0.6, 0.4) # 混合比例
data <- matrix(0, nr = 2, nc = n) # 初始化数据矩阵
阅读全文
相关推荐


















