R/model_anything_byMonteCarloSimulation.R

In msxakk89/dat: Data Analyst Toolkit (dat): a library of tools for data analysts and scientists

#' Simulate anything by Monte Carlo assuming normal or skew-normal distributions and linear relationship
#' 
#' By relying on Law of Large Numbers estimate the statistical distribution by Monte Carlo simulation. 
#' 
#' The assumption is that the object to be modeled can be so by linear function with single or multiple terms (tasks) where the terms are summed-up. In the case of computing the time required to complete a project this means each term (task) is dependent on the completion of other sub-task.
#' An underlying statistical distribution for the individual terms has to be assumed. In this case it is either normal or skew-normal distribution.
#' For example the model can be TIME_TOTAL = TIME_TASK_1 + TIME_TASK_2 + TIME_TASK_3 which would have its P(TIME_TOTAL) = P(TIME_TASK_1) + P(TIME_TASK_2) + P(TIME_TASK_3) described by parameters of normal os skewed-normal distribution.
#' An inverse of cumulative density function (iCDF) is used to compute the simulated values corresponding to randomly generated cumulative probability (within the range of 0 and 1).
#' 
#' @param mod A list whose elements are vectors describing the normal or skew normal probability distributions which will be used in Monte Carlo simulation. For normal distribution the length of vector must be 2 with mean and st.dev being 1st and 2nd element respectively. For skew-normal distribution the length of vector must be 3 with mean (mu), st.dev (omega), and skewness (alpha)  being 1st, 2nd and 3rd element respectively.  
#' @param spins A numeric of count of spins to be performed for the Monte Carlo simulation.
#' @param verbose A boolean indicating whether to execute in verbose mode.
#' @return A numeric vector of simulated values.
#' @export
#' 

model_anything_byMonteCarloSimulation = function(mod,spins,verbose=T){
  out = c()
  for(i in 1:spins){
    if(verbose){
      print(paste("Spin count: ",toString(i)))
    }
    # Initiating a new instance of simulated sum of terms to be appended to output
    simulated_sum_of_terms = c()
    # Performing simulation for each term in the overall linear model
    # Initating a sample file 
    sample_file = c()
    for(i in 1:length(mod)){
      # Access the probability distribution of the term 
      term_prob_dist = mod[[i]]
      # Generating a random number within interval 0-to-1 (sampled from uniform distribution) for the considered term
      random = runif(1)
      # Check if the term is characterised by normal distribution
      if(length(term_prob_dist)==2){
        if(verbose){
          print("Simulating for normal distribution term")
        }
        mean = term_prob_dist[1]
        sd = term_prob_dist[2]
        # Computing the output for inverse cumulative distribution function
        simulated = qnorm(p = random, mean = mean, sd = sd)
        sample_file = c(sample_file, simulated)
      }
      # Else it must be a skew-normal distribution
      else{
        if(verbose){
          print("Simulating for skewed-normal distribution term")
        }
        mean = term_prob_dist[1]
        sd = term_prob_dist[2]
        skewness = term_prob_dist[3]
        simulated = qsn(p = random, xi = mean, omega = sd, alpha = skewness)
        sample_file = c(sample_file, simulated)
      }
    }
    simulated_sum_of_terms = sum(sample_file)
    out = c(out, simulated_sum_of_terms)
  }
  return(out)
}