-
Notifications
You must be signed in to change notification settings - Fork 0
/
run_eBASCS_full_marginal_fixedk.R
168 lines (118 loc) · 5.13 KB
/
run_eBASCS_full_marginal_fixedk.R
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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
# Extended BASCS algorithm: MCMC for separating overlapping sources
# David Jones, Harvard University
# Luis Campos, Harvard University
####################################################################################
# FIXED K
# NO RJMCMC VERSION - ignore RJMCMC related notes
####################################################################################
####################################################################################
# MARGINAL ENERGY DISTRIBUTION
# Energy Distribution does not vary with time bins
####################################################################################
args=(commandArgs(TRUE))
eval(parse(text = args))
file.name = gsub('./sim_data/', '', file.name)
start <- proc.time()
############################
# Required R packages
############################
library("MASS")
library("MCMCpack")
library("cubature")
library("Rcpp")
#############################
# Load data
#############################
home = getwd()
results <- paste(home,"/results/",sep="")
data_dir = '/sim_data/'
initial_run = NULL
out.file.name = gsub('simulation', 'analysis_eBASCS', file.name)
load(paste(home, data_dir, file.name, sep = ''))
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# (1) Split data into component pieces and calculate image parameters
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
source('./additional_functions/_clean_data.R')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# (2) Set model options, number of breaks, psf, etc.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
source('./additional_functions/_model_options.R')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# (3) Prior parameter settings
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
source('./additional_functions/_prior_options.R')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# (4) MCMC seetings
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
source('./additional_functions/_mcmc_options.R')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# (5) Simulation Initialization, no burnin
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
source('./additional_functions/_init.time.sim.R')
#############################
# Code required
#############################
setwd(paste(home, "/additional_functions/", sep=""))
source(paste("mcmc.full",function_load_name,sep=""))
source(paste("mcmc.spatial.time.energy.full.marginal",function_load_name,sep=""))
source(paste("mcmc.extended.full",function_load_name,sep=""))
source(paste("mcmc.spatial.time.energy.extended.full.marginal",function_load_name,sep=""))
source(paste("extended_full_log_posterior",function_load_name,sep=""))
source(paste("full_log_posterior",function_load_name,sep=""))
source("logit.R")
source("inv.logit.R")
if (spectral_model=="full"){
fixedk_mcmc <- mcmc.full
fixedk_eBASCS_mcmc <- mcmc.spatial.time.energy.full
source("spectral_post_full.R")
} else {
fixedk_mcmc <- mcmc.extended.full
fixedk_eBASCS_mcmc <- mcmc.spatial.time.energy.extended.full.marginal
source("spectral_post_extended_full.R")
}
setwd(home)
#############################
# Main MCMC run
#############################
online_ordering <- "reference"
no_guess <- k_curr
mu_guess <- mu_curr[,order(w[-mix_num],decreasing=TRUE)]
adapt_end <- 0 # Iteration number on which to stop adaptive MCMC location updates
paras <- vector('list', main_mcmc_iters)
store_log_post <- vector('list', main_mcmc_iters)
store_log_like <- vector('list', main_mcmc_iters)
store_alloc <- allocate_curr*0
for (tt in 1:main_mcmc_iters){
# Print iteration at intervals
if (tt/print_interval == round(tt/print_interval)){
print(paste("Iteration number: ",tt,sep=""))
}
main_run_time <- fixedk_eBASCS_mcmc(mcmc_runs,online_ordering,tt,w,allocate_curr,mu_curr,eparas_all,ewt_all,k_curr,mix_num,bk,num_time_breaks, lambda, time_bin)
# extract allocation and add it to the overall allocation to get allocation rates
store_alloc = store_alloc + main_run_time[[2]]
# Extract parameters and allocation matrix after burnin
new_paras <- main_run_time[[1]]
k_curr <- new_paras[1]
mix_num <- k_curr+1
mu_curr <- matrix(new_paras[2:(2*k_curr+1)],2,k_curr)
w <- new_paras[(2*(k_curr+1)):(3*k_curr+2)]
eparas_all <- main_run_time[[3]]
ewt_all <- main_run_time[[4]]
eparas_all <- main_run_time[[3]]
ewt_all <- main_run_time[[4]]
lambda <- main_run_time[[5]]
time_bin <- main_run_time[[6]]
bk <- main_run_time[[7]]
num_time_breaks <- main_run_time[[8]]
allocate_curr <- main_run_time[[2]]
predicted = factor((allocate_curr %*% c(1, 2, 3))[,1], levels = 1:3)
confusion = table(predicted, truth = dat[,5])
paras[[tt]] <- list(par = new_paras, eparas_all = eparas_all, ewt_all = ewt_all,lambda = lambda, confusion = confusion)
}
main_run_time <- list(paras, allocate_curr,
alloc_rate = store_alloc/main_mcmc_iters)
# Save output
save(list=ls(),file=paste(results,out.file.name,sep=""))
# Show time taken
end <- proc.time()
print(end-start)