############################################################# ## CEoptim: Cross-Entropy R package for Optimization ############################################################# library("CEoptim") ## 5.1 Unconstrained Minimization of the Griewank function griewank <- function(X) { return(1 + sum(X^2)/4000 - prod(cos(X/sqrt(1:length(X))))) } set.seed(1234) mu0 <- c(5, 5, 5, 5, 5) sigma0 <- c(20, 20, 20, 20, 20) res <- CEoptim(griewank, continuous = list(mean = mu0, sd = sigma0), rho = 0.1, N = 1000L, verbose = TRUE, noImproveThr = Inf) res res$optimum res$optimizer$continuous print(res, optimum = TRUE) ## 5.2 Constrained Minimization of the Griewank function set.seed(123) griewank <- function(X) { return(1 + sum(X^2)/4000 - prod(cos(X/sqrt(1:length(X))))) } A <- rbind(c(0, 1), c(-1, -1), c(1, -1)) b <- c(4, -4, 4) res <- CEoptim(griewank, continuous = list(mean = c(0, 0), sd = c(10, 10), conMat = A, conVec = b), rho = 0.1, N = 200L, verbose = TRUE, noImproveThr = Inf) cat("direct optimizer =", res$optimizer$continuous, "\n") cat("direct minimum =", res$optimum, "\n") griewank.penalty <- function(X, A, b) { fn <- griewank(X) if (any(A %*% as.vector(X) > b)) { penalty <- norm(A %*% as.vector(X) - b) fn <- fn + 100 * penalty } return(fn) } set.seed(123) res.pen <- CEoptim(griewank.penalty, f.arg = list(A, b), continuous = list(mean = c(0, 0), sd = c(10, 10)), rho = 0.01, N = 2000L, verbose = TRUE, noImproveThr = Inf) cat("penalty minimizer =", res.pen$optimizer$continuous, "\n") cat("penalty minimum =", griewank(res.pen$optimizer$continuous), "\n") griewank.contour <- function(x, y) { return(1 + (x^2 + y^2)/4000 - cos(x/sqrt(1)) * cos(y/sqrt(2))) } x <- seq(-2, 10, length.out = 100) y <- seq(-1, 5, length.out = 100) z <- outer(x, y, FUN = griewank.contour) par(mai = c(0.5, 0.5, 0.3, 0.3), cex.axis = 0.8) contour(x, y, z) K <- dim(A)[1] for (i in 1:K) { if (A[i, 2] == 0) den <- 1e-100 else den <- A[i, 2] abline(b[i]/den, -A[i, 1]/den, lwd = 2) } points(res$optimizer$continuous[1], res$optimizer$continuous[2], pch = 4, lwd = 3) polygon(x = c(0, 4, 8), y = c(4, 0, 4), density = c(18, 0), lwd = 0.8) ## 5.3 Non-linear regression library("deSolve") FN <- function(t, state, parameters) { with(as.list(c(state, parameters)), { dV <- c * (V - V^3/3 + R) dR <- -1/c * (V - a + b * R) list(c(dV, dR)) }) } ssres <- function(x, fundf, times, y) { parameters <- c(a = x[1], b = x[2], c = x[3]) state <- c(V = x[4], R = x[5]) out <- ode(y = state, times = times, func = fundf, parms = parameters) return(sum((out[, 2] - y)^2)) } set.seed(123405) times <- seq(0, 20, by = 0.05) data("FitzHugh", package = "CEoptim") X <- c(0.2, 0.2, 3, -1, 1) parameters <- c(a = X[1], b = X[2], c = X[3]) state <- c(V = X[4], R = X[5]) out <- ode(y = state, times = times, func = FN, parms = parameters) plot(times, FitzHugh, cex = 1, xlab = "", ylab = "") lines(times, out[, 2], col = "red", lwd = 1) res <- CEoptim(ssres, f.arg = list(fundf = FN, times = times, y = FitzHugh), continuous = list(mean = c(0, 0, 5, 0, 0), sd = c(1, 1, 1, 1, 1), smoothMean = 0.9, smoothSd = 0.5), verbose = TRUE) res mu <- res$states[, "mean1"] sig <- res$states[, "maxSd"] library("scatterplot3d") t <- length(mu) a <- 15 x <- seq(from = mu[a] - 3 * sig[a], to = mu[a] + 3 * sig[a], by = 6 * sig[a]/100) s3d <- scatterplot3d(x, y = rep(a, 101), z = dnorm(x, mean = mu[a], sd = sig[a]), type = "l", angle = 15, ylim = c(a, t), zlim = c(0, 380), box = F, xlab = "First parameter", ylab = "Iteration", zlab = "Density", mar = c(3, 2.5, 0, 2.5), y.margin.add = 0.5) for (i in (a + 1):t) { x <- seq(from = mu[i] - 3 * sig[i], to = mu[i] + 3 * sig[i], by = 6 * sig[i]/100) s3d$points3d(x, y = rep(i, 101), z = dnorm(x, mean = mu[i], sd = sig[i]), type = "l") } ## 5.4 Max-Cut Problem library("sna") data("lesmis", package = "CEoptim") gplot(lesmis, gmode = "graph") fmaxcut <- function(x, costs) { v1 <- which(x == 1) v2 <- which(x == 0) return(sum(costs[v1, v2])) } set.seed(5) p0 <- list() for (i in 1:77) { p0 <- c(p0, list(rep(0.5, 2))) } p0[[1]] <- c(0, 1) res <- CEoptim(fmaxcut, f.arg = list(costs = lesmis), maximize = TRUE, discrete = list(probs = p0), N = 3000L, verbose = TRUE) res ind <- res$optimizer$discrete group1 <- colnames(lesmis)[which(ind == FALSE)] group2 <- colnames(lesmis)[which(ind == TRUE)] group1 group2 probs <- res$states.probs X <- matrix(NA, nrow = length(probs), ncol = 77) prob0 <- cbind(1, t(rep(0.5, 76))) for (i in 1:length(probs)) { for (j in 1:77) { X[i, j] <- res$states.probs[[i]][[j]][2] } } X <- rbind(prob0, X) par(mfcol = c(5, 2), mar = c(1, 1.5, 1, 1.5), oma = c(1, 1, 1, 1)) for (i in 1:5) { plot(X[i, ], type = "h", lwd = 4, col = "blue", ylim = c(0, 1), xaxt = "n", yaxt = "n", ylab = "", main = paste("t=", i - 1, sep = "")) axis(2, at = 0.5, labels = 0.5) } for (i in 1:5) { plot(X[1 + 4 * i, ], type = "h", lwd = 4, col = "blue", ylim = c(0, 1), xaxt = "n", yaxt = "n", ylab = "", main = paste("t=", 1 + 4 * i, sep = "")) axis(2, at = 0.5, labels = 0.5) } # a<-c() for (i in 1:1000){ # res<-CEoptim(fmaxcut,f.arg=list(costs=lesmis),maximize=T,discrete=list(probs=p0),N=3000L) # a[i]<-res$optimum } repeating the CEoptim 1000 times the following vector of # results a is obtained: a <- c(534, 534, 535, 534, 531, 534, 531, 534, 534, 533, 534, 535, 535, 534, 534, 535, 535, 535, 532, 533, 534, 535, 535, 534, 535, 535, 535, 534, 533, 533, 534, 532, 529, 535, 533, 534, 535, 535, 533, 535, 533, 533, 531, 532, 532, 533, 535, 535, 534, 534, 533, 533, 528, 532, 534, 533, 535, 532, 531, 535, 534, 533, 535, 535, 534, 535, 535, 530, 535, 533, 535, 533, 534, 535, 532, 535, 533, 534, 531, 534, 535, 533, 534, 534, 533, 534, 532, 535, 530, 532, 532, 535, 534, 533, 535, 534, 533, 535, 535, 535, 535, 533, 531, 533, 532, 535, 533, 535, 533, 530, 531, 533, 535, 533, 531, 531, 531, 532, 532, 532, 533, 534, 531, 532, 533, 532, 532, 535, 533, 532, 531, 534, 534, 529, 534, 535, 534, 535, 534, 534, 531, 535, 535, 535, 533, 535, 532, 535, 535, 534, 531, 532, 535, 534, 532, 535, 531, 535, 535, 533, 531, 530, 535, 534, 531, 531, 534, 533, 534, 532, 535, 534, 531, 535, 532, 535, 535, 531, 534, 533, 534, 534, 530, 535, 535, 533, 532, 534, 535, 533, 534, 532, 535, 534, 532, 534, 535, 533, 532, 532, 533, 531, 533, 534, 535, 531, 535, 531, 532, 535, 532, 535, 535, 534, 533, 533, 532, 534, 534, 534, 533, 535, 535, 535, 535, 535, 534, 533, 535, 535, 534, 533, 534, 534, 532, 533, 530, 531, 530, 532, 535, 534, 533, 535, 535, 534, 535, 535, 533, 533, 534, 532, 533, 535, 532, 531, 533, 529, 533, 534, 534, 534, 534, 532, 534, 534, 534, 535, 533, 533, 535, 534, 532, 531, 535, 535, 535, 531, 533, 535, 533, 533, 533, 531, 535, 535, 533, 535, 535, 532, 535, 534, 534, 535, 535, 533, 533, 532, 532, 533, 534, 535, 534, 532, 535, 532, 534, 535, 535, 531, 532, 532, 533, 535, 533, 535, 532, 534, 533, 531, 531, 535, 530, 533, 531, 535, 535, 535, 530, 530, 530, 534, 534, 534, 531, 533, 532, 532, 535, 535, 535, 533, 535, 532, 532, 533, 535, 533, 535, 535, 532, 533, 535, 533, 532, 533, 534, 531, 533, 533, 533, 533, 532, 533, 534, 534, 535, 534, 533, 533, 535, 535, 533, 533, 533, 535, 534, 535, 535, 535, 535, 535, 533, 535, 534, 531, 532, 534, 533, 535, 535, 533, 533, 533, 534, 530, 534, 534, 535, 532, 535, 530, 535, 535, 533, 535, 534, 535, 533, 533, 533, 529, 534, 534, 533, 535, 534, 534, 535, 533, 534, 533, 535, 533, 533, 534, 535, 533, 534, 533, 531, 535, 534, 531, 532, 534, 528, 528, 535, 534, 532, 532, 535, 530, 531, 535, 535, 532, 528, 532, 535, 535, 535, 534, 535, 535, 531, 535, 532, 534, 532, 535, 530, 531, 532, 532, 533, 532, 532, 534, 532, 535, 534, 535, 534, 533, 533, 534, 532, 532, 531, 532, 534, 534, 533, 532, 534, 534, 534, 534, 535, 533, 534, 535, 534, 533, 532, 535, 535, 532, 533, 533, 534, 535, 535, 532, 533, 533, 533, 534, 534, 532, 535, 534, 535, 532, 534, 532, 534, 533, 535, 533, 535, 533, 535, 529, 532, 533, 532, 533, 531, 534, 534, 535, 535, 535, 535, 534, 535, 533, 534, 534, 534, 533, 534, 533, 535, 534, 534, 535, 534, 533, 532, 532, 535, 535, 534, 533, 533, 534, 533, 534, 533, 533, 533, 531, 535, 534, 534, 534, 535, 534, 535, 533, 535, 535, 535, 532, 533, 533, 533, 534, 528, 533, 535, 533, 535, 533, 533, 533, 534, 533, 535, 534, 534, 535, 535, 533, 535, 533, 534, 535, 535, 534, 533, 533, 533, 535, 535, 532, 535, 535, 535, 534, 533, 533, 535, 534, 533, 534, 531, 535, 534, 529, 532, 534, 534, 535, 534, 533, 533, 530, 535, 534, 531, 535, 535, 530, 535, 535, 532, 533, 530, 531, 535, 534, 532, 532, 534, 535, 535, 534, 534, 535, 534, 532, 533, 531, 532, 535, 530, 531, 533, 535, 535, 532, 535, 535, 532, 534, 532, 534, 535, 531, 532, 535, 535, 535, 535, 535, 531, 531, 532, 533, 535, 535, 535, 535, 535, 532, 530, 533, 531, 533, 535, 535, 533, 533, 534, 531, 534, 535, 532, 535, 535, 529, 535, 534, 534, 533, 534, 533, 533, 530, 533, 533, 532, 534, 532, 535, 535, 533, 533, 533, 533, 533, 532, 530, 535, 533, 532, 534, 530, 535, 535, 535, 533, 535, 531, 533, 532, 535, 534, 534, 534, 533, 533, 534, 535, 534, 531, 531, 535, 534, 533, 533, 535, 535, 535, 534, 535, 534, 534, 533, 535, 535, 535, 531, 534, 535, 533, 534, 535, 535, 533, 534, 535, 534, 534, 531, 535, 533, 535, 533, 535, 535, 533, 533, 532, 533, 534, 534, 534, 535, 533, 529, 533, 531, 534, 535, 534, 532, 532, 532, 534, 533, 535, 532, 534, 535, 535, 535, 535, 532, 535, 533, 534, 533, 531, 531, 534, 533, 532, 531, 535, 535, 530, 535, 535, 533, 535, 535, 535, 535, 532, 531, 535, 533, 533, 533, 534, 532, 534, 529, 535, 535, 532, 535, 535, 534, 532, 532, 535, 535, 529, 531, 534, 529, 534, 535, 535, 531, 534, 535, 532, 535, 533, 535, 534, 535, 535, 532, 531, 530, 534, 535, 532, 533, 532, 533, 534, 533, 535, 531, 535, 535, 534, 533, 535, 532, 535, 534, 533, 534, 535, 535, 535, 535, 534, 535, 534, 530, 535, 533, 535, 532, 533, 535, 531, 535, 533, 532, 535, 532, 533, 534, 534, 535, 534, 533, 535, 532, 531, 533, 532, 534, 535, 533, 531, 534, 534, 533, 535, 535, 531, 531, 534, 534, 535, 533, 534, 535, 533, 535, 533, 534, 535, 535, 533, 535, 533, 534, 530, 532, 534, 529, 534, 533, 532, 535, 534, 535, 535, 535, 533, 528, 533, 535, 535, 533, 532, 531, 532, 533, 530, 535, 532, 535, 533, 533, 535, 534, 534, 529, 534, 533, 533, 534, 531, 534, 533, 534, 534, 535, 535, 532, 533, 534, 530, 531) breaks <- 527.5:535.5 histgr <- hist(a, breaks, plot = FALSE) par(mfcol = c(1, 1), mai = c(0.9, 1.2, 0.5, 0.5), xpd = TRUE) plot(histgr, border = "dark blue", col = "light blue", main = "", xlab = "", labels = TRUE, xaxt = "n") axis(1, at = 527.5:534.5, labels = 528:535) ## 5.5. Dirichlet data library("stats") set.seed(12345) a <- 1:5 K <- length(a) - 1 n <- 100 y <- dirichletrnd(a, n) dirichletLoglike <- function(alpha, Y, n, K) { t <- apply(Y, MARGIN = 1, function(y) { sum((alpha[1:K] - 1) * log(y[1:K])) + (alpha[K + 1] - 1) * log(1 - sum(y[1:K])) }) out <- n * (lgamma(sum(alpha)) - sum(lgamma(alpha))) + sum(t) return(out) } mu0 <- rep(0, times = K + 1) sigma0 <- rep(10, times = K + 1) A <- matrix(rep(0, times = 25), nrow = 5) diag(A) <- rep(-1, times = 5) b <- rep(0, times = 5) res <- CEoptim(dirichletLoglike, f.arg = list(Y = y, n = 100, K = 4), maximize = TRUE, continuous = list(mean = mu0, sd = sigma0, conMat = A, conVec = b, smoothSd = 0.5), N = 10000L, verbose = TRUE) res par(mai = c(0.6, 1, 0.5, 0.2), oma = c(0, 0, 0, 1), mfrow = c(1, 1)) plot(res$states[, "iter"], res$states[, "gammat"], type = "s", col = "blue", xlab = "", ylab = "") lines(res$states[, "optimum"], type = "s", col = "red") ## Minka 2000 dirichletFP <- function(data, K) { logpdata <- colMeans(log(data)) afp <- rep(1, times = K) afpold <- -Inf * afp while (sqrt(sum((afp - afpold)^2) > 10^(-12))) { afpold <- afp s <- sum(afpold) for (k in 1:K) { y <- psigamma(s) + logpdata[k] if(y >= -2.22) { ak <- exp(y) + 0.5 } else { ak <- -1/(y-psigamma(1)) } akold <- -Inf while (abs(ak-akold) > 10^(-12)) { akold <- ak ak <- akold - ((psigamma(akold) - y)/psigamma(akold, 1)) } afp[k] <- ak } } return(afp) } est_alpha <- dirichletFP(y, 5) est_likelihood <- dirichletLoglike(est_alpha, y, n, K) est_alpha est_likelihood ## 5.6 Lasso Regression library("glmnet") set.seed(10) n <- 150 p <- 60 beta <- c(runif(10, 0.5, 1), rep(0, 50)) X <- matrix(rnorm(n * p), ncol = 60) Y <- X %*% matrix(beta, ncol = 1) + rnorm(n) res.glmnet <- glmnet(X, Y) ## Find the lambda value to get a model with a sparsity=10 sparsity.10 <- which(res.glmnet$df == 10) (lambda.10 <- res.glmnet$lambda[sparsity.10[1]]) beta.glmnet <- res.glmnet$beta[, sparsity.10[1]] (ind.beta <- which(res.glmnet$beta[, sparsity.10[1]] != 0)) (beta.glmnet.NZ <- res.glmnet$beta[ind.beta, sparsity.10[1]]) RSS.penalized <- function(x, X, Y, lambda) { out <- (1/2) * mean((Y - X %*% matrix(x, ncol = 1, nrow = dim(X)[2], byrow = TRUE))^2) + lambda * sum(abs(x)) return(out) } mu0 <- rep(0, times = p) sigma0 <- rep(5, times = p) N <- 1000 set.seed(1212) res <- CEoptim(RSS.penalized, f.arg = list(X = X, Y = Y, lambda = lambda.10), continuous = list(mean = mu0, sd = sigma0, sdThr = 1e-05), N = N) beta.CEoptim <- res$optimizer$continuous ## Index of the non-zero coefficient (ind.beta.CEoptim.NZ <- which(abs(beta.CEoptim) > 1e-06)) beta.CEoptinm.NZ <- beta.CEoptim[ind.beta.CEoptim.NZ] (compare.beta.NZ <- rbind(beta.glmnet.NZ, beta.CEoptinm.NZ)) (RSS.penalized(beta.CEoptim, X = X, Y = Y, lambda = lambda.10)) (RSS.penalized(beta.glmnet, X = X, Y = Y, lambda = lambda.10)) ## For all sequence of lambda ## RSS.glmnet <- NULL ## RSS.CEoptim <- NULL ## i <- 0 ## for (lambda in res.glmnet$lambda) { ## i <- i+1 ## print(i) ## res <- CEoptim(f = RSS.penalized, ## continuous = list(mean = mu0, sd = sigma0, sdThr = 0.00001), ## discrete = NULL, N = 1000, f.arg = list(X = X, Y = Y, lambda = lambda)) ## RSS.CEoptim <- c(RSS.CEoptim, RSS.penalized(res$optimizer$continuous, ## X = X, Y = Y, lambda = lambda)) ## RSS.glmnet <- c(RSS.glmnet, RSS.penalized(res.glmnet$beta[, i], ## X = X, Y = Y, lambda = lambda)) ## print(c(RSS.CEoptim[i], RSS.glmnet[i])) ## } ## RSS.CEoptim and RSS.glmnet are the result of all sequence of lambda above. RSS.CEoptim <- c(3.25696291498332, 3.23607605242235, 3.19989619332353, 3.14918372040524, 3.08549881263831, 3.01019999577768, 2.92436607357019, 2.83019168887936, 2.72874062989224, 2.6231176322645, 2.51543159005404, 2.40747624843747, 2.30019688664026, 2.19409020461103, 2.09046503612616, 1.99026788362109, 1.89310770671471, 1.79906951118657, 1.70881136353426, 1.62274733628999, 1.54093726271561, 1.46348866030959, 1.39020954968511, 1.32049713132573, 1.2543234554238, 1.19164731970086, 1.13245037503487, 1.07670841032135, 1.02434108679836, 0.975281675443305, 0.929429784779865, 0.886668107545139, 0.846840233662103, 0.809684512285167, 0.775019801350901, 0.742709070973167, 0.71232755012457, 0.684107402714839, 0.657891003234518, 0.633881760661952, 0.610931793468332, 0.590058333520513, 0.570774027296325, 0.553520393763572, 0.536646313623411, 0.521587205971055, 0.507693489877075, 0.496848331808305, 0.483179439066598, 0.473058883655357, 0.463848142902946, 0.453444874856478, 0.445923086774139, 0.43792893948665, 0.43057322385097, 0.425293956206956, 0.41876568073007, 0.413838707118664, 0.408553373195811, 0.404591206901926, 0.399359798342677, 0.395990059714885, 0.392120294864786, 0.389610966771488, 0.386368304457584, 0.384967404861846, 0.381631703128341, 0.379329478544319, 0.377232778186633, 0.376392868089297, 0.373565990873167, 0.372760569600075, 0.370488553814504) RSS.glmnet <- c(3.26181154040858, 3.24556278824333, 3.21454320884016, 3.16144530298561, 3.09677135308556, 3.0231749323876, 2.93532632512027, 2.83960007586826, 2.73674741980047, 2.63018106089729, 2.5215658484197, 2.41337120812514, 2.3052267503836, 2.19839330112557, 2.09416080613105, 1.99362173142844, 1.89594783375495, 1.80147982835252, 1.71086307421586, 1.62452642320072, 1.54272972922794, 1.46542917522047, 1.39209430798758, 1.32215066446162, 1.25587581663187, 1.1929373454032, 1.1335035258305, 1.07753894288917, 1.02500315056843, 0.975812585148281, 0.929901485867709, 0.887048819600468, 0.847162378011546, 0.809975663457998, 0.775258070195694, 0.742863559524035, 0.712581328872051, 0.684348554676498, 0.658102864983643, 0.633742754544896, 0.61118215570154, 0.590321988746143, 0.571047367564942, 0.553268290158364, 0.536899209761347, 0.521842055319522, 0.508000953895088, 0.495269912785656, 0.483572547967427, 0.472822666336802, 0.462957019869853, 0.45392405947405, 0.445646551572839, 0.438089132328549, 0.431151909401976, 0.424802781205431, 0.418970773612479, 0.413640814474774, 0.408788933884092, 0.404351116586549, 0.400282512456757, 0.39656648308618, 0.393170053075825, 0.390065996786484, 0.387229559487527, 0.38462321793441, 0.382251210494903, 0.380088245845547, 0.37811492510563, 0.376312385342856, 0.374666324207335, 0.373164198392223, 0.371793667415879) par(mai = c(1.2, 1.2, 0.6, 0.2), xpd = FALSE) plot((RSS.glmnet - RSS.CEoptim) ~ res.glmnet$lambda, xlab = expression(lambda), type = "l", col = "red", ylab = "RSS(glmnet)-RSS(CEoptim)", cex.lab = 1.2) abline(h = 0, col = "blue") sum(RSS.glmnet < RSS.CEoptim) ## 5.7 AR(1) Model with Regime Switching set.seed(123) sumsqrs <- function(theta, rm1, x) { N <- length(x) #without x[0] r <- 1 + sort(rm1) # internal end points of regimes if (r[1] == r[2]) { # test for dupes -> invalid regime return(Inf) } thetas <- rep(theta, times = c(r, N) - c(1, r + 1) + 1) xhat <- c(0, head(x, -1)) * thetas ## Compute sum of squared errors sum((x - xhat)^2) } ## Read the data from CEoptim package data("yt", package = "CEoptim") xt <- yt - c(0, yt[-300]) A <- rbind(diag(3), -diag(3)) b <- rep(1, 6) res <- CEoptim(sumsqrs, f.arg = list(xt), continuous = list(mean = c(0, 0, 0), sd = rep(1, 3), conMat = A, conVec = b), discrete = list(categories = c(298L, 298L), smoothProb = 0.5), N = 10000, rho = 0.001) res (est.r <- sort(res$optimizer$discrete) + 1) est.theta <- res$optimizer$continuous est.thetas <- rep(est.theta, times = c(est.r, 300) - c(1, est.r + 1) + 1) xfit <- c(0, head(xt, -1)) * est.thetas # prepend x(0)=0, discard x(N) t <- 1:300 ## Plot for Xt par(mai = c(0, 1.2, 0.2, 0.2), mfrow = c(2, 1)) plot(xt ~ t, type = "l", col = "blue", xlab = "", ylab = "xt", cex.lab = 1.2, xaxt = "n") lines(xfit, col = "red") abline(v = c(100, 200)) legend("topleft", legend = c("Observed data", "Fit from CEoptim"), col = c("blue", "red"), lty = c(1, 1)) par(mai = c(0.5, 1.2, 0.2, 0.2)) plot(yt, type = "l", col = "blue", , xlab = "t", ylab = "yt", cex.lab = 1.2) lines(cumsum(xfit), col = "red") abline(v = c(100, 200)) legend("bottomleft", legend = c("Observed data", "Fit from CEoptim"), col = c("blue", "red"), lty = c(1, 1)) dev.new() par(cex.main = 1.5, cex.axis = 1.5, cex.lab = 1.5, mai = c(1.2, 1.2, 1.2, 0.2), mfrow = c(1, 2)) resid <- xfit - xt plot(resid, ylab = " ", xlab = "t", main = "Residuals") qqnorm(resid, ylab = "residuals")