5

N(m1,(s1)^2) と N(m2, (s2)^2) の混合物から n 個のサンプルを生成する次の MATLAB 関数と同等のものを (r で) 作成しようとしています。 、アルファ、最初のガウスから。

まず始めに、MATLAB と R の間で結果が著しく異なります (つまり、MATLAB の結果では時折 +-8 の値が返されますが、R バージョンでは +-5 の値が返されることはありません)。ここで何が問題なのかを整理するのを手伝ってください。ありがとう ：-）

例: N(0,1) と N(0,36) の混合から 1000 個のサンプルをプロットし、最初のガウス分布からのサンプルの 95% を使用します。ゼロと標準偏差 1 を意味するようにサンプルを正規化します。

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])

結果のプロット

MATLAB の 2 つのガウス分布からランダムに生成されたサンプルのプロット

結果の履歴

MATLAB の 2 つのガウス分布からランダムに生成されたサンプルのヒストグラム

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))

結果のプロット

R の 2 つのガウス分布からランダムに生成されたサンプルのプロット

結果の履歴

R の 2 つのガウス分布からランダムに生成されたサンプルのヒストグラム

いつもありがとうございます！

解決

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))
4

4 に答える 4

6

2つの問題があると思います...（1）Rコードは、標準偏差が 1 と 37の正規分布の混合を作成しています。(2)呼び出しでprobalpha に等しく設定することにより、最初のモードではなく2 番目rbinom()のモードで分数アルファを取得します。したがって、sd 6 のガウスの 5% 混合物によって汚染された sd 1 のガウスではなく、sd 1 のガウスの 5% 混合物によって汚染された、ほとんどが sd 37 のガウスである分布が得られます。 . 混合物の標準偏差 (約 36.6) でスケーリングすると、基本的に、原点近くでわずかに隆起した標準ガウス分布に縮小されます ...

（ここに投稿された他の回答は問題を完全に解決しますが、診断に興味があるかもしれないと思いました...）

Matlab 関数のよりコンパクトな (そしておそらくより慣用的な) バージョン(よりもわずかに効率的gaussmixだと思います)runif(n)<alpharbinom(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)
于 2012-09-16T19:55:43.340 に答える
1

このタスクを実行するコードは次のとおりです。

「例: N(0,1) と N(0,36) の混合物から 1000 個のサンプルをプロットし、最初のガウス分布からのサンプルの 95% を使用します。サンプルを正規化して平均ゼロと標準偏差 1 にします。」

 plot(multG <- c( rnorm(950), rnorm(50, 0, 36))[sample(1000)] , type="h")
 scmulG <- scale(multG)
 summary(scmulG)
 #-----------    
   V1          
 Min.   :-9.01845  
 1st Qu.:-0.06544  
 Median : 0.03841  
 Mean   : 0.00000  
 3rd Qu.: 0.13940  
 Max.   :12.33107

ここに画像の説明を入力

于 2012-09-16T20:39:55.017 に答える
1

私は最近、正規分布の多項混合の密度とサンプリング関数を書きました。

dmultiNorm <- function(x,means,sds,weights)
{
  if (length(means)!=length(sds)) stop("Length of means must be equal to length of standard deviations")
  N <- length(x)
  n <- length(means)
  if (missing(weights))
  {
    weights <- rep(1,n)  
  }
  if (length(weights)!=n) stop ("Length of weights not equal to length of means and sds")
  weights <- weights/sum(weights)
  dens <- numeric(N)
  for (i in 1:n)
  {
    dens <- dens + weights[i] * dnorm(x,means[i],sds[i])
  }
  return(dens)
}

rmultiNorm <- function(N,means,sds,weights,scale=TRUE)
{
  if (length(means)!=length(sds)) stop("Length of means must be equal to length of standard deviations")
  n <- length(means)
  if (missing(weights))
  {
    weights <- rep(1,n)  
  }
  if (length(weights)!=n) stop ("Length of weights not equal to length of means and sds")

  Res <- numeric(N)
  for (i in 1:N)
  {
    s <- sample(1:n,1,prob=weights)
    Res[i] <- rnorm(1,means[s],sds[s])  
  }
  return(Res)
}

means平均のベクトルでsdsあり、標準偏差のベクトルでありweights、各分布からサンプリングする比例確率を持つベクトルです。これは役に立ちますか？

于 2012-09-16T19:48:54.037 に答える