2016-04-11 6 views
0

に私が1.48を超えると0.67を下回っている拒絶領域にあるどのように私の F値の数を確認するために1000回のシミュレーションをしようとしています。ループ内でインクリメントしますか? R

私はこれを持っているが、彼らが必要として変数は増加しません。最終目標は

は、私も試してみましたBの合計は見つけることです

for (k in 1:1000){ 
    Adata = rnorm(100, mean = 30, sd = 10) 
    Bdata = rnorm(100, mean = 45, sd = 10) 
    f = (sd(Bdata)^2)/(sd(Adata)^2) 
    if (f > 1.48){ 
    a = 0 
    a <- a + 1} 
    if (f < .67){ 
    b = 0 
    b <- b + 1} 
} 

a 
[1] 1 
b 
[1] 1 

for (k in 1000){ 
    Adata = rnorm(100, mean = 30, sd = 10) 
    Bdata = rnorm(100, mean = 45, sd = 10) 
    f = (sd(Bdata)^2)/(sd(Adata)^2) 
    a = f > 1.48 
    b = f < .67 
    } 
y = sum(a)+sum(b) 
y 
[1] 0 

f 'のは拒絶領域にありますか?

+0

代わりにグローバル変数を作成してみてください。汚いですが、問題を解決します。 –

答えて

1

最初の例では、if文がtrueになるたびにabをゼロにリセットします。したがって、最大値は常に1になります。

は、修正し、それらの行を並べ替えるには:

a = 0 #initialize outside of the loop 
b = 0 #initialize outside of the loop 
set.seed(1) # added for SO as you are using rnorm, remove this when you run your simulations 
for (k in 1:1000){ 
    Adata = rnorm(100, mean = 30, sd = 10) 
    Bdata = rnorm(100, mean = 45, sd = 10) 
    f = (sd(Bdata)^2)/(sd(Adata)^2) 
    if (f > 1.48){ 
    a <- a + 1} 
    if (f < .67){ 
    b <- b + 1} 
} 

私は今取得= 13、B = はRでこのような変数をインクリメントしていない、と述べた

29。行列やベクトル演算を利用することができます。

まず作成したシミュレーション行列

set.seed(1) 
Adata = matrix(data = rnorm(100*1000, mean = 30, sd = 10), nrow = 1000, ncol = 100) 
Bdata = matrix(data = rnorm(100*1000, mean = 30, sd = 10), nrow = 1000, ncol = 100) 

は次に、各ラインのためのあなたのfスコアを計算:

f <- apply(Bdata,1,function(x){sd(x)^2})/apply(Adata,1,function(x){sd(x)^2})

は今、あなたは、単に使用することができます。

sum(f > 1.48) 
[1] 15 

とを:

sum(f < .67) 
[1] 25 
1

最初のコードブロックでは、aとbを0回ごとに0にリセットし、おそらく1を追加します(次の反復で0に設定されるため、1になります)。

2番目のブロックでは、aとbをTRUEまたはFALSEに設定しますが、値を上書きしているため、最終反復の値のみが表示されます(実際には、 1:1000の場合は最後の反復しか表示されません)。

単純な解決策は、ループの外側でa=0b=0(またはそれ以上の場合はa <- 0b <- 0)を移動することです。

より良いアプローチは、機能のapply家族で何かを使用することです。

私のような何かを示唆している:

out <- replicate(1000, { 
    Adata = rnorm(100, mean = 30, sd = 10) 
    Bdata = rnorm(100, mean = 45, sd = 10) 
    (sd(Bdata)^2)/(sd(Adata)^2) 
    }) 

sum(out > 1.48) 
sum(out < 0.67) 

sum(out > 1.48 | out < 0.67) 
関連する問題