I have the following non-linear DDE (delay differential equation) model solved with deSolve:
model <- function(parms, times = seq(0, 60, 1)) {
}
The known parameters are:
G <- structure(c(-1, 1, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0,
0, 1, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0,
0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, -1, 0,
0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0, -1, 0,
1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, -1,
0, 0, 1, 0, -1, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0,
0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0,
0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0,
0, -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, -1), dim = c(8L, 22L))
U <- structure(c(9.6, 0, 0, 1.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.01,
0, 0, 0, 0, 0, 0, 0, 0, 11.2, 0, 0, 0, 0.333333333333333, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0.01, 0, 0, 0, 0, 0, 0, 0, 0, 0.6, 0,
0, 0, 0, 11.8, 0, 0, 0, 0, 0, 0, 0, 0, 0.01, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0.045, 0, 0, 0, 0, 0, 0, 0, 0, 0.01, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.8, 0, 0, 0, 0, 0, 0,
0, 0.01, 0, 0, 0, 0, 0, 0, 0, 0.006, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0.01, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.01, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.01), dim = c(22L, 8L))
forcing_func <- approxfun(x = c(0:60),
y = c(3.789978749,5.316037452,5.185396143,6.558284554,5.549130548,5.712633241,
4.152581514,5.942285729,5.814900551,3.968558097,4.276266213,5.616072402,
6.407433514,1.79539372,0.633871866,5.303194797,4.099478638,4.137054321,
4.175807584,1.877168312,1.93214366,0.393710033,2.64911332,4.674492235,
3.698606684,3.442866381,3.73217674,4.638280031,3.754229767,2.865937731,
1.459804928,2.069766183,1.808161085,3.305388336,3.092063238,3.593206788,
3.8777349,2.819175366,3.55283987,4.51495667,2.009203206,2.304827672,
3.109626536,4.471240675,4.004597144,4.4204254,4.381900771,2.069175771,
-1.341448145,4.540838287,3.310323684,2.70840268,2.807892132,3.090647561,
3.077577537,2.805257996,3.38624182,3.286864016,2.591419861,-3.114741981,
5.874225293)/100)
The starting parameter values are (based on visual comparison between data and outputs):
parms <- list(
G = G,
U = U,
alpha_9 = 0.045,
alpha_11 = 1.8,
gamma = 1/3,
tau_1 = 8,
tau_2 = 2,
z_e = 0.75,
c = 2,
i = 3,
k_1 = 0.5,
k_2 = 0.6,
n = 0.1,
q = 0.02,
delta = 0.2)
This is my data:
data <- structure(list(label = c("u1", "u1", "u1", "u1", "u1", "u1",
"u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1",
"u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1",
"u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1",
"u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1",
"u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1",
"u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3",
"u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3",
"u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3",
"u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3",
"u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3", "u3",
"u3", "u3", "u3", "u3", "u3", "u4", "u4", "u4", "u4", "u4", "u4",
"u4", "u4", "u4", "u4", "u4", "s4", "s4", "s4", "s4", "s4", "s4",
"s4", "s4", "s4", "s4", "s4", "s8", "s8", "s8", "s8", "s8", "s8",
"s8", "s8", "s8", "s8", "s8"), time = c(0, 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, 0, 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,
0, 4, 9, 14, 19, 24, 29, 34, 39, 44, 49, 0, 4, 9, 14, 19, 24,
29, 34, 39, 44, 49, 0, 4, 9, 14, 19, 24, 29, 34, 39, 44, 49),
value = c(3.24406229750676, 3.24529921770188, 3.30999392307195,
4.90064495158511, 6.47113483169738, 7.39318318121391, 5.43023912454796,
5.99536558710986, 6.72344971662373, 6.09831566465937, 4.42937629378807,
3.97134606875211, 5.78036541991601, 5.4689247746114, 2.95374636834702,
3.31021629913642, 2.86113772997592, 2.5675164398153, 3.34789915159254,
3.34551666519448, 2.66765839921454, 2.33677898884761, 2.58119782921292,
2.28347874367683, 2.37466930175288, 2.0011187459536, 2.36939791414476,
3.24749629574441, 3.5769155243051, 3.21823760580411, 2.85565628815629,
2.73520083932854, 2.21677481653148, 2.75677120669056, 3.19434802321364,
2.54418699186992, 2.65048506014746, 2.01223641524736, 1.95144023986766,
2.27822759631491, 2.06042863359443, 2.06007595772787, 2.23510657453936,
3.37162931372549, 4.19484198023565, 7.47734649610678, 7.46145283018868,
6.76642007133593, 5.33651666666667, 7.53477958333333, 7.95332207207207,
7.22536031108826, 6.68377436416087, 6.34325092421442, 5.62865849506299,
5.17861436170213, 6.34115107913669, 6.40804546941446, 6.04034673366834,
6.30620105549881, 8.51649908592322, 3.979, 4.11, 4.179, 4.331,
4.531, 4.51, 4.755, 4.891, 5.568, 5.671, 5.958, 6.391, 6.58,
6.758, 6.841, 6.937, 6.937, 7.143, 7.206, 7.064, 7.327, 7.552,
7.669, 7.72, 7.819, 8.079, 8.397, 8.575, 9.163, 9.061, 8.921,
9.376, 9.381, 9.906, 9.906, 10.903, 11.34, 12.046, 12.562,
12.989, 13.532, 13.358, 13.419, 14.12, 14.682, 14.75, 15.218,
15.288, 15.631, 16.1, 16.295, 16.711, 18.272, 18.488, 19.133,
20.356, 20.238, 20.626, 20.552, 20.576, 2.3, 3, 3.5, 3.75,
4.75, 4.5, 5.1, 6.5, 7, 7.5, 8.75, 75, 90, 100, 125, 150,
175, 200, 225, 275, 300, 350, 40, 45, 50, 55, 60, 70, 75,
90, 100, 125, 140)), class = "data.frame", row.names = c("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", "510", "1010", "151", "201", "251",
"301", "351", "401", "451", "501", "123", "511", "1011", "152",
"202", "252", "302", "352", "402", "452", "502", "124", "512",
"1012", "153", "203", "253", "303", "353", "403", "453", "503"
))
Now, I have tried to fit some parameters based on the available data using the adaptive and delayed rejection MCMC algorithm from the FME package. However, I cannot seem to get it to converge:
Objective <- function(x, parset = names(x)) {
parms[parset] <- x
out <- model(parms) #,tout
## Model cost
return(modCost(obs = data, model = out, y = "value"))
}
Fit <- modMCMC(f = Objective,
p = c(alpha_9 = 0.045,
alpha_11 = 1.8,
z_e = 0.75,
c = 2,
i = 3,
k_1 = 0.5,
k_2 = 0.6,
n = 0.1,
q = 0.02,
delta = 0.2),
lower = c(0.01,
0.01,
0.01,
0.1,
0.1,
0.1,
0.1,
0.01,
0.001,
0),
upper = c(1,
10,
1,
10,
10,
1,
1,
3,
10,
10),
niter = 2000,
updatecov = 100,
ntrydr = 2)
How can the routine be modified to make it converge?