4

私は小さなサンプルサイズのデータ​​を扱っています:

>dput(dat.demand2050.unique)  
c(79, 56, 69, 61, 53, 73, 72, 86, 75, 68, 74.2, 80, 65.6, 60, 54)

密度分布は次のようになります。
データのpdf

値が 2 つの領域 (低と高) からのものであることはわかっており、基になるプロセスが正常であると仮定して、mixtoolsパッケージを使用してバイモーダル分布に適合 させました。

set.seed(99)  
dat.demand2050.mixmdl <- normalmixEM(dat.demand2050.unique, lambda=c(0.3,0.7), mu=c(60,70), k=2)

これにより、次の結果が得られます
ここに画像の説明を入力
（実線は曲線に適合し、破線は元の密度です）。 

# get the parameters of the mixture
dat.demand2050.mixmdl.prop <- dat.demand2050.mixmdl$lambda    #mix proportions
dat.demand2050.mixmdl.means <- dat.demand2050.mixmdl$mu    #modal means
dat.demand2050.mixmdl.dev <- dat.demand2050.mixmdl$sigma   #modal std dev

混合パラメーターは次のとおりです。 

>dat.demand2050.mixmdl.prop  #mix proportions  
[1] 0.2783939 0.7216061  
>dat.demand2050.mixmdl.means  #modal means  
[1] 56.21150 73.08389  
>dat.demand2050.mixmdl.dev  #modal std dev  
[1] 3.098292 6.413906

次の質問があります。

  1. 基礎となる分布を近似する新しい値のセットを生成するには、私のアプローチは正しいですか、それともより良いワークフローがありますか?
  2. 私のアプローチが正しい場合、この結果を使用して、この混合分布からランダムな値のセットを生成するにはどうすればよいですか?
4

2 に答える 2

6

あなたのサンプルサイズは、混合物に適合するかどうか少し疑わしいですが、気にしないでください。当てはめた混合物から次のようにサンプリングできます。

probs <- dat.demand2050.mixmdl$lambda
m <- dat.demand2050.mixmdl$mu
s <- at.demand2050.mixmdl$sigma

N <- 1e5
grp <- sample(length(probs), N, replace=TRUE, prob=probs)
x <- rnorm(N, m[grp], s[grp])
于 2013-07-29T14:29:40.817 に答える
4

あなたのアプローチは正しいです。

混合分布の各サンプルについて、2 つの成分ガウス分布のどちらからサンプルを取得するかを選択し、その分布からサンプルを抽出するだけです。

見つけた混合比率を使用して、2 つの分布から選択できます。0 と 1 の間の乱数をシミュレートし、乱数が最初の比率より小さい場合は最初の分布からサンプリングし、それ以外の場合は 2 番目の分布からサンプリングします。

最後に、rnorm 関数を使用して、関連するガウス分布からサンプリングします。

dat.demand2050.mixmdl.prop=c(0.2783939,0.7216061)
dat.demand2050.mixmdl.means=c(56.21150,73.08389)
dat.demand2050.mixmdl.dev=c(3.098292,6.413906)

sampleMixture=function(prop,means,dev){
    # Generate a uniformly distributed random number between 0 and 1
    # in order to choose between the two component distributions
    distTest=runif(1)
    if(distTest<prop[1]){
        # Then sample from the first component of the mixture
        sample=rnorm(1,mean=means[1],sd=dev[1])
    }else{
        # Sample from the second component of the mixture
        sample=rnorm(1,mean=means[2],sd=dev[2])
    }
    return(sample)
}

# Generate a single sample
sampleMixture(dat.demand2050.mixmdl.prop,dat.demand2050.mixmdl.means,dat.demand2050.mixmdl.dev)

# Generate 100 samples and plot resulting distribution
samples=replicate(100,sampleMixture(dat.demand2050.mixmdl.prop,dat.demand2050.mixmdl.means,dat.demand2050.mixmdl.dev))
plot(density(samples))
于 2013-07-29T13:34:12.137 に答える