-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path049b_ID_Inf_Fig12_TST9d_TW_&_WL.R
206 lines (153 loc) · 5.74 KB
/
049b_ID_Inf_Fig12_TST9d_TW_&_WL.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
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
#============
# 9 Mar. 2024
#============
# Aim:
# use neg_logL_Matern function in 049 to do optimization
# Method:
source("049_ID_Inf_Fig12_TST9d_neg_logL_Matern.R")
#==========
# Settings
#==========
ini <- c(0.5, 0.5, 1, 2) # A, dlt, sig2, kappa
Vals <- c()
for (i in 1:length(all_pars_lst)){
value <- rep(ini[i], sum(is.na(all_pars_lst[[i]])))
Vals <- c(Vals, value)
}
## lower bound for each parameters,
# NA: no lower bound
# kappa > 0
lower_bound <- c(rep(NA, sum(is.na(all_pars_lst[[1]]))),
rep(0.05, sum(is.na(all_pars_lst[[2]]))),
rep(0.001, sum(is.na(all_pars_lst[[3]]))),
rep(0.001, sum(is.na(all_pars_lst[[4]]))),
rep(0.001, p))
# [1] NA NA NA 0.050
# [5] 0.050 0.050 0.001 0.001
# [9] 0.001 0.001 0.001 0.001
# [13] 0.001 0.001 0.001
#================
# Optim: Tir-Wave
#================
optm_pars_Matern_TW <- optim(par = c(Vals, rep(1, p)), # ini guess
fn = neg_logL_Matern,
p = p, data_str = hierarchy_data6,
all_pars_lst = all_pars_lst_Matern_6, df = df_TW,
fit_indx = fit_indx, b = "Tri-Wave",
method = "L-BFGS-B",
lower = lower_bound,
control = list(trace = 0,
maxit = 200,
pgtol = 1e-4))
# 1st run
# final value 1717.052764
#stopped after 3 iterations
optm_pars_Matern_TW$message
# [1] "NEW_X"
optm_pars_Matern_TW$counts
# function gradient
# 12 12
# 2nd run
#final value 1563.972442
#stopped after 85 iterations
optm_pars_Matern_TW$message
# [1] "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"
# initial value
optm_pars_Matern_TW$counts
# function gradient
# 256 256
optm_pars_Matern_TW$par
# [1] 0.48628902 0.54297915
#[3] 0.55838092 0.53285438
#[5] 0.49051338 0.50782285
#[7] 0.50001489 0.49953224
#[9] 0.49947753 0.08658082
#[11] 0.42727072 2.22067814
#[13] 1.05183670 0.42125199
#[15] 0.47794978 0.19150125
#[17] 0.56442973 0.65038259
#[19] 1.63223177 1.15320121
#[21] 2.07270914 1.55084153
#[23] 1.40549876 0.67715591
#[25] 2.57817159 1.36318412
#[27] 2.96140610 2.47055638
#[29] 1.55639986 1.17298566
#[31] 1.24305039 1.42695151
#[33] 1.63231830 1.44185847
#[35] 1.40587505 1.58684220
optm_pars_Matern_TW$convergence
# [1] 52
# 3rd run:
# use the returned pars as 3rd run initial value
# meanwhile change pgtol = 1e-4
optm_pars_Matern_TW$message
#[1] "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"
optm_pars_Matern_TW$convergence
#[1] 52
optm_pars_Matern_TW$counts
#function gradient
# 103 103
optm_pars_Matern_TW$value
# [1] 1648.297
optm_pars_Matern_TW$par
# [1] 0.4829771 0.4803281 0.5261407
#[4] 0.4673359 0.5177298 0.5000071
#[7] 0.4995936 0.5015581 0.4936938
#[10] 2.2549959 0.3976075 2.6997030
#[13] 0.5816490 0.6544571 0.4532069
#[16] 0.5492003 0.2490240 0.7504331
#[19] 1.3748299 1.1990764 1.2761730
#[22] 0.9253619 1.3559191 0.3015382
#[25] 1.3712188 1.2371381 1.2213733
#[28] 1.2601855 1.1369781 1.0350282
#[31] 1.3310886 1.2168976 1.6480171
#[34] 1.7890516 1.7061369 1.8597839
#============
# Conclusion
#============
# For non-cross MRF, Matern + DAG
# it's very difficulty to get the converged optimization results
# due to the kappa in Matern is highly unlikely to converge
## from this point of view, cross-MRF evaude the kappa
# and has a set parameters that are much easier to fit.
# During a line search, the objective function may be evaluated multiple times
# to find the step size that minimizes the objective function along the search direction
# it may evaluate the gradient multiple times within each iteration
# to determine the search direction
# The optimization algorithm may perform additional function evaluations
# or convergence criteria, such as changes in the objective function values
# or parameter value
# Q: what's the difference with iteration?
# a simplified overview of what typically happens in an iteration
# of an optimization algorithm
# 1. Current Solution
# 2. Compute Gradient: If the algorithm uses gradients (first-order derivatives)
# to guide the search, it computes the gradient of the objective function
# with respect to the parameters at the current solution point.
# This gradient provides information about the direction of steepest ascent
# (for maximization) or descent (for minimization) of the objective function.
# 3. Update Parameters: Based on the gradient information the algorithm updates
# the parameters (or decision variables) in the direction that is expected to improve the objective function value.
# 4. Convergence Check
# 5. Repeat or Terminate: If the convergence criteria are not met,
# the algorithm repeats the process by starting a new iteration
# with the updated parameters. Otherwise, if the convergence criteria
# are satisfied, the algorithm terminates, and the current solution
# is considered the optimal or near-optimal solution.
##================================
# Need to find a good inital start
##================================
# Grid search
#================
# Optim: Wendland
#================
optm_pars_Matern_TW <- optim(par = c(Vals, rep(1, p)), # ini guess
fn = neg_logL_Matern,
p = p, data_str = hierarchy_data6,
all_pars_lst = all_pars_lst_Matern_6, df = df_WL,
fit_indx = fit_indx,
method = "L-BFGS-B",
lower = lower_bound,
control = list(trace = 1,
maxit = 1,
pgtol = 1e-5))