mcSimulation: Perform a Monte Carlo simulation.

Description

This function generates a random sample of an output distribution defined as the transformation of an input distribution by a mathematical model, i.e. a mathematical function. This is called a Monte Carlo simulation. For details cf. below.

Usage

mcSimulation(
  estimate,
  model_function,
  ...,
  numberOfModelRuns,
  randomMethod = "calculate",
  functionSyntax = "data.frameNames",
  relativeTolerance = 0.05,
  verbosity = 0
)

Value

An object of class mcSimulation, which is a list with elements:

$x

data.frame containing the sampled $x -$ (input) values which are generated from estimate.

$y

data.frame containing the simulated $y -$ (output) values, i.e. the model function values for x.The return of the decision model function may include the assignment of names for the output variables, e.g. like this:


           profit <- function(x){
            revenue - costs
            return(list(Revenue = revenue,
                   Costs = cost))
         }

Arguments

estimate: estimate: estimate of the joint probability distribution of the input variables. This can be read from a csv file and calculated with the estimate_read_csv function.
model_function: function: The function that transforms the input distribution. It has to return a single numeric value or a list with named numeric values.
...: Optional arguments of model_function.
numberOfModelRuns: The number of times running the model function.
randomMethod: character: The method to be used to sample the distribution representing the input estimate. For details see option method in random.estimate.
functionSyntax: character: The syntax which has to be used to implement the model function. Possible values are "data.frameNames", "matrixNames" or "plainNames". Details are given below.
relativeTolerance: numeric: the relative tolerance level of deviation of the generated confidence interval from the specified interval. If this deviation is greater than relativeTolerance a warning is given.
verbosity: integer: if 0 the function is silent; the larger the value the more verbose is output information.

Details

This method solves the following problem. Given a multivariate random variable $x = (x_1,\ldots,x_k)$ with joint probability distribution $P$, i.e. $$x \sim P.$$ Then the continuous function $$f:R^k \rightarrow R^l, y = f(x)$$ defines another random variable with distribution $$y \sim f(P).$$ Given a probability density $\rho$ of x that defines $P$ the problem is the determination of the probability density $\phi$ that defines $f(P)$. This method samples the probability density $\phi$ of $y$ as follows: The input distribution $P$ is provided as estimate. From estimate a sample x with numberOfModelRuns is generated using random.estimate. Then the function values $y=f(x)$ are calculated, where $f$ is model_function.

functionSyntax defines the syntax of model_function, which has to be used, as follows:

"data.frameNames"

The model function is constructed, e.g. like this:


         profit<-function(x){
           x[["revenue"]]-x[["costs"]]
         }

or like this:


         profit<-function(x){
           x$revenue-x$costs
         }

"matrixNames"

The model function is constructed, e.g. like this:


           profit<-function(x){
             x[,"revenue"]-x[,"costs"]
           }

"plainNames"

model_function is constructed, e.g. like this:


         profit<-function(x){
           revenue-costs
         }

Note: this is the slowest of the possibilities for functionSyntax.

Examples

Run this code

 #############################################################
 # Example 1 (Creating the estimate from the command line):
 #############################################################
 # Create the estimate object:
 variable=c("revenue","costs")
 distribution=c("norm","norm")
 lower=c(10000,  5000)
 upper=c(100000, 50000)
 costBenefitEstimate<-as.estimate(variable, distribution, lower, upper)
 # (a) Define the model function without name for the return value:
 profit1<-function(x){
   x$revenue-x$costs
 }
 # Perform the Monte Carlo simulation:
 predictionProfit1<-mcSimulation( estimate=costBenefitEstimate, 
                                 model_function=profit1, 
                                 numberOfModelRuns=10000,
                                 functionSyntax="data.frameNames")
 # Show the simulation results:
 print(summary(predictionProfit1))
 hist(predictionProfit1,xlab="Profit")
 #############################################################
 # (b) Define the model function with a name for the return value:
 profit1<-function(x){
   list(Profit=x$revenue-x$costs)
 } 
 # Perform the Monte Carlo simulation:
 predictionProfit1<-mcSimulation( estimate=costBenefitEstimate, 
                                 model_function=profit1, 
                                 numberOfModelRuns=10000,
                                 functionSyntax="data.frameNames")
 # Show the simulation results:
 print(summary(predictionProfit1, classicView=TRUE))
 hist(predictionProfit1) 
 #########################################################
 # (c) Using plain names in the model function syntax
 profit1<-function(){
   list(Profit=revenue-costs)
 } 
 # Perform the Monte Carlo simulation:
 predictionProfit1<-mcSimulation( estimate=costBenefitEstimate, 
                                 model_function=profit1, 
                                 numberOfModelRuns=1000,
                                 functionSyntax="plainNames")
 # Show the simulation results:
 print(summary(predictionProfit1, probs=c(0.05,0.50,0.95)))
 hist(predictionProfit1) 
 #########################################################
 # (d) Using plain names in the model function syntax and
 #     define the model function without name for the return value:
 profit1<-function() revenue-costs
 # Perform the Monte Carlo simulation:
 predictionProfit1<-mcSimulation( estimate=costBenefitEstimate,
                                  model_function=profit1,
                                  numberOfModelRuns=1000,
                                  functionSyntax="plainNames")
 # Show the simulation results:
 print(summary(predictionProfit1, probs=c(0.05,0.50,0.95)))
 hist(predictionProfit1, xlab="Profit")
 #############################################################
 # Example 2(Reading the estimate from file):
 #############################################################
 # Define the model function:
 profit2<-function(x){
   Profit<-x[["sales"]]*(x[["productprice"]] - x[["costprice"]])
   list(Profit=Profit)
 }  
 # Read the estimate of sales, productprice and costprice from file:
 inputFileName=system.file("extdata","profit-4.csv",package="decisionSupport")
 parameterEstimate<-estimate_read_csv(fileName=inputFileName)
 print(parameterEstimate)
 # Perform the Monte Carlo simulation:
 predictionProfit2<-mcSimulation( estimate=parameterEstimate, 
                                 model_function=profit2, 
                                 numberOfModelRuns=10000,
                                 functionSyntax="data.frameNames")
 # Show the simulation results:
 print(summary(predictionProfit2))
 hist(predictionProfit2)