2013-06-28 59 views
2

我想模擬每月值的公式: V(d)= V(m)*(1 +噪聲),噪聲通常以均值0和標準差0.18分佈)R中的兩個for循環

我有在另一個12X 1矩陣每個月12X 1矩陣稱爲wind.m和天月值稱爲day.I使用兩個for循環:

for (i in (1:12)) { 
    for (j in (1:12)){ 
    A<-wind.m[i,]*(1+rnorm(day[j,],0,0.18)) 
    } 
    print(A) 
} 

這個代碼的結果模擬12套31每日價值,這是錯誤的(2月有28天,4月,6月有30天等)我不知道如何解決我的代碼。

以下是所使用的數據:

> wind.m 
     [,1] 
[1,] 2.78 
[2,] 2.93 
[3,] 3.09 
[4,] 3.11 
[5,] 3.44 
[6,] 3.44 
[7,] 3.71 
[8,] 3.86 
[9,] 4.05 
[10,] 4.08 
[11,] 4.22 
[12,] 4.30 
> day 
     [,1] 
[1,] 31 
[2,] 28 
[3,] 31 
[4,] 30 
[5,] 31 
[6,] 30 
[7,] 31 
[8,] 31 
[9,] 30 
[10,] 31 
[11,] 30 
[12,] 31 
+0

我什至不能看到你是如何產生31塊..它看起來像我在做你12天的塊。 –

+0

這很沒有意義,尤其是沒有你的'day'和'wind.m'數據......你如何從任何事物得到31?它是在'天'嗎? 「 – Justin

+0

」以及另一個12x1矩陣每個月的日期稱爲日。「你的天矩陣不應該是12x31嗎? –

回答

2

你可以離開循環並解決您的問題,並瞭解了一下關於R矢量三個簡單的步驟。

設置:

wind.m <- c(2.78, 2.93, 3.09, 3.11, 3.44, 3.44, 
      3.71, 3.86, 4.05, 4.08, 4.22, 4.30) 
day <- c(31, 28, 31, 30, 31, 30, 
     31, 31, 30, 31, 30, 31) 

然後使用rep擴大你的月均值爲一年每日值。

day_means <- rep(wind.m, times = day) 

然後生成您的模擬值。

A <- rnorm(length(day_means), mean = day_means , sd=0.18) 

A是您的模擬數據。

+0

非常感謝! – user2472273

+0

+1不錯的使用R向量化 – dickoa

+0

不客氣,謝謝! – Kikapp

1

無需嵌套循環。我認爲你可以在這種情況下使用mapply

set.seed(42) 
df <- data.frame(month = month.abb, 
       numb_day = c(31, 28, 31, 30, 31, 30, 31, 
       31, 30, 31, 30 , 31), 
       wind_m = rnorm(12, mean = 20, sd = 5)) 

addNoise <- function(n, x) x * (1 + rnorm(n, 0, 0.18)) 

wind_d <- mapply(addNoise, df$numb_day, df$wind_m) 
str(wind_d) 
## List of 12 
## $ : num [1:31] 21.5 27.5 35.5 38.6 21.6 ... 
## $ : num [1:28] 21.7 11.3 15 16.2 12 ... 
## $ : num [1:31] 13.7 20.8 17.2 27.5 27.1 ... 
## $ : num [1:30] 26.3 14.2 20.2 23.1 17.1 ... 
## $ : num [1:31] 17.3 22.5 22 26.1 25.6 ... 
## $ : num [1:30] 20.8 12.8 13.1 15.5 18.3 ... 
## $ : num [1:31] 32.1 19.7 30.5 28 32.4 ... 
## $ : num [1:31] 19.5 19.5 20.1 21.6 19.1 ... 
## $ : num [1:30] 35.5 34.1 26.7 32.2 25.3 ... 
## $ : num [1:31] 22.7 19.5 15.8 21.7 24.3 ... 
## $ : num [1:30] 20.2 32.2 23.7 32.3 24.3 ... 
## $ : num [1:31] 41.2 31.1 28.4 35.7 28.7 ... 

編輯:使用OP的數據

set.seed(123) 
wind_d <- mapply(addNoise, day, wind.m) 
## List of 12 
## $ : num [1:31] 2.5 2.66 3.56 2.82 2.84 ... 
## $ : num [1:28] 2.77 3.4 3.39 3.36 3.29 ... 
## $ : num [1:31] 3.21 3.3 2.81 2.9 2.52 ... 
## $ : num [1:30] 3.67 3.42 3.24 2.76 3.87 ... 
## $ : num [1:31] 3.51 2.85 3.14 3.28 4.58 ... 
## $ : num [1:30] 3.92 3.65 2.82 3.37 3.27 ... 
## $ : num [1:31] 4.55 3.48 3.13 3.55 3.58 ... 
## $ : num [1:31] 4.72 3.5 3.17 5.02 3.55 ... 
## $ : num [1:30] 3.22 3.92 5.44 3.98 3.06 ... 
## $ : num [1:31] 3.06 4.7 3.75 4.21 4.13 ... 
## $ : num [1:30] 3.89 4.47 2.69 4.38 5.16 ... 
## $ : num [1:31] 4.59 2.71 3.24 4.14 4.97 ... 
+0

非常感謝! – user2472273

+0

無論如何,我可以將結果保存到365x1向量? – user2472273

+0

@ user2472273是的,你可以。只要做'unlist(wind_d)' – dickoa