如何在R语言中构造一个易于扩展的montecarlo模型

Question

我有一个简单的模型，公司有两个农场，每个农场种植两种作物（苹果和梨）。第一步就是把树的数量乘以每棵树上的果实数量。在

每棵树上的果实数量是模拟的（跨农场和作物）。在

在R中建模时，我至少要面对三个决策：

• 如何构造变量
• 如何模拟
• 如何将模拟变量与非模拟变量相乘

即使我添加了另一种作物和/或农场，我也希望它能起作用——理想情况下，即使我添加了另一个维度，例如作物品种（Granny Smith等）。我还想用名字来指农场和农作物，而不是用索引号。在

这是我想出的方法。它是可行的，但是很难添加另一个维度，而且需要大量的代码。有更整洁的方法吗？在

要构建变量：

farms <- c('Farm 1', 'Farm 2');
crops <- c('Pear', 'Apple');
params <- c('mean','sd');

numTrees <- array(0, dim=c(length(farms), length(crops)), dimnames=list(farms,crops));
fruitPerTree <- array(0, dim=c(length(farms), length(varieties), length(params)), 
                      dimnames=list(farms,varieties,params));

# input data e.g.
numTrees['Farm 1', 'Pear'] = 100
# and
fruitPerTree['Farm 1', 'Pear', 'mean'] = 50

要模拟：

^{pr2}$

要将simFruitPerTree和numTrees相乘，我发现首先必须执行手动广播：

simNumTrees <- array(numTrees, dim=c(length(dims[[1]]), length(dims[[2]]), numSims), 
                     dimnames=list(dims[[1]], dims[[2]], c(1:numSims)));
simTotalFruit <- simFruitPerTree * simNumTrees;

为了进行比较，在Python中（我比R更了解它），我可以在几行代码中执行所有这些步骤，方法是使用元组索引字典，以及两种字典理解，例如：

fruit_per_tree = {}
fruit_per_tree[('Farm 1', 'Pear')]  = (50, 15) # normal params
sim_fruit_per_tree = {key: random.normal(*params, size=num_sims) 
                      for key, params in fruit_per_tree.items() }
sim_total_fruit = {key: sim_fruit_per_tree[key]*num_trees[key] for key in num_trees }

在R里也有简单的方法吗？谢谢！在

Answer 1

如果我没弄错的话，你是在建模n个农场的水果总产量，每个农场有k种作物（这里，n=k=2）。每个农场都有一定数量的每种树木，每个农场的生产力（果实/树）是分布在N（μ，σ）上的随机变量，其中μ和σ取决于农场和品种。在

所以对于输入，我们构造一个数据帧，params，有5列：farm, crop, trees, mean, and sd。然后每行包含树的数量，每棵树的平均生产力，和每棵树在给定的农场/作物组合中的生产力变化。这些是输入。在

如果我们在树的水平上建模，那么给定农场中给定品种的每棵树的果实产量为：

rnorm(trees,mean,sd)

也就是说，输出是长度=树的随机样本，平均值和标准差适用于给定的品种和农场。那么这个品种/农场的所有树木的总产量就是上面向量的和，总产量就是所有农场/作物的总和。在

所有这些都给了我们蒙特卡罗模型的1迭代。为了确定总产出的分布，我们必须重复这个过程若干次。幸运的是，在R中这是相当简单的：

^{pr2}$

这个代码对于农场或品种的数量是不可知的；只需在开头更改farms和{}向量。如果不是所有农场都有所有的品种，请为缺失的品种设置params$trees <- 0。在

我们可以检查n.iterations的影响，如下所示。这段代码只运行完整的模拟100、1000和10000次，并使用ggplot绘制分布图。在

gg     <- do.call(rbind,
                  lapply(c(100,1000,10000),
                         function(n)cbind(n=n,total=sapply(1:n,output,params))))
gg     <- data.frame(gg)
library(ggplot2)
ggplot(gg)+
  geom_histogram(aes(x=total, y=..density.., fill=factor(n)))+
  scale_fill_discrete("Iterations")+
  facet_wrap(~n)

最后，我建议您考虑每棵树的输出比正常情况下更可能是泊松分布。如果使用rpois(...)而不是rnorm(...)重新运行模拟，则总体sd略低（约150而不是~200）。在

Answer 2

这是我问题的一般解决办法。我从罗兰的方法开始，并对其进行了更新，使分布、参数和尺寸都可以轻松更改。在

