Running ‘complex_rlm.R’ [1s/2s]
Running the tests in ‘tests/complex_rlm.R’ failed.
Complete output:
  > ####
  > #### William Ryan
  > #### 12 September, 2021
  > #### Test script for complex robust linear fit (rlm).
  > #### Modified from rlm in MASS.
  > #### Test script based on Wuber's answer on stackexchange https://stats.stackexchange.com/questions/66088/analysis-with-complex-data-anything-different#
  > ####
  > 
  > 
  > library(complexlm)
  Loading required package: MASS
  
  Attaching package: 'complexlm'
  
  The following objects are masked from 'package:MASS':
  
      psi.bisquare, psi.hampel, psi.huber, rlm
  
  The following objects are masked from 'package:stats':
  
      cor, cov, lm, lm.fit, lm.wfit, mad, mahalanobis, var
  
  The following object is masked from 'package:base':
  
      range
  
  > library(MASS, include.only = "mvrnorm")
  > library(ggplot2)
  > # Synthesize data.
  > # (1) the independent variable `w`.
  > #
  > w.max <- 6 # Max extent of the independent values
  > w <- expand.grid(seq(-w.max,w.max, 3), seq(-w.max,w.max, 3))
  > w <- complex(real=w[[1]], imaginary=w[[2]])
  > w <- w[Mod(w) <= w.max] ### This drops any element of w that has a modulus greater than w.max
  > n <- length(w)
  > #
  > # (2) the dependent variable `z`.
  > #
  > beta <- c(-20+5i, complex(argument=2*pi/3, modulus=3/2)) ### The fist is the intercept, the 2nd the slope.
  > print(beta)
  [1] -20.00+5.000000i  -0.75+1.299038i
  > sigma <- 1.5; rho <- 0.7 # Parameters of the error distribution
  > set.seed(17)
  > e <- mvrnorm(n, c(0,0), matrix(c(1,rho,rho,1)*sigma^2, 2)) ### A bunch of random numbers to add errors to the example data.
  > e <- complex(real=e[,1], imaginary=e[,2]) ### Make those random numbers complex.
  > z <- as.vector((desX <- cbind(rep(1,n), w)) %*% beta + e) ### The design matrix desX is defined in this line.
  > z.pure <- as.vector((desX <- cbind(rep(1,n), w)) %*% beta) # z data without the noise.
  > # 
  > # Add m outliers to z. Draw the outliers from a 
  > #
  > z.clean <- z # Preserve outlier free z.
  > m <- 3
  > outlier.pos <- sample(1:n, m) # Positions of outliers in z.
  > outliers <- mvrnorm(m, c(5.8, -7.32),  matrix(c(1,rho,rho,1)*sigma^2, 2)) # The real and imaginary values of the outliers.
  > outliers <- complex(real=outliers[,1], imaginary=outliers[,2]) # Make the outliers comlex numbers.
  > z <- replace(z, outlier.pos, outliers)
  > #
  > # Collect everything into a dataframe.
  > #
  > fitframe <- data.frame(w, z.pure, z.clean, z)
  > #
  > # Whuber's ordinary complex linear fit.
  > #
  > print(beta.hat <- solve(Conj(t(desX)) %*% desX, Conj(t(desX)) %*% z), digits=4)
                   [,1]
  [1,] -13.6046+2.0829i
  [2,]  -0.1656+0.7124i
  > z.whuber <- beta.hat[2] * w + beta.hat[1]
  > fitframe$z.whuber <- z.whuber
  > fitframe$res.whuber <- as.vector(z - desX %*% beta.hat)
  > #
  > # Robust complex linear fit.
  > #
  > print(mytestfit <- rlm(x = w, y = z, interc = TRUE)) # Uses default psi=psi.huber, the Huber objective function.
  Call:
  rlm(x = w, y = z, interc = TRUE)
  Converged in 18 iterations
  
  Coefficients:
            (intercept)                     w 
  -18.9125204+5.294413i  -0.5504508+1.098552i
  
  Degrees of freedom: 13 total; 11 residual
  Scale estimate: 2.26 
  > rbeta.hat <- mytestfit$coefficients
  > fitframe$z.robust <- mytestfit$coefficients[2] * w + mytestfit$coefficients[1]
  > fitframe$res.robust <- mytestfit$residuals
  > # Robust complex linear fit, with Hampel objective function.
  > #
  > print(mytestfitHam <- rlm(x = w, y = z, psi = psi.hampel, interc = TRUE)) # Uses psi=psi.hampel, the Hampel objective function.
  Call:
  rlm(x = w, y = z, psi = psi.hampel, interc = TRUE)
  Converged in 9 iterations
  
  Coefficients:
             (intercept)                      w 
  -13.6739082+2.1406753i  -0.1870532+0.7302799i
  
  Degrees of freedom: 13 total; 11 residual
  Scale estimate: 13.1 
  > rHambeta.hat <- mytestfitHam$coefficients
  > fitframe$z.robustHam <- mytestfitHam$coefficients[2] * w + mytestfitHam$coefficients[1]
  > fitframe$res.robustHam <- mytestfitHam$residuals
  > # Robust complex linear fit, with Tukey's bisquare objective function.
  > #
  > print(mytestfitBi <- rlm(x = w, y = z, psi = psi.bisquare, interc = TRUE)) # Uses psi=psi.bisquare, Tukey's bisquare objective function.
  Call:
  rlm(x = w, y = z, psi = psi.bisquare, interc = TRUE)
  Converged in 6 iterations
  
  Coefficients:
            (intercept)                     w 
  -19.7175660+5.751964i  -0.5898994+1.139501i
  
  Degrees of freedom: 13 total; 11 residual
  Scale estimate: 2.58 
  > rBibeta.hat <- mytestfitBi$coefficients
  > fitframe$z.robustBi <- mytestfitBi$coefficients[2] * w + mytestfitBi$coefficients[1]
  > fitframe$res.robustBi <- mytestfitBi$residuals
  > 
  > library(reshape2)
  Error in library(reshape2) : there is no package called 'reshape2'
  Execution halted