R Under development (unstable) (2024-02-26 r85990 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(mvord)
Loading required package: minqa
Loading required package: BB
Loading required package: ucminf
Loading required package: dfoptim
> 
> #check z2r
> mvord:::check(identical(mvord:::z2r(355),1))
> mvord:::check(identical(mvord:::z2r(0),(exp(0)-1)/(1+exp(0))))
> mvord:::check(identical(mvord:::z2r(2),(exp(4)-1)/(1+exp(4))))
> 
> error.structure <- cor_ar1(~ 1)
> ndim <- 5
> covar_error <- matrix(rep(1,10), ncol = 1)
> 
> attr(error.structure, "ndim") <- ndim
> attr(error.structure, "covariate") <- covar_error
> attr(error.structure, "npar") <- NCOL(attr(error.structure, "covariate"))
> 
> par.sigma <- 2
> sigma <- diag(ndim)
> sigma[lower.tri(sigma)]   <- mvord:::z2r(par.sigma)^sequence((ndim-1):1)
> sigma <- sigma + t(sigma) - diag(diag(sigma))
> ## make list with two elemets:
> ## rVec: for each row the sigma[lower.tri(sigma)]
> ## sdVec: for correlation a vector of ones
> l <- list(rVec = t(sapply(1:10, function(i) sigma[lower.tri(sigma)])),
+ 	 sdVec = rep(1, ndim))
> 
> mvord:::check(identical(mvord:::build_error_struct.cor_ar1(error.structure,
+                                                            par.sigma), l))
> 
> # rho$ncor.levels <- 2
> # rho$npar.cor <- 1
> # par.sigma <- c(2,0)
> #
> # sigma <- diag(rho$ndim)
> # sigma[lower.tri(sigma)]   <- z2r(par.sigma[1])^sequence((rho$ndim-1):1)
> # sigma <- sigma + t(sigma) - diag(diag(sigma))
> #
> # sigma2 <- diag(rho$ndim)
> # sigma2[lower.tri(sigma2)]   <- z2r(par.sigma[2])^sequence((rho$ndim-1):1)
> # sigma2 <- sigma2 + t(sigma2) - diag(diag(sigma2))
> #
> # check(identical(transf.sigmas.corAR1(par.sigma, rho), list(sigma, sigma2)))
> 
> error.structure <- cor_equi(~ 1)
> ndim <- 5
> covar_error <- cbind(1:10,11:20,21:30)
> 
> attr(error.structure, "ndim") <- ndim
> attr(error.structure, "covariate") <- covar_error
> attr(error.structure, "npar") <- NCOL(attr(error.structure, "covariate"))
> 
> 
> par.sigma <- c(1,2,-3)
> sigma <- diag(ndim)
> sigma[lower.tri(sigma)]   <- mvord:::z2r(covar_error %*% par.sigma)
> sigma <- sigma + t(sigma) - diag(diag(sigma))
> sigma
     [,1] [,2] [,3] [,4] [,5]
[1,]    1   -1   -1   -1   -1
[2,]   -1    1   -1   -1   -1
[3,]   -1   -1    1   -1   -1
[4,]   -1   -1   -1    1   -1
[5,]   -1   -1   -1   -1    1
> ## make list with two elemets:
> ## rVec: for each row the sigma[lower.tri(sigma)]
> ## sdVec: for correlation a vector of ones
> l <- list(rVec = t(sapply(1:10, function(i) sigma[lower.tri(sigma)])),
+ 	 sdVec = rep(1, ndim))
> 
> mvord:::check(identical(mvord:::build_error_struct.cor_equi(error.structure, par.sigma,
+                                                             rveclen = ndim*(ndim-1)/2), l))
> 
> 
> par.sigma <- c(0.1,0.2,-0.3)
> sigma <- diag(ndim)
> sigma[lower.tri(sigma)]   <- mvord:::z2r(covar_error %*% par.sigma)
> sigma <- sigma + t(sigma) - diag(diag(sigma))
> l <- list(rVec = t(sapply(1:10, function(i) sigma[lower.tri(sigma)])),
+ 	 sdVec = rep(1, ndim))
> 
> mvord:::check(identical(mvord:::build_error_struct.cor_equi(error.structure, par.sigma,
+                                                             rveclen = ndim*(ndim-1)/2), l))
> 
> 
> # spherical param
> error.structure <- cor_general(~ 1)
> ndim <- 3
> covar_error <- cbind(rep(1, 10))
> 
> attr(error.structure, "ndim") <- ndim
> attr(error.structure, "covariate") <- covar_error
> npar1 <- attr(error.structure, "ndim") * (attr(error.structure, "ndim") - 1)/2
> attr(error.structure, "npar") <- npar1 * NCOL(attr(error.structure, "covariate"))
> 
> par.sigma <- c(0.5, 1, 2)
> l <-  mvord:::build_error_struct.cor_general(error.structure, par.sigma)
> sigma <- diag(ndim)
> sigma[lower.tri(sigma)]   <- l$rVec[1,]
> sigma <- sigma + t(sigma) - diag(diag(sigma))
> 
> mvord:::check(all.equal(mvord:::backtransf_sigmas(sigma), par.sigma))
> 
> #check for cor_general(~f)
> error.structure <- cor_general(~ factor(c(rep("a", 5), rep("b", 5))))
> ndim <- 3
> covar_error <- cbind(c(rep(1, 5), rep(0,5)), c(rep(0, 5), rep(1, 5)))
> 
> attr(error.structure, "ndim") <- ndim
> attr(error.structure, "covariate") <- covar_error
> npar1 <- attr(error.structure, "ndim") * (attr(error.structure, "ndim") - 1)/2
> attr(error.structure, "npar") <- npar1 * NCOL(attr(error.structure, "covariate"))
> 
> par.sigma <- c(0.5, 1, 2, -0.5, -1, -2)
> l <-  mvord:::build_error_struct.cor_general(error.structure, par.sigma)
> sigma1 <- diag(ndim)
> sigma1[lower.tri(sigma1)]   <- l$rVec[1,]
> sigma1 <- sigma1 + t(sigma1) - diag(ndim)
> 
> sigma2 <- diag(ndim)
> sigma2[lower.tri(sigma2)]   <- l$rVec[6,]
> sigma2 <- sigma2 + t(sigma2) - diag(ndim)
> 
> 
> mvord:::check(all.equal(mvord:::backtransf_sigmas(sigma1), par.sigma[1:3]))
> mvord:::check(all.equal(mvord:::backtransf_sigmas(sigma2), par.sigma[4:6]))
> 
> 
> #check for cov_general(~f)
> error.structure <- cov_general(~ factor(c(rep("a", 5), rep("b", 5))))
> ndim <- 3
> covar_error <- cbind(c(rep(1, 5), rep(0,5)), c(rep(0, 5), rep(1, 5)))
> 
> attr(error.structure, "ndim") <- ndim
> attr(error.structure, "covariate") <- covar_error
> npar1 <- attr(error.structure, "ndim") * (attr(error.structure, "ndim") - 1)/2
> attr(error.structure, "npar") <- (npar1 + ndim)* NCOL(attr(error.structure, "covariate"))
> attr(error.structure, "npar.cor") <- npar1* NCOL(attr(error.structure, "covariate"))
> attr(error.structure, "npar.sd") <-  ndim * NCOL(attr(error.structure, "covariate"))
> 
> ## new
> attr(error.structure, "npar.cor") <- ndim * (ndim - 1)/2 * NCOL(covar_error)
> ##
> par.sigma <- c(0.5, 1, 2, -0.5, -1, -2, c(-1, 0.2, 0.5), c(1, -0.2, -0.5))
> l <-  mvord:::build_error_struct.cov_general(error.structure, par.sigma)
> l
$rVec
            [,1]       [,2]       [,3]
 [1,] -0.3752971 -0.6638028 -0.3960711
 [2,] -0.3752971 -0.6638028 -0.3960711
 [3,] -0.3752971 -0.6638028 -0.3960711
 [4,] -0.3752971 -0.6638028 -0.3960711
 [5,] -0.3752971 -0.6638028 -0.3960711
 [6,]  0.3752971  0.6638028  0.8943176
 [7,]  0.3752971  0.6638028  0.8943176
 [8,]  0.3752971  0.6638028  0.8943176
 [9,]  0.3752971  0.6638028  0.8943176
[10,]  0.3752971  0.6638028  0.8943176

$sdVec
           [,1]      [,2]      [,3]
 [1,] 0.3678794 1.2214028 1.6487213
 [2,] 0.3678794 1.2214028 1.6487213
 [3,] 0.3678794 1.2214028 1.6487213
 [4,] 0.3678794 1.2214028 1.6487213
 [5,] 0.3678794 1.2214028 1.6487213
 [6,] 2.7182818 0.8187308 0.6065307
 [7,] 2.7182818 0.8187308 0.6065307
 [8,] 2.7182818 0.8187308 0.6065307
 [9,] 2.7182818 0.8187308 0.6065307
[10,] 2.7182818 0.8187308 0.6065307

> sigma1 <- diag(ndim)
> sigma1[lower.tri(sigma1)]   <- l$rVec[1,]
> sigma1 <- (sigma1 + t(sigma1) - diag(ndim))
> sigma1 <- t(l$sdVec[1, ] * sigma1) * l$sdVec[1, ]
> sigma2 <- diag(ndim)
> sigma2[lower.tri(sigma2)]   <- l$rVec[6,]
> sigma2 <- sigma2 + t(sigma2) - diag(ndim)
> sigma2 <- t(l$sdVec[6, ] * sigma2) * l$sdVec[6, ]
> 
> sigma1
           [,1]       [,2]       [,3]
[1,]  0.1353353 -0.1686319 -0.4026167
[2,] -0.1686319  1.4918247 -0.7975892
[3,] -0.4026167 -0.7975892  2.7182818
> sigma2
          [,1]      [,2]      [,3]
[1,] 7.3890561 0.8352392 1.0944258
[2,] 0.8352392 0.6703200 0.4441050
[3,] 1.0944258 0.4441050 0.3678794
> 
> mvord:::check(all.equal(c(mvord:::backtransf_sigmas(cov2cor(sigma1)), diag(sigma1)),
+ 	c(par.sigma[1:3], exp(2 * par.sigma[7:9]))))
> mvord:::check(all.equal(c(mvord:::backtransf_sigmas(cov2cor(sigma2)), diag(sigma2)),
+ 	c(par.sigma[4:6], exp(2 * par.sigma[10:12]))))
> 
> 
> proc.time()
   user  system elapsed 
   1.29    0.17    1.40