Non central Chi Squared implemented (issue #13).

config-i1 · Aug 31, 2018 · db86e01 · db86e01
1 parent 8882790
commit db86e01
Show file tree

Hide file tree

Showing 5 changed files with 39 additions and 27 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: greybox
 Type: Package
 Title: Toolbox for Model Building and Forecasting
-Version: 0.3.1.41023
-Date: 2018-08-30
+Version: 0.3.1.41024
+Date: 2018-08-31
 Authors@R: person("Ivan", "Svetunkov", email = "ivan@svetunkov.ru", role = c("aut", "cre"),
                   comment="Lecturer at Centre for Marketing Analytics and Forecasting, Lancaster University, UK")
 URL: https://github.com/config-i1/greybox

diff --git a/NEWS b/NEWS
@@ -34,6 +34,8 @@ Changes:
 * dpois and dnbinom distributions in alm. alm() allows producing prediction intervals for both of them. But covariance matrix of parameters for dnbinom might be hard to calculate...
 * The dispersion parameter of dnbinom in alm() is now estimated separately, which now solves a lot of problems.
 * Renamed parameter A into B for alm(). Very serious thing!
+* distribution="dchisq" in alm() now estimates the non-central Chi Squared distribution. The returned scale corresponds to the estimated number of degrees of freedom, while mu is the squared expectation.
+* rmc() now colours the lines depending on the number of groups. If there's only one, then there's one group and the differences are not significant.
 
 Bugfixes:
 * Fixed a bug with the style="line" in rmc(), where the grouping would be wrong in cases, when one method significantly differs from the others.

diff --git a/R/alm.R b/R/alm.R
@@ -327,7 +327,7 @@ alm <- function(formula, data, subset, na.action,
 
     fitter <- function(B, distribution, y, matrixXreg){
         mu[] <- switch(distribution,
-                       "dchisq" =,
+                       "dchisq" = abs(matrixXreg %*% B[-1])^2,
                        "dpois" = exp(matrixXreg %*% B),
                        "dnbinom" = exp(matrixXreg %*% B[-1]),
                        "dnorm" =,
@@ -346,9 +346,9 @@ alm <- function(formula, data, subset, na.action,
                         "dlaplace" = meanFast(abs(y-mu)),
                         "dlogis" = sqrt(meanFast((y-mu)^2) * 3 / pi^2),
                         "ds" = meanFast(sqrt(abs(y-mu))) / 2,
-                        "dchisq" = 2*mu,
+                        "dchisq" = abs(B[1]),
                         "dpois" = mu,
-                        "dnbinom" = B[1],
+                        "dnbinom" = abs(B[1]),
                         "pnorm" = sqrt(meanFast(qnorm((y - pnorm(mu, 0, 1) + 1) / 2, 0, 1)^2)),
                         "plogis" = sqrt(meanFast(log((1 + y * (1 + exp(mu))) / (1 + exp(mu) * (2 - y) - y))^2)) # Here we use the proxy from Svetunkov et al. (2018)
         );
@@ -366,10 +366,9 @@ alm <- function(formula, data, subset, na.action,
                            "dlaplace" = dlaplace(y, mu=fitterReturn$mu, b=fitterReturn$scale, log=TRUE),
                            "dlogis" = dlogis(y, location=fitterReturn$mu, scale=fitterReturn$scale, log=TRUE),
                            "ds" = ds(y, mu=fitterReturn$mu, b=fitterReturn$scale, log=TRUE),
-                           "dchisq" = dchisq(y, df=fitterReturn$mu, log=TRUE),
+                           "dchisq" = dchisq(y, df=fitterReturn$scale, ncp=fitterReturn$mu, log=TRUE),
                            "dpois" = dpois(y, lambda=fitterReturn$mu, log=TRUE),
-                           "dnbinom" = ifelseFast(any(fitterReturn$scale<=0),-1E+300,
-                                                  dnbinom(y, mu=fitterReturn$mu, size=fitterReturn$scale, log=TRUE)),
+                           "dnbinom" = dnbinom(y, mu=fitterReturn$mu, size=fitterReturn$scale, log=TRUE),
                            "pnorm" = c(pnorm(fitterReturn$mu[ot], mean=0, sd=1, log.p=TRUE),
                                        pnorm(fitterReturn$mu[!ot], mean=0, sd=1, lower.tail=FALSE, log.p=TRUE)),
                            "plogis" = c(plogis(fitterReturn$mu[ot], location=0, scale=1, log.p=TRUE),
@@ -387,8 +386,11 @@ alm <- function(formula, data, subset, na.action,
 
     #### Estimate parameters of the model ####
     if(is.null(B)){
-        if(any(distribution==c("dlnorm","dchisq"))){
-            B <- .lm.fit(matrixXreg,log(y))$coefficients
+        if(distribution=="dlnorm"){
+            B <- .lm.fit(matrixXreg,log(y))$coefficients;
+        }
+        else if(distribution=="dchisq"){
+            B <- .lm.fit(matrixXreg,sqrt(y))$coefficients;
         }
         else{
             B <- .lm.fit(matrixXreg,y)$coefficients;
@@ -397,6 +399,9 @@ alm <- function(formula, data, subset, na.action,
         if(distribution=="dnbinom"){
             B <- c(var(y), B);
         }
+        else if(distribution=="dchisq"){
+            B <- c(1, B);
+        }
 
         # Although this is not needed in case of distribution="dnorm", we do that in a way, for the code consistency purposes
         res <- nloptr(B, CF,
@@ -418,6 +423,10 @@ alm <- function(formula, data, subset, na.action,
     if(distribution=="dnbinom"){
         names(B) <- c("size",variablesNames);
     }
+    else if(distribution==c("dchisq")){
+        scale <- abs(B[1]);
+        names(B) <- c("df",variablesNames);
+    }
     else{
         names(B) <- variablesNames;
     }
@@ -432,9 +441,9 @@ alm <- function(formula, data, subset, na.action,
                        "dlaplace" =,
                        "dlogis" =,
                        "ds" =,
-                       "dchisq" =,
                        "dpois" =,
                        "dnbinom" = mu,
+                       "dchisq" = mu + scale,
                        "dlnorm" = exp(mu),
                        "pnorm" = pnorm(mu, mean=0, sd=1),
                        "plogis" = plogis(mu, location=0, scale=1)
@@ -476,7 +485,7 @@ alm <- function(formula, data, subset, na.action,
                                   distribution=distribution, y=y, matrixXreg=matrixXreg);
         }
 
-        if(distribution=="dnbinom"){
+        if(any(distribution==c("dchisq","dnbinom"))){
             vcovMatrix <- vcovMatrix[-1,-1];
         }
 
@@ -555,7 +564,7 @@ alm <- function(formula, data, subset, na.action,
         dataWork <- eval(mf, parent.frame());
     }
 
-    if(distribution=="dnbinom"){
+    if(any(distribution==c("dchisq","dnbinom"))){
         B <- B[-1];
     }
 

diff --git a/R/methods.R b/R/methods.R
@@ -196,7 +196,7 @@ pointLik.alm <- function(object, ...){
                         "ds" = ds(y, mu=mu, b=scale, log=TRUE),
                         "dpois" = dpois(y, lambda=mu, log=TRUE),
                         "dnbinom" = dnbinom(y, mu=mu, size=scale, log=TRUE),
-                        "dchisq" = dchisq(y, df=mu, log=TRUE),
+                        "dchisq" = dchisq(y, df=scale, ncp=mu, log=TRUE),
                         "plogis" = c(plogis(mu[ot], location=0, scale=1, log.p=TRUE),
                                      plogis(mu[!ot], location=0, scale=1, lower.tail=FALSE, log.p=TRUE)),
                         "pnorm" = c(pnorm(mu[ot], mean=0, sd=1, log.p=TRUE),
@@ -479,16 +479,17 @@ predict.alm <- function(object, newdata, interval=c("none", "confidence", "predi
                                      greyboxForecast$mean*(1-2*pnorm(-greyboxForecast$mean/sqrt(greyboxForecast$variance))));
     }
     else if(object$distribution=="dchisq"){
-        greyboxForecast$mean <- exp(greyboxForecast$mean);
-        if(interval!="n"){
-            greyboxForecast$lower <- qchisq(levelLow,greyboxForecast$mean);
-            greyboxForecast$upper <- qchisq(levelUp,greyboxForecast$mean);
+        greyboxForecast$mean <- greyboxForecast$mean^2;
+        if(interval=="p"){
+            greyboxForecast$lower <- qchisq(levelLow,df=object$scale,ncp=greyboxForecast$mean);
+            greyboxForecast$upper <- qchisq(levelUp,df=object$scale,ncp=greyboxForecast$mean);
         }
         else if(interval=="c"){
-            greyboxForecast$lower <- exp(greyboxForecast$lower);
-            greyboxForecast$upper <- exp(greyboxForecast$upper);
+            greyboxForecast$lower <- (greyboxForecast$lower)^2;
+            greyboxForecast$upper <- (greyboxForecast$upper)^2;
         }
-        greyboxForecast$scale <- 2*greyboxForecast$mean;
+        greyboxForecast$mean <- greyboxForecast$mean + object$scale;
+        greyboxForecast$scale <- object$scale;
     }
     else if(object$distribution=="dlnorm"){
         if(interval=="p"){
@@ -1268,7 +1269,7 @@ vcov.alm <- function(object, ...){
         newCall$subset <- object$subset;
         newCall$distribution <- object$distribution;
         newCall$B <- coef(object);
-        if(object$distribution=="dnbinom"){
+        if(any(object$distribution==c("dchisq","dnbinom"))){
             newCall$B <- c(object$scale, newCall$B);
         }
         newCall$vcovProduce <- TRUE;

diff --git a/R/rmc.R b/R/rmc.R
@@ -259,14 +259,14 @@ plot.rmc <- function(x, ...){
     parDefault <- par(no.readonly=TRUE);
     parMar <- parDefault$mar;
 
-    if(x$p.value > 1-x$level){
-        pointCol <- "darkgrey";
-        lineCol <- "grey";
-    }
-    else{
+    if(ncol(x$groups)>1){
         pointCol <- "#0C6385";
         lineCol <- "#0DA0DC";
     }
+    else{
+        pointCol <- "darkgrey";
+        lineCol <- "grey";
+    }
 
     if(("style" %in% argsNames)){
         style <- substr(args$style,1,1);