我正在尝试创建(在 r 中)与以下 MATLAB 函数等效的函数,该函数将从 N(m1,(s1)^2) 和 N(m2, (s2)^2) 的混合中生成 n 个样本带有来自第一个高斯分布的分数 alpha。
我有一个开始,但 MATLAB 和 R 之间的结果明显不同(即,MATLAB 结果偶尔给出 +-8 的值,但 R 版本甚至从来没有给出 +-5 的值)。 请帮我解决这里的问题。谢谢:-)
例如: 从 N(0,1) 和 N(0,36) 的混合中绘制 1000 个样本,其中 95% 的样本来自第一个高斯分布。将样本归一化为均值零和标准差一。
MATLAB
函数
function y = gaussmix(n,m1,m2,s1,s2,alpha)
y = zeros(n,1);
U = rand(n,1);
I = (U < alpha)
y = I.*(randn(n,1)*s1+m1) + (1-I).*(randn(n,1)*s2 + m2);
实现
P = gaussmix(1000,0,0,1,6,.95)
P = (P-mean(P))/std(P)
plot(P)
axis([0 1000 -15 15])
hist(P)
axis([-15 15 0 1000])
结果图
结果历史记录
R
yn <- rbinom(1000, 1, .95)
s <- rnorm(1000, 0 + 0*yn, 1 + 36*yn)
sn <- (s-mean(s))/sd(s)
plot(sn, xlim=range(0,1000), ylim=range(-15,15))
hist(sn, xlim=range(-15,15), ylim=range(0,1000))
结果图
结果历史记录
一如既往,谢谢!
解决方案
gaussmix <- function(nsim,mean_1,mean_2,std_1,std_2,alpha){
U <- runif(nsim)
I <- as.numeric(U<alpha)
y <- I*rnorm(nsim,mean=mean_1,sd=std_1)+
(1-I)*rnorm(nsim,mean=mean_2,sd=std_2)
return(y)
}
z1 <- gaussmix(1000,0,0,1,6,0.95)
z1_standardized <- (z1-mean(z1))/sqrt(var(z1))
z2 <- gaussmix(1000,0,3,1,1,0.80)
z2_standardized <- (z2-mean(z2))/sqrt(var(z2))
z3 <- rlnorm(1000)
z3_standardized <- (z3-mean(z3))/sqrt(var(z3))
par(mfrow=c(2,3))
hist(z1_standardized,xlim=c(-10,10),ylim=c(0,500),
main="Histogram of 95% of N(0,1) and 5% of N(0,36)",
col="blue",xlab=" ")
hist(z2_standardized,xlim=c(-10,10),ylim=c(0,500),
main="Histogram of 80% of N(0,1) and 10% of N(3,1)",
col="blue",xlab=" ")
hist(z3_standardized,xlim=c(-10,10),ylim=c(0,500),
main="Histogram of samples of LN(0,1)",col="blue",xlab=" ")
##
plot(z1_standardized,type='l',
main="1000 samples from a mixture N(0,1) and N(0,36)",
col="blue",xlab="Samples",ylab="Mean",ylim=c(-10,10))
plot(z2_standardized,type='l',
main="1000 samples from a mixture N(0,1) and N(3,1)",
col="blue",xlab="Samples",ylab="Mean",ylim=c(-10,10))
plot(z3_standardized,type='l',
main="1000 samples from LN(0,1)",
col="blue",xlab="Samples",ylab="Mean",ylim=c(-10,10))
最佳答案
我认为有两个问题... (1) 您的 R 代码正在创建正态分布的混合,标准差为 1 和 37。 (2) 通过设置prob
等于你的 alpha rbinom()
调用,您将在 second 模式而不是第一个模式中获得分数 alpha。所以你得到的是一个分布,它主要是一个 sd 37 的高斯分布,被 5% 的高斯和 sd 1 的混合污染,而不是一个 sd 1 的高斯分布,它被 5% 的高斯和 sd 6 的混合污染. 按混合的标准偏差(大约 36.6)缩放基本上将其降低为标准高斯分布,在原点附近有轻微的凸起 ...
(此处发布的其他答案确实很好地解决了您的问题,但我认为您可能对诊断感兴趣...)
您的 Matlab 的更紧凑(也许更惯用)版本 gaussmix
函数(我认为 runif(n)<alpha
比 rbinom(n,size=1,prob=alpha)
稍微更有效率)
gaussmix <- function(n,m1,m2,s1,s2,alpha) {
I <- runif(n)<alpha
rnorm(n,mean=ifelse(I,m1,m2),sd=ifelse(I,s1,s2))
}
set.seed(1001)
s <- gaussmix(1000,0,0,1,6,0.95)
关于r - 从 r 中的双高斯混合生成样本(MATLAB 中给出的代码),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/12450007/