distSim <- function(nSims, simName, distFn, dimList, paramList, constList) {
    #
    # Simulate from a distribution across all the dimensions.
    #
    # Args:
    #   nSims:     integer, e.g. 10000
    #   simName:   name of the output column, e.g. 'harvestPerTree'
    #   distFn:    distribution function, e.g. rnorm
    #   dimList:   list of dimensions, 
    #              e.g. list(farms=c('farm A', 'farm B'), crops=c('apple', 'pear', 'durian'))
    #   paramList: list of parameters, each of length = product(length(d) for d in dimList),
    #              to be passed to the distribution function,
    #              e.g. list(mean=c(10,20,30,5,10,15), sd=c(2,4,6,1,2,3))
    #   constList: optional vector of length = product(length(d) for d in dimList)
    #              these are included in the output dataframe
    #              e.g. list(nTrees=c(10,20,30,1,2,3))
    #
    # Returns:
    #   a dataframe with one row per iteration x product(d in dimList)
    #

    # expand out the dimensions into a dataframe grid - one row per combination
    df <- do.call(expand.grid, dimList);
    nRows <- nrow(df);
    # add the parameters, and constants, if present
    df <- cbind(df, paramList);
    if (!missing(constList)) {
        df <- cbind(df, constList);
    }
    # repeat this dataframe for each iteration of the simulation
    df <- do.call("rbind",replicate(nSims, df, simplify=FALSE));
    # add a new column giving the iteration number ('simId')
    simId <- sort(rep(seq(1:nSims),nRows));
    df <- cbind(simId, df);
    # simulate from the distribution
    df[simName] <- do.call(distFn, c(list(n=nrow(df)), df[names(paramList)]))
    return(df);
}

示例用法（仅为简单起见，使用正态分布）：

^{pr2}$

样本输出：

   simId  farms  crops mean sd numTrees harvestPerTree
1      1 farm A  apple   10  2       10       9.602849
2      1 farm B  apple   20  4       20      20.153225
3      1 farm A   pear   30  6       30      25.839825
4      1 farm B   pear    5  1        1       1.733120
5      1 farm A durian   10  2        2      13.506135
6      1 farm B durian   15  3        3      11.690910
7      2 farm A  apple   10  2       10       7.678611
8      2 farm B  apple   20  4       20      22.119560
9      2 farm A   pear   30  6       30      31.488002
10     2 farm B   pear    5  1        1       5.366725
11     2 farm A durian   10  2        2       6.333747
12     2 farm B durian   15  3        3      13.918085
13     3 farm A  apple   10  2       10      10.387194
14     3 farm B  apple   20  4       20      21.086388
15     3 farm A   pear   30  6       30      34.076926
16     3 farm B   pear    5  1        1       6.159415
17     3 farm A durian   10  2        2       8.322902
18     3 farm B durian   15  3        3       9.458085

还请注意，您还可以以一种良好的索引方式定义输入值；例如，如果您定义

numTrees2 <- array(0, dim=c(length(farms), length(crops)), dimnames=tree_dimList);
numTrees2['Farm A','apple'] = 200; 
# etc

然后c(numTrees)对其输出排序的方式与expand.grid相匹配，因此可以将constList = list(numTrees=c(numTrees2)传递给distSim。在

Answer 3

下面是我如何设置这样的模拟：

#for reproducibility
set.seed(42)

#data
farms <- data.frame(farm=rep(1:2, each=2),
                    trees=sample(100, 4),
                    crop=rep(c("pear", "apple")),
                    mean=c(100, 200, 70, 120),
                    sd=c(30, 15, 25, 20))

#n
n <- 100

#simulation
fruits <- t(matrix(rnorm(n*nrow(farms), farms$mean, farms$sd), ncol=n))

#check 
colMeans(fruits)
#[1] 101.10215 200.06649  68.01185 120.05096

library(reshape2)
fruits <- melt(fruits, value.name="harvest_per_tree")
farms$i <- seq_len(nrow(farms))

farm_sim <- merge(farms, fruits, by.x="i", by.y="Var2", all=TRUE)
names(farm_sim)[7] <- "sim_i"

#multiply with number of trees
farm_sim$harvest_total <- farm_sim$harvest_per_tree * farm_sim$trees
head(farm_sim)
#   i farm trees crop mean sd sim_i harvest_per_tree harvest_total
# 1 1    1    92 pear  100 30     1        110.89385     10202.234
# 2 1    1    92 pear  100 30     2        145.34566     13371.801
# 3 1    1    92 pear  100 30     3        139.14609     12801.440
# 4 1    1    92 pear  100 30     4         96.00036      8832.033
# 5 1    1    92 pear  100 30     5         26.78599      2464.311
# 6 1    1    92 pear  100 30     6         94.84248      8725.508

library(ggplot2)
ggplot(farm_sim, aes(x=sim_i, y=harvest_total, colour=factor(farm))) +
  geom_line() +
  facet_wrap(~crop)