{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "| [Table of Contents](#table_of_contents) | [Data and model](#data_and_model) | [Natural estimators](#natural_estimators) | [NN-DOOLSE, MLE](#doolse) | [NN-MDOOLSE, REMLE](#mdoolse) | [References](#references) |" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Authors:** Andrej Gajdoš, Jozef Hanč, Martina Hančová
*[Faculty of Science](https://www.upjs.sk/en/faculty-of-science/?prefferedLang=EN), P. J. Šafárik University in Košice, Slovakia*
email: [andrej.gajdos@student.upjs.sk](mailto:andrej.gajdos@student.upjs.sk)\n", "***\n", "** EBLUP-NE for cyber attacks** " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " R-based computational tools - ** CVXR**\n", "\n", "### Table of Contents \n", "* [Data and model](#data_and_model) - data and model description of empirical data\n", "* [Natural estimators](#natural_estimators) - EBLUPNE based on NE\n", "* [NN-DOOLSE, MLE](#doolse) - EBLUPNE based on nonnegative DOOLSE (same as MLE)\n", "* [NN-MDOOLSE, REMLE](#mdoolse) - EBLUPNE based on nonnegative MDOOLSE (same as REMLE)\n", "* [References](#references)\n", "\n", "**To get back to the contents, use the Home key.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## R packages" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## CVXR: R-Embedded Object-Oriented Modeling Language for Convex Optimization \n", "\n", "* _Purpose_: R package for solving convex optimization tasks\n", "* _Version_: 0.99-5, 2019\n", "* _URL_: https://cvxr.rbind.io/\n", "* _Computational parameters_ of CVXR:\n", "> * *solver* - the convex optimization solver ECOS, OSQP, and SCS chosen according to the given problem\n", " * **OSQP** for convex quadratic problems\n", " * `max_iter` - maximum number of iterations (default: 4000).\n", " * `eps_abs` - absolute accuracy (default: 1e-3).\n", " * `eps_rel` - relative accuracy (default: 1e-4).\n", " * **ECOS** for convex second-order cone problems \n", " * `maxit` - maximum number of iterations (default: 100).\n", " * `abstol` - absolute accuracy (default: 1e-8).\n", " * `reltol` - relative accuracy (default: 1e-8).\n", " * `feastol` - tolerance for feasibility conditions (default: 1e-8).\n", " * `abstol_inacc` - absolute accuracy for inaccurate solution (default: 5e-5).\n", " * `reltol_inacc` - relative accuracy for inaccurate solution (default: 5e-5).\n", " * `feastol_inacc` - tolerance for feasibility condition for inaccurate solution (default: 1e-4).\n", " * **SCS** for large-scale convex cone problems\n", " * `max_iters` - maximum number of iterations (default: 5000).\n", " * `eps` - convergence tolerance (default: 1e-5).\n", " * `alpha` - relaxation parameter (default: 1.5).\n", " * `scale` - factor by which the data is rescaled, only used if `normalize` is TRUE (default: 1.0).\n", " * `normalize` - whether the heuristic data rescaling should be used (default: TRUE)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Important note:* \n", "**After our testing, we found that the standard CVXR package installation, unlike a local installation using Anaconda R distribution or CRAN distribution, did not work in a Binder repository. We are intensively communicating with authors of the package to fix it. At this moment, you have to use your local installation of Jupyter with R kernel or you can see our live Binder Python notebooks using CVXPY with the same results and very similar commands.** " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:03.456000Z", "start_time": "2019-05-19T17:02:01.290Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Warning message:\n", "\"package 'CVXR' was built under R version 3.5.3\"\n", "Attaching package: 'CVXR'\n", "\n", "The following object is masked from 'package:stats':\n", "\n", " power\n", "\n" ] } ], "source": [ "library(CVXR)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "# Data and Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This FDSLRM application describes the real time series data set representing total weekly number of cyber attacks against a honeynet -- an unconventional tool which mimics real systems connected to Internet like business or school computers intranets to study methods, tools and goals of cyber attackers. \n", "\n", "Data, taken from from _Sokol, 2017_ were collected from November 2014 to May 2016 in CZ.NIC honeynet consisting of Kippo honeypots in medium-interaction mode. The number of time series observations is \n", "$n=72$. \n", "\n", "The suitable FDSLRM, after a preliminary logarithmic transformation of \n", "data $Z(t) = \\log X(t)$, is Gaussian orthogonal:\n", "\n", "$ \n", "\\begin{array}{rl}\n", "& Z(t) & \\! = \\! &\\beta_1+\\beta_2\\cos\\left(\\tfrac{2\\pi t\\cdot 3}{72}\\right)+\\beta_3\\sin\\left(\\tfrac{2\\pi t\\cdot 3}{72}\\right)+\\beta_4\\sin\\left(\\tfrac{2\\pi t\\cdot 4}{72}\\right) + \\\\\n", "& & & +Y_1\\sin\\left(\\tfrac{2\\pi\\ t\\cdot 6}{72}\\right)+Y_2\\sin\\left(\\tfrac{2\\pi\\cdot t\\cdot 7}{72}\\right)+w(t), \\, t\\in \\mathbb{N}, \n", "\\end{array}\n", "$ \n", "\n", "where $\\boldsymbol{\\beta}=(\\beta_1,\\,\\beta_2,\\,\\beta_3,\\,\\beta_4)' \\in \\mathbb{R}^4, \\mathbf{Y} = (Y_1,Y_2)' \\sim \\mathcal{N}_2(\\boldsymbol{0}, \\mathrm{D}), w(t) \\sim \\mathcal{iid}\\, \\mathcal{N} (0, \\sigma_0^2), \\boldsymbol{\\nu}= (\\sigma_0^2, \\sigma_1^2, \\sigma_2^2) \\in \\mathbb{R}_{+}^3$.\n", "\n", "We identified the given and most parsimonious structure of the FDSLRM using an iterative process of the model building and selection based on exploratory tools of *spectral analysis* and *residual diagnostics* (for details see our Jupyter notebook `cyberattacks.ipynb`)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:03.535000Z", "start_time": "2019-05-19T17:02:01.294Z" } }, "outputs": [], "source": [ "# data - time series observation\n", "x <- read.csv2('cyberattacks.csv', header = FALSE, sep = \",\")\n", "x <- as.numeric(as.vector(x[-1,]))\n", "x <- log(x)\n", "t <- 1:length(x)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:03.574000Z", "start_time": "2019-05-19T17:02:01.296Z" } }, "outputs": [], "source": [ "# auxilliary functions to create design matrices F, V and the projection matrix M_F\n", "makeF <- function(times, freqs) {\n", "\n", " n <- length(times)\n", " f <- length(freqs)\n", " c1 <- matrix(1, n)\n", " F <- matrix(c1, n, 1)\n", "\n", " for (i in 1:f) {\n", " F <- cbind(F, cos(2 * pi * freqs[i] * times))\n", " F <- cbind(F, sin(2 * pi * freqs[i] * times))\n", " }\n", "\n", " return(F)\n", "}\n", "\n", "makeV <- function(times, freqs) {\n", "\n", " V <- makeF(times, freqs)\n", " V <- V[, -1]\n", "\n", " return(V)\n", "}\n", "\n", "makeM_F <- function(F) {\n", "\n", " n <- nrow(F)\n", " c1 <- rep(1, n)\n", " I <- diag(c1)\n", " M_F <- I - F %*% solve((t(F) %*% F)) %*% t(F)\n", "\n", " return(M_F)\n", "}" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:03.633000Z", "start_time": "2019-05-19T17:02:01.299Z" } }, "outputs": [], "source": [ "# model parameters\n", "n <- 72 \n", "k <- 4 \n", "l <- 2 \n", "\n", "# model - design matrices F, V\n", "F <- makeF(t, c(3/72, 4/72))\n", "F <- F[,-4]\n", "\n", "V <- makeV(t, c(6/72, 7/72))\n", "V <- V[,-c(1,3)]\n", "\n", "# columns vj of V and their squared norm ||vj||^2\n", "nv2 = colSums(V ^ 2)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:03.700000Z", "start_time": "2019-05-19T17:02:01.301Z" } }, "outputs": [], "source": [ "# auxiliary matrices and vectors\n", "\n", "# Gram matrices GF, GV\n", "GF <- t(F) %*% F\n", "GV <- t(V) %*% V\n", "InvGF <- solve(GF)\n", "InvGV <- solve(GV)\n", "\n", "# projectors PF, MF, PV, MV\n", "In <- diag(n)\n", "PF <- F %*% InvGF %*% t(F)\n", "PV <- V %*% InvGV %*% t(V)\n", "MF <- In - PF\n", "MV <- In - PV\n", "\n", "# residuals e\n", "e <- MF %*% x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "# Natural estimators\n", "\n", "## ANALYTICALLY \n", "using formula (4.1) from _Hancova et al 2019_\n", "\n", ">$\n", "\\renewcommand{\\arraystretch}{1.4}\n", "\\breve{\\boldsymbol{\\nu}}(\\mathbf{e}) =\n", "\\begin{pmatrix}\n", "\\tfrac{1}{n-k-l}\\,\\mathbf{e}'\\,\\mathrm{M_V}\\,\\mathbf{e} \\\\\n", "(\\mathbf{e}'\\mathbf{v}_1)^2/||\\mathbf{v}_1||^4 \\\\\n", "\\vdots \\\\\n", "(\\mathbf{e}'\\mathbf{v}_l)^2/||\\mathbf{v}_l||^4\n", "\\end{pmatrix} \n", "$\n", "\n", "## $\\boldsymbol{1^{st}}$ stage of EBLUP-NE " ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:03.734000Z", "start_time": "2019-05-19T17:02:01.304Z" } }, "outputs": [], "source": [ "# auxilliary function to calculate the norm of a vector\n", "norm_vec <- function(x) sqrt(sum(x^2))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:03.782000Z", "start_time": "2019-05-19T17:02:01.307Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 0.05934201 0.02547468 0.01549533\n", "[1] 0.06641189\n" ] } ], "source": [ "# NE according to formula (4.1)\n", "NE0 <- 1/(n-k-l) * t(e) %*% MV %*% (e)\n", "NEj <- (t(e) %*% V) ^ 2 / nv2 ^ 2\n", "NE <- c(NE0, NEj) \n", "print(NE)\n", "print(norm_vec(NE))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## CVXR\n", "\n", "\n", "NE as a convex optimization problem\n", "\n", ">$\n", "\\begin{array}{ll} \n", "\\textit{minimize} & \\quad \n", "f_0(\\boldsymbol{\\nu})=||\\mathbf{e}\\mathbf{e}' - \\mathrm{VDV'}||^2+||\\mathrm{M_V}\\mathbf{e}\\mathbf{e}'\\mathrm{M_V}-\\nu_0\\mathrm{M_F}\\mathrm{M_V}||^2 \\\\[6pt]\n", "\\textit{subject to} & \\quad \\boldsymbol{\\nu} = \\left(\\nu_0, \\ldots, \\nu_l\\right)'\\in [0, \\infty)^{l+1} \n", "\\end{array}\n", "$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:07.271000Z", "start_time": "2019-05-19T17:02:01.309Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "NEcvxr = 0.05934215 0.02547457 0.01549518\n", "norm = 0.06641194" ] } ], "source": [ "# the optimization variable, objective function\n", "v <- Variable(l+1) \n", "fv <- sum_squares(e%*%t(e)-V%*%diag(v[2:(l+1)])%*%t(V)) + sum_squares(MV%*%e%*%t(e)%*%MV-v[1]%*%MF%*%MV)\n", "\n", "# the optimization problem for NE\n", "objective <- Minimize(fv)\n", "constraints <- list(v >= 0)\n", "prob <- Problem(objective, constraints)\n", "\n", "# solve the NE problem\n", "sol <- solve(prob) \n", "cat(\"NEcvxr =\",as.vector(sol$getValue(v)))\n", "cat(\"\\n\")\n", "cat(\"norm =\", norm_vec(as.vector(sol$getValue(v))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## $\\boldsymbol{2^{nd}}$ stage of EBLUP-NE\n", "using formula (3.10) from _Hancova et al 2019_.\n", ">$\n", "\\mathring{\\nu}_j = \\rho_j^2 \\breve{\\nu}_j; j = 0,1 \\ldots, l\\\\\n", "\\rho_0 = 1, \\rho_j = \\dfrac{\\hat{\\nu}_j||\\mathbf{v}_j||^2}{\\hat{\\nu}_0+\\hat{\\nu}_j||\\mathbf{v}_j||^2} \n", "$\n", ">\n", ">where $\\boldsymbol{\\breve{\\nu}}$ are NE, $\\boldsymbol{\\hat{\\nu}}$ are initial estimates for EBLUP-NE" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:07.305000Z", "start_time": "2019-05-19T17:02:01.314Z" } }, "outputs": [], "source": [ "# EBLUP-NE based on formula (3.9)\n", "rho2 <- function(est) {\n", " result <- c(1)\n", " for(j in 2:length(est)) {\n", " result <- c(result, (est[j]*nv2[j-1]/(est[1]+est[j]*nv2[j-1])) ^ 2)\n", " }\n", " return(result)\n", "}" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:07.340000Z", "start_time": "2019-05-19T17:02:01.316Z" } }, "outputs": [], "source": [ "EBLUPNE <- function(est) {\n", " result <- NE * rho2(est)\n", " return(result)\n", " }" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:07.382000Z", "start_time": "2019-05-19T17:02:01.319Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 1.0000000 0.8821446 0.8169426\n" ] } ], "source": [ "# numerical results\n", "print(rho2(NE))" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:07.426000Z", "start_time": "2019-05-19T17:02:01.321Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 0.05934201 0.02247235 0.01265880\n", "[1] 0.06470492\n" ] } ], "source": [ "print(EBLUPNE(NE))\n", "print(norm_vec(EBLUPNE(NE)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "# NN-DOOLSE or MLE\n", "using the result of the KKT algorithm (tab.3, _Hancova et al 2019_) from `PY-estimation-cyberattacks-SciPy.ipynb`\n", "\n", "## $\\boldsymbol{1^{st}}$ stage of EBLUP-NE " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## KKT algorithm\n", "using the the KKT algorithm (tab.3, _Hancova et al 2019_) \n", " \n", "\n", "$~$\n", ">$\n", "\\qquad \\mathbf{q} = \n", "\\left(\\begin{array}{c}\n", "\\mathbf{e}' \\mathbf{e}\\\\\n", "(\\mathbf{e}' \\mathbf{v}_{1})^2 \\\\\n", "\\vdots \\\\\n", "(\\mathbf{e}' \\mathbf{v}_{l})^2\n", "\\end{array}\\right)\n", "$\n", ">\n", "> $\\qquad\\mathrm{G} = \\left(\\begin{array}{ccccc}\n", "\\small\n", "n^* & ||\\mathbf{v}_{1}||^2 & ||\\mathbf{v}_{2}||^2 & \\ldots & ||\\mathbf{v}_{l}||^2 \\\\\n", "||\\mathbf{v}_{1}||^2 & ||\\mathbf{v}_{1}||^4 & 0 & \\ldots & 0 \\\\\n", "||\\mathbf{v}_{2}||^2 & 0 & ||\\mathbf{v}_{2}||^4 & \\ldots & 0 \\\\\n", "\\vdots & \\vdots & \\vdots & \\ldots & \\vdots \\\\\n", "||\\mathbf{v}_{l}||^2 & 0 & 0 & \\ldots & ||\\mathbf{v}_{l}||^4\n", "\\end{array}\\right)\n", "$" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:07.457000Z", "start_time": "2019-05-19T17:02:01.325Z" } }, "outputs": [], "source": [ "NNMDOOLSE_kkt <- function(X, F, V, method = \"NNDOOLSE\") {\n", " n_star <- length(X)\n", " k <- ncol(F)\n", " l <- ncol(V)\n", " if(method == \"NNMDOOLSE\") {\n", " n_star <- n_star - k\n", " }\n", " \n", " MF <- makeM_F(F)\n", " u <- diag(t(V) %*% V)\n", " G <- rbind(c(n_star,u),cbind(u, diag(u^2)))\n", " b_comb <- expand.grid(rep(list(0:1), l))\n", " \n", " eps <- c(MF %*% X)\n", " q <- c(eps %*% eps, (eps %*% V)^2)\n", " K <- G\n", " s <- vector()\n", " \n", " for(i in 1:nrow(b_comb)) {\n", " K_inv <- matrix()\n", " b <- as.vector(unlist(b_comb[i,]))\n", " for(j in 1:length(b)) {\n", " if(b[j] == 0) {\n", " K[1,j+1] <- 0\n", " K[j+1,j+1] <- -1\n", " }\n", " K_inv <- solve(K)\n", " }\n", " beta <- c(K_inv %*% q)\n", " if(all(beta >= 0)) {\n", " s <- beta[1]\n", " for(m in 1:l) {\n", " if(b[m] == 0) {\n", " s <- c(s, 0)\n", " } else {\n", " s <- c(s, beta[m+1])\n", " }\n", " }\n", " break\n", " } else {\n", " K <- G\n", " }\n", " }\n", " return(list(\"estimates\" = s, \"b\" = b))\n", "}" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:07.545000Z", "start_time": "2019-05-19T17:02:01.327Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 0.05595104 0.02392048 0.01394113\n", "[1] 0.06242647\n", "[1] 1 1\n" ] } ], "source": [ "NN_DOOLSE <- NNMDOOLSE_kkt(x, F, V)\n", "print(NN_DOOLSE$estimates)\n", "print(norm_vec(NN_DOOLSE$estimates))\n", "print(NN_DOOLSE$b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## CVXR\n", "\n", "nonnegative DOOLSE as a convex optimization problem\n", "\n", ">$\n", "\\begin{array}{ll} \n", "\\textit{minimize} & f_0(\\boldsymbol{\\nu})=||\\mathbf{e}\\mathbf{e}'-\\Sigma_\\boldsymbol{\\nu}||^2 \\\\[6pt]\n", "\\textit{subject to} & \\boldsymbol{\\nu} = \\left(\\nu_0, \\ldots, \\nu_l\\right)'\\in [0, \\infty)^{l+1} \n", "\\end{array}\n", "$" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:07.584000Z", "start_time": "2019-05-19T17:02:01.330Z" }, "scrolled": true }, "outputs": [], "source": [ "NNDOOLSE_CVXR <- function(X, F, V) {\n", " n <- length(X)\n", " k <- ncol(F)\n", " l <- ncol(V)\n", "\n", " # GF <- t(F) %*% F\n", " # InvGF <- solve(GF)\n", " I <- diag(n)\n", " # PF <- F %*% InvGF %*% t(F)\n", " #MF <- I - PF\n", " MF <- makeM_F(F)\n", " MFV <- MF %*% V\n", "\n", " SX <- MF %*% X %*% t(X) %*% MF\n", " s <- Variable(l+1)\n", " p_obj <- Minimize(sum_squares(SX - (s[1] %*% I) - (V %*% diag(s[2:(l+1)]) %*% t(V))))\n", " constr <- list(s >= 0)\n", " prob <- Problem(p_obj, constr)\n", "\n", " sol <- solve(prob)\n", "\n", " return(as.vector(sol$getValue(s)))\n", "}" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:08.938000Z", "start_time": "2019-05-19T17:02:01.335Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 0.05595096 0.02392044 0.01394109\n", "[1] 0.06242637\n" ] } ], "source": [ "NN_DOOLSEcvxr <- NNDOOLSE_CVXR(x, F, V)\n", "print(NN_DOOLSEcvxr)\n", "print(norm_vec(NN_DOOLSEcvxr))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## CVXR\n", "\n", "using equivalent (RE)MLE convex problem (proposition 5, _Hancova et al 2019_)\n", "\n", "\n", ">$\n", "\\begin{array}{ll} \n", "\\textit{minimize} & \\quad f_0(\\mathbf{d})=-(n^*\\!-l)\\ln d_0 - \\displaystyle\\sum\\limits_{j=1}^{l} \n", "\t\t\\ln(d_0-d_j||\\mathbf{v}_j||^2+d_0\\mathbf{e}'\\mathbf{e}-\\mathbf{e}'\\mathrm{V}\\,\\mathrm{diag}\\{d_j\\}\\mathrm{V}'\\mathbf{e} \\\\[6pt]\n", "\\textit{subject to} & \\quad d_0 > \\max\\{d_j||\\mathbf{v}_j||^2, j = 1, \\ldots, l\\} \\\\\n", " & \\quad d_j \\geq 0, j=1,\\ldots l \\\\\n", " & \\\\\n", "& \\quad\\text{for MLE: } n^* = n, \\text{ for REMLE: } n^* = n-k \\\\\n", "\\textit{back transformation:} & \\quad \\nu_0 = \\dfrac{1}{d_0}, \\nu_j = \\dfrac{d_j}{d_0\\left(d_0 -d_j||\\mathbf{v}_j||^2\\right)}\n", "\\end{array}\n", "$" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:08.973000Z", "start_time": "2019-05-19T17:02:01.339Z" } }, "outputs": [], "source": [ "MLE_CVXR <- function(X, F, V){\n", "\n", " n <- length(X)\n", " l <- ncol(V)\n", "\n", " MF <- makeM_F(F)\n", " GV <- t(V) %*% V\n", " p <- n - l\n", "\n", " e <- as.vector(MF %*% X)\n", " ee <- as.numeric(t(e) %*% e)\n", " eV <- t(e) %*% V\n", " Ve <- t(V) %*% e\n", " d <- Variable(l+1)\n", " logdetS <- p * log(d[1]) + sum(log(d[1] - GV %*% d[2:(l+1)]))\n", " obj <- Maximize(logdetS - ((d[1] * ee) - (eV %*% diag(d[2:(l+1)]) %*% Ve)))\n", " constr <- list(d[2:(l+1)] >= 0, d[1] >= max_entries(GV %*% d[2:(l+1)]))\n", " p_MLE <- Problem(obj, constr)\n", "\n", " sol <- solve(p_MLE)\n", "\n", " s <- 1 / sol$getValue(d)[1]\n", " sj <- vector()\n", " for(j in 2:(l+1)) {\n", " sj <- c(sj, sol$getValue(d)[j]/(sol$getValue(d)[1] * (sol$getValue(d)[1] - sol$getValue(d)[j] * GV[j-1,j-1])))\n", " }\n", " result <- c(s, sj)\n", "\n", " return(as.vector(result))\n", "\n", "}" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:13.110000Z", "start_time": "2019-05-19T17:02:01.342Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 0.05595104 0.02392050 0.01394115\n", "[1] 0.06242648\n" ] } ], "source": [ "MLEcvxr <- MLE_CVXR(x, F, V)\n", "print(MLEcvxr)\n", "print(norm_vec(MLEcvxr))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## $\\boldsymbol{2^{nd}}$ stage of EBLUP-NE" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:13.156000Z", "start_time": "2019-05-19T17:02:01.345Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 1.0000000 0.8817033 0.8094585\n" ] } ], "source": [ "# numerical results\n", "print(rho2(NN_DOOLSE$estimates))" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:13.198000Z", "start_time": "2019-05-19T17:02:01.348Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 0.05934201 0.02246111 0.01254283\n", "[1] 0.06467842\n" ] } ], "source": [ "print(EBLUPNE(NN_DOOLSE$estimates)) \n", "print(norm_vec(EBLUPNE(NN_DOOLSE$estimates)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "# NN-DOOLSE or REMLE\n", "using the KKT algorithm (tab.3, _Hancova et al 2019_)\n", "\n", "## $\\boldsymbol{1^{st}}$ stage of EBLUP-NE " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## KKT algorithm" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:13.237000Z", "start_time": "2019-05-19T17:02:01.351Z" } }, "outputs": [], "source": [ "NN_MDOOLSE <- NNMDOOLSE_kkt(x, F, V, method = \"NNMDOOLSE\")" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:13.279000Z", "start_time": "2019-05-19T17:02:01.354Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 0.05934201 0.02382629 0.01384694\n", "[1] 0.06542862\n", "[1] 1 1\n" ] } ], "source": [ "print(NN_MDOOLSE$estimates)\n", "print(norm_vec(NN_MDOOLSE$estimates))\n", "print(NN_MDOOLSE$b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## CVXR\n", "\n", "nonnegative DOOLSE as a convex optimization problem\n", "\n", ">$\n", "\\begin{array}{ll} \n", "\\textit{minimize} & f_0(\\boldsymbol{\\nu})=||\\mathbf{e}\\mathbf{e}'-\\mathrm{M_F}\\Sigma_\\boldsymbol{\\nu}\\mathrm{M_F}||^2 \\\\[6pt]\n", "\\textit{subject to} & \\boldsymbol{\\nu} = \\left(\\nu_0, \\ldots, \\nu_l\\right)'\\in [0, \\infty)^{l+1} \n", "\\end{array}\n", "$" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:13.315000Z", "start_time": "2019-05-19T17:02:01.356Z" } }, "outputs": [], "source": [ "NNMDOOLSE_CVXR <- function(X, F, V) {\n", " n <- length(X)\n", " k <- ncol(F)\n", " l <- ncol(V)\n", "\n", " # GF <- t(F) %*% F\n", " # InvGF <- solve(GF)\n", " I <- diag(n)\n", " # PF <- F %*% InvGF %*% t(F)\n", " #MF <- I - PF\n", " MF <- makeM_F(F)\n", " MFV <- MF %*% V\n", "\n", " SX <- MF %*% X %*% t(X) %*% MF\n", " s <- Variable(l+1)\n", " p_obj <- Minimize(sum_squares(SX - (s[1] %*% MF) - (MFV %*% diag(s[2:(l+1)]) %*% t(MFV))))\n", " constr <- list(s >= 0)\n", " prob <- Problem(p_obj, constr)\n", "\n", " sol <- solve(prob)\n", "\n", " return(as.vector(sol$getValue(s)))\n", "}" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:14.668000Z", "start_time": "2019-05-19T17:02:01.359Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 0.05934192 0.02382624 0.01384690\n", "[1] 0.06542851\n" ] } ], "source": [ "NN_MDOOLSEcvxr <- NNMDOOLSE_CVXR(x, F, V)\n", "print(NN_MDOOLSEcvxr)\n", "print(norm_vec(NN_MDOOLSEcvxr))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## CVXR\n", "\n", "using equivalent (RE)MLE convex problem (proposition 5, _Hancova et al 2019_)\n", "\n", "\n", ">$\n", "\\begin{array}{ll} \n", "\\textit{minimize} & \\quad f_0(\\mathbf{d})=-(n^*\\!-l)\\ln d_0 - \\displaystyle\\sum\\limits_{j=1}^{l} \n", "\t\t\\ln(d_0-d_j||\\mathbf{v}_j||^2+d_0\\mathbf{e}'\\mathbf{e}-\\mathbf{e}'\\mathrm{V}\\,\\mathrm{diag}\\{d_j\\}\\mathrm{V}'\\mathbf{e} \\\\[6pt]\n", "\\textit{subject to} & \\quad d_0 > \\max\\{d_j||\\mathbf{v}_j||^2, j = 1, \\ldots, l\\} \\\\\n", " & \\quad d_j \\geq 0, j=1,\\ldots l \\\\\n", " & \\\\\n", "& \\quad\\text{for MLE: } n^* = n, \\text{ for REMLE: } n^* = n-k \\\\\n", "\\textit{back transformation:} & \\quad \\nu_0 = \\dfrac{1}{d_0}, \\nu_j = \\dfrac{d_j}{d_0\\left(d_0 -d_j||\\mathbf{v}_j||^2\\right)}\n", "\\end{array}\n", "$" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:14.712000Z", "start_time": "2019-05-19T17:02:01.361Z" } }, "outputs": [], "source": [ "REMLE_CVXR <- function(X, F, V){\n", "\n", " n <- length(X)\n", " k <- ncol(F)\n", " l <- ncol(V)\n", "\n", " MF <- makeM_F(F)\n", " GV <- t(V) %*% V\n", " p <- n - l - k\n", "\n", " e <- as.vector(MF %*% X)\n", " ee <- as.numeric(t(e) %*% e)\n", " eV <- t(e) %*% V\n", " Ve <- t(V) %*% e\n", " d <- Variable(l+1)\n", " logdetS <- p * log(d[1]) + sum(log(d[1] - GV %*% d[2:(l+1)]))\n", " obj <- Maximize(logdetS - ((d[1] * ee) - (eV %*% diag(d[2:(l+1)]) %*% Ve)))\n", " constr <- list(d[2:(l+1)] >= 0, d[1] >= max_entries(GV %*% d[2:(l+1)]))\n", " p_remle <- Problem(obj, constr)\n", "\n", " sol <- solve(p_remle)\n", "\n", " s <- 1 / sol$getValue(d)[1]\n", " sj <- vector()\n", " for(j in 2:(l+1)) {\n", " sj <- c(sj, sol$getValue(d)[j]/(sol$getValue(d)[1] * (sol$getValue(d)[1] - sol$getValue(d)[j] * GV[j-1,j-1])))\n", " }\n", "\n", " result <- c(s, sj)\n", "\n", " return(as.vector(result))\n", "}" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:18.653000Z", "start_time": "2019-05-19T17:02:01.364Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 0.05934201 0.02382629 0.01384695\n", "[1] 0.06542862\n" ] } ], "source": [ "REMLEcvxr <- REMLE_CVXR(x, F, V)\n", "print(REMLEcvxr)\n", "print(norm_vec(REMLEcvxr))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## $\\boldsymbol{2^{nd}}$ stage of EBLUP-NE" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:18.689000Z", "start_time": "2019-05-19T17:02:01.366Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 1.0000000 0.8747730 0.7985572\n" ] } ], "source": [ "# numerical results\n", "print(rho2(NN_MDOOLSE$estimates))" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "ExecuteTime": { "end_time": "2019-05-19T17:02:18.724000Z", "start_time": "2019-05-19T17:02:01.369Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] 0.05934201 0.02228456 0.01237391\n", "[1] 0.06458475\n" ] } ], "source": [ "print(EBLUPNE(NN_MDOOLSE$estimates)) \n", "print(norm_vec(EBLUPNE(NN_MDOOLSE$estimates)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "# References \n", "This notebook belongs to suplementary materials of the paper submitted to Statistical Papers and available at .\n", "\n", "* Hančová, M., Vozáriková, G., Gajdoš, A., Hanč, J. (2019). [Estimating variance components in time series\n", "\tlinear regression models using empirical BLUPs and convex optimization](https://arxiv.org/abs/1905.07771), https://arxiv.org/, 2019.\n", "\n", "### Abstract of the paper\n", "\n", "We propose a two-stage estimation method of variance components in time series models known as FDSLRMs, whose observations can be described by a linear mixed model (LMM). We based estimating variances, fundamental quantities in a time series forecasting approach called kriging, on the empirical (plug-in) best linear unbiased predictions of unobservable random components in FDSLRM. \n", "\n", "The method, providing invariant non-negative quadratic estimators, can be used for any absolutely continuous probability distribution of time series data. As a result of applying the convex optimization and the LMM methodology, we resolved two problems $-$ theoretical existence and equivalence between least squares estimators, non-negative (M)DOOLSE, and maximum likelihood estimators, (RE)MLE, as possible starting points of our method and a \n", "practical lack of computational implementation for FDSLRM. As for computing (RE)MLE in the case of $ n $ observed time series values, we also discovered a new algorithm of order $\\mathcal{O}(n)$, which at the default precision is $10^7$ times more accurate and $n^2$ times faster than the best current Python(or R)-based computational packages, namely CVXPY, CVXR, nlme, sommer and mixed. \n", "\n", "We illustrate our results on three real data sets $-$ electricity consumption, tourism and cyber security $-$ which are easily available, reproducible, sharable and modifiable in the form of interactive Jupyter notebooks." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Sokol P., Gajdoš, A. (2017). [Prediction of Attacks Against Honeynet Based on Time Series Modeling](https://www.springer.com/gp/book/9783319676203). Silhavy, R., Silhavy, P., & Prokopova, Z. (Eds.). (2017). _Applied Computational Intelligence and Mathematical Methods (Vol. 662)_. Cham: Springer International Publishing, pp. 360-371" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| [Table of Contents](#table_of_contents) | [Data and model](#data_and_model) | [Natural estimators](#natural_estimators) | [NN-DOOLSE, MLE](#doolse) | [NN-MDOOLSE, REMLE](#mdoolse) | [References](#references) |" ] } ], "metadata": { "kernelspec": { "display_name": "R", "language": "R", "name": "ir" }, "language_info": { "codemirror_mode": "r", "file_extension": ".r", "mimetype": "text/x-r-source", "name": "R", "pygments_lexer": "r", "version": "3.5.1" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 1 }