-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path03a_estim_matern.R
144 lines (112 loc) · 5.99 KB
/
03a_estim_matern.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
# initalizations
source('00_system.R')
args <- commandArgs(trailingOnly = TRUE)
loadstring <- as.character(args[1])
print(loadstring)
# estimation
tt <- sum(.Internal(gc(FALSE, TRUE, TRUE))[13:14]) # memory report
elapsed_time <- system.time({
# get the data and simulationspecs
load(file = paste0(CommonPath, '/fields/', loadstring))
if(this_field@gaussian){
# Estimate betas (Oracle approach) -----;
distmat <- as.matrix(dist(this_field@locs))
Sigma_oracle <- cov.mat2(h = distmat, theta = this_field@truetheta)
Rstruct <- chol(Sigma_oracle)
z <- this_field@fulldata
n <- length(z)
X <- cbind(rep(1, length(this_field@fulldata)), cos(this_field@locs))
colnames(X) <- NULL
xsigma <- t(solve(Sigma_oracle, X))
right_side <- xsigma %*% matrix(z, ncol = 1)
beta0_init <- solve(xsigma %*% X, right_side)
z_nomean <- z - ( X %*% beta0_init )
neg2loglikelihood <- function(theta) {
Sigma <- cov.mat2(h = distmat, c(theta[1:2], this_field@truetheta[3], theta[3]))
cholS <- chol(Sigma, Rstruct = Rstruct)
logdet <- sum(log(diag(cholS)))
return(n * log(2 * pi) + 2 * logdet + sum(z_nomean * backsolve(cholS,
forwardsolve(cholS, z_nomean, transpose = TRUE, upper.tri = TRUE),
n)))
}
theta_lower <- c(1e-8, 1e-8, 1e-8)
theta_upper <- c(3 * sqrt(2), 20, 1)
theta_0 <- c(0.1, 0.9*var(z_nomean), 0.1*var(z_nomean))
cl <- makeCluster(ifelse(ncores > length(theta_0)*2 +1, length(theta_0)*2 + 1, 1 + floor(ncores/3)),
useXDR = F)
setDefaultCluster(cl = cl)
clusterExport(cl, c('lib.loc'))
clusterEvalQ(cl, .libPaths(new = lib.loc))
clusterEvalQ(cl, library('spam', lib.loc = lib.loc))
clusterEvalQ(cl, library('spam64', lib.loc = lib.loc))
clusterExport(cl, c('neg2loglikelihood', 'z_nomean', 'distmat', 'Rstruct', 'n', 'this_field', 'cov.mat2'))
output <- optimParallel::optimParallel(fn = neg2loglikelihood,
method = 'L-BFGS-B',
lower = theta_lower,
upper = theta_upper,
par = theta_0,
parallel = list(forward = FALSE,
loginfo = TRUE),
control = list(factr = 1e8, # approx. 1e-8
trace = TRUE, maxit = 1000))
stopCluster(cl)
} else {
z <- this_field@fulldata
X <- cbind(rep(1, length(this_field@fulldata)), cos(this_field@locs))
colnames(X) <- NULL
distmat <- as.matrix(dist(this_field@locs))
Sigma_oracle <- cov.mat2(h = distmat, theta = this_field@truetheta)
Rstruct <- chol(Sigma_oracle)
n <- length(z)
neg2logliksasbetas <- function(theta){
Sigma <- cov.mat2(h = distmat,
theta = c(theta[1:2], this_field@truetheta[3], theta[3]))
corr_matrixd <- Sigma/Sigma[1, 1]
cholS <- chol(corr_matrixd, Rstruct = Rstruct)
logdet <- 2 * sum(log(diag(cholS)))
stdata <- (z - X %*% matrix(theta[6:8], nrow = 3)) / sqrt(theta[2])
skew <- theta[4]
kurt <- theta[5]
Z <- sinh(kurt * asinh(stdata) - skew)
C <- sqrt(1+Z^2)
return( - 2 * (n * log(kurt) - (n/2) * log(2 * pi * (theta[2])) - 0.5 * logdet +
sum( log(C) - 0.5 * log(1 + stdata^2)) - 0.5 * sum(Z * backsolve(cholS,
forwardsolve(cholS, Z, transpose = TRUE, upper.tri = TRUE),
n)) )
)
}
theta_lower <- c(1e-8, 1e-8, 1e-8, -5, 1e-8, -20, - 20, -20)
theta_upper <- c(3 * sqrt(2), 20, 1, 5, 5, 20, 20, 20)
theta_0 <- c(0.1, 0.9 * var(this_field@fulldata), 0.1 * var(this_field@fulldata), 0, 1, mean(this_field@fulldata), 0, 0)
cl <- makeCluster(ifelse(ncores > length(theta_0)*2 +1, length(theta_0)*2 + 1, 1 + floor(ncores/3)),
useXDR = F)
setDefaultCluster(cl = cl)
clusterExport(cl, c('lib.loc'))
clusterEvalQ(cl, .libPaths(new = lib.loc))
clusterEvalQ(cl, library('spam', lib.loc = lib.loc))
clusterEvalQ(cl, library('spam64', lib.loc = lib.loc))
clusterExport(cl, c('neg2logliksasbetas', 'z', 'distmat', 'Rstruct', 'n', 'X', 'this_field', 'cov.mat2'))
output <- optimParallel::optimParallel(fn = neg2logliksasbetas,
method = 'L-BFGS-B',
lower = theta_lower,
upper = theta_upper,
par = theta_0,
parallel = list(forward = FALSE,
loginfo = TRUE),
control = list(factr = 1e8, # approx. 1e-8
trace = TRUE, maxit = 1000))
stopCluster(cl)
}
})[3]
ifelse(this_field@gaussian, betas_hat <- beta0_init,
betas_hat <- output$par[6:8])
memory_used <- sum(.Internal(gc(FALSE, FALSE, TRUE))[13:14]) - tt
this_estimates <- new('matern_estimate',
specsnr = this_field@specsnr,
fielditer = this_field@fielditer,
time = elapsed_time,
optimout = output,
betas_hat = betas_hat,
memory_used = memory_used,
loadstring = loadstring)
save(this_estimates, file =paste0(CommonPath, '/estimates/estimates.matern.', this_field@specsnr, '.iter', this_field@fielditer, '.RData'))