2

我想根据以下公式从每月值模拟每日值：V(d)=V(m)*(1+noise)，噪声正态分布，均值为 0，标准差为 0.18)

我在名为 wind.m 的 12x1 矩阵中有每月值，在另一个名为 day 的 12x1 矩阵中有每个月的天数。我使用两个 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
4

2 回答 2

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 是您的模拟数据。

于 2013-06-28T21:00:20.287 回答
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 ...
于 2013-06-28T20:06:40.083 回答