Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#######################################################################################
# Author: b.arendarczyk@sms.ed.ac.uk
# Purpose: Sample posterior distribution of parameters in PLUMv2 using
# Approximate Bayesian Sampling with Markov Chain Monte Carlo (ABC-MCMC)
#######################################################################################
# LAND_CHANGE_COST = function() {runif(1, 0, 1)}
# CROP_INCREASE_FACTOR = function() {runif(1, 0, 2)}
# CROP_DECREASE_FACTOR = function() {runif(1, 0, 2)}
EPSILON <- 0.5 # proposal acceptance threshold
SIGMA <- 0.1 # SD of proposal kernel
land_covers <- c("cropland", "pasture", "timberForest", "carbonForest", "natural")
start_values <- data.frame(fromLC = character(), toLC = character(), Cost = numeric())
# Generate land cover change pairs
n_lc <- length(land_covers)
for (i in 1:n_lc) {
for (j in 1:n_lc) {
start_values[(i - 1) * n_lc + j, 1:2] <- c(land_covers[i], land_covers[j])
start_values[(i - 1) * n_lc + j, 3] <- runif(1, 0, 1)
}
}
# Write land cover change file
write_param_file <- function(proposal, out_folder) {
prop_vals <- start_values
prop_vals$Cost <- proposal
#colnames(prop_vals) <- c("fromLC", "toLC", "Cost")
out_dir <- paste0("/exports/csce/eddie/geos/groups/LURG/models/PLUM/output/", out_folder ,"/conversion_costs.csv")
write.csv(prop_vals, out_dir, row.names = F)
}
# Calculate output statistics
calc_stats <- function(out_dir) {
}
# Calculate distance of stats between generated and observed
calc_dist <- function(stats) {
}
run_ABC_MCMC <- function(start_values, iterations) {
chain <- array(dim = c(iterations + 1, nrow(start_values)))
chain[1, ] <- start_values[, 3]
for (i in 1:iterations) {
# Generate new proposal for params
proposal <- rnorm(nrow(start_values), mean = chain[i,], sd = SIGMA)
# Make new output dir
out_folder <- paste0("abc_mcmc/iter", i)
out_dir <- paste0("/exports/csce/eddie/geos/groups/LURG/models/PLUM/output/", out_folder)
mkdir_comm <- paste0("mkdir ", out_dir)
system(mkdir_comm)
return()
# Write proposed params to file
write_param_file(proposal, out_folder)
# Run PULMv2
run_plum_str <- paste0("qsub /exports/csce/eddie/geos/groups/LURG/models/PLUM/plumv2/scripts/runPlum.sh -s ", out_dir, " -r y -p n -b Bart")
system(run_plum_str)
# Generate summary statistics {s1, s2, ...}
stat <- calc_stats(out_dir)
# Calculate distance between observed and simulated statistics d
d <- calc_dist(stat)
# Accept or reject proposal
if (d < EPSILON) {
chain[i+1, ] <- proposal
} else {
chain[i+1, ] <- chain[i, ]
}
}
return(chain)
}
run_ABC_MCMC(start_values, 1)