老师要求的图


代码
memory<- 5
nFarm <- 20    # Number of farms in the system
nAct  <- 72     # Number of activities in the system
nYear <- 11    # Number of year in total (1999 - 2009)
region_names<- read.csv('region_names.csv', header = FALSE)
farm_names  <- as.character(region_names[1:nFarm,])
act_names   <- paste('activity', formatC(1:nAct, width=nchar(nAct), flag='0'), sep='')

####---------- I. Set up ranges of means for generating input parameters ----------####

input_values <- list(
  ref_income_mu  = list(min=100, max=400, randomFunc="qunif"),
  #ref_income_ratio_mu= list(min=0.4,  max=0.9),   # Mean of aspiration levels
  tolincome_mu =       list(min=0.01, max=0.3,randomFunc="qunif"),  # Mean of tolerance of dissimilarity in income change
  tolactivi_mu =       list(min=0.1,  max=0.75,randomFunc="qunif"),  # Mean of tolerance of dissimilarity in activity
  lambda_mu    =       list(min=1.50, max=4.00,randomFunc="qunif"),  # Mean of lambda in satisfaction calculation using CPT
  alpha_plus_mu=       list(min=0.50, max=1,randomFunc="qunif"),  # Mean of alpha_plus in satisfaction calculation using CPT
  alpha_minus_mu=      list(min=0.50, max=1,randomFunc="qunif"),  # Mean of alpha_plus in satisfaction calculation using CPT
  phi_plus_mu  =       list(min=0.50, max=1,randomFunc="qunif"),  # Mean of alpha_minus in satisfaction calculation using CPT
  phi_minus_mu =       list(min=0.50, max=1,randomFunc="qunif"),   # Mean of alpha_minus in satisfaction calculation using CPT
  price_mu           = list(min=300, max=400, randomFunc="qunif")         # Mean of price
)


# TODO: number of parameter sets (for LHS)
sample_count <- 200

# TODO: give number of bootstraps in SRC/SRRC
src_nboot <- 100

# TODO: SRC or SRRC (on ranks)
on_rank <- FALSE

# TODO: names of output values
output_names  <- c("opt-out","imitation","optimization","repetition") # Corresponding to c(1,2,3,4)
input_names   <- names(input_values)

# how many repetitions for each input factor set should be run (to control stochasticity)?
# TODO: adapt the number of repititions, set to 1 if deterministic model
#no.repeated.sim <- 10
num_repeated_simu <- 3  

# TODO: should R report the progress
trace_progress = FALSE
output_names  <- c("opt-out","imitation","optimization","repetition") # Corresponding to c(1,2,3,4)
input_names   <- names(input_values)


#### Load "sim_results_lhs_orig" and "lhs_design"

# transform the data
sim_results_lhs <- t(sim_results_lhs_orig)
output_names  <- c("opt-out","imitation","optimization","repetition") # Corresponding to c(1,2,3,4)
colnames(sim_results_lhs) <- output_names
sim_results_lhs <- cbind(as.data.frame(lhs_design), sim_results_lhs)

#-------------------------------------------------------------------------------------
# V. Run of SRC/SRRC
#-------------------------------------------------------------------------------------
require(sensitivity)

# iterate over different evaluation criteria in simulation results
# calculate SRC
src_list <- list()
for (o in output_names) {
  src_list[[o]] <- src(X=sim_results_lhs[,1:length(input_values)], y=sim_results_lhs[o], nboot = src_nboot, rank = FALSE)
}

# calculate SRRC
srrc_list <- list()
for (o in output_names) {
  srrc_list[[o]] <- src(X=sim_results_lhs[,1:length(input_values)], y=sim_results_lhs[o], nboot = src_nboot, rank = TRUE)
}


#-------------------------------------------------------------------------------------
# VI. Calculation of R? for the original data
#-------------------------------------------------------------------------------------
## Function calculating R squared
get.rsquare <- function(x, y, rank) {
  data <- data.frame(Y = y, x)
  if (rank) {
    for (i in 1:ncol(data)) {
      data[,i] <- rank(data[,i])
    }
  }
  i = 1:nrow(data)
  d <- data[i, ]
  lm.Y <- lm(formula(paste(colnames(d)[1], "~", paste(colnames(d)[-1], collapse = "+"))), data = d)
  return(summary(lm.Y)$r.squared)
}

## Calculate R squared for SRC
r_square_src <- list()
for (o in output_names) {
  r_square_src[[o]] <- get.rsquare(x=sim_results_lhs[,1:length(input_values)], y=sim_results_lhs[o], rank=FALSE)
}

print(r_square_src)

## Calculate R squared for SRRC
r_square_srrc <- list()
for (o in output_names) {
  r_square_srrc[[o]] <- get.rsquare(x=sim_results_lhs[,1:length(input_values)], y=sim_results_lhs[o], rank=TRUE)
}

print(r_square_srrc)

#-------------------------------------------------------------------------------------
# IV. Analysis of the results (postprocessing)
#-------------------------------------------------------------------------------------

# plot of package sensitivitiy (not shown in the paper)
# (package sesnitivity must be loaded)
for (o in output_names)
{   
  plot(src_list[[o]])
  title(sub=o)

  plot(srrc_list[[o]])
  title(sub=o)
}
library(sensitivity)

老师给的数据