R version 4.5.0 RC (2025-04-04 r88113 ucrt) -- "How About a Twenty-Six"
Copyright (C) 2025 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(PCDimension)
Loading required package: ClassDiscovery
Loading required package: cluster
Loading required package: oompaBase
> # simulate one unstructured dataset
> suppressWarnings( RNGversion("3.5.0") )
> set.seed(992500)
> NC <- 15
> NS <- 200
> ranData <- matrix(rnorm(NS*NC, 6), ncol=NC)
> # perform PCA and get the set of variances/eigenvectors
> spca <- SamplePCA(ranData)
> lam <- spca@variances[1:(NC-1)] # remove 0 eigenvalue
> rg <- 0:(NC-2)
> 
> ######################################################
> # Test some internal routines
> 
> # Theory says that fhat is non-decreasing
> fh <- sapply(rg, PCDimension:::fhat, Lambda=lam)
> all( diff(fh) >= 0 ) # should be 'TRUE'
[1] TRUE
> 
> plot( rg, fh, type='b', pch=16 )
> 
> # Theory says that dhat is non-increasing
> thetas <- seq(0, 0.2, by=0.01)
> dset <- sapply(thetas, PCDimension:::dhat, Lambda=lam)
> all( diff(dset) <= 0 ) # should be 'TRUE'
[1] TRUE
> 
> plot(thetas, dset, type='l')
> 
> ######################################################
> # test external interface on random data
> 
> # Auer-Gervini Bayesian method
> obj <- AuerGervini(spca)
> summary(obj)
An 'AuerGervini' object that estimates the number of principal components to be 0.
> agDimension(obj)
[1] 0
> 
> plot(obj)
> 
> # broken-stick model
> bs <- brokenStick(1:NC, NC)
> bs0 <- brokenStick(1:(NC-1), (NC-1))
> 
> pts <- screeplot(spca, ylim=c(0, 0.2))
> lines(pts, bs, type='b', col='blue', lwd=2, pch=16)
> lines(pts[-NC], bs0, type='b', col='red', lwd=2, pch=16)
> 
> # randomization based model
> rndLambdaF(ranData)
rndLambda      rndF 
        0         0 
> 
> rm(spca)
> ######################################################
> # get the previously computed simulated SamplePCA object
> data(spca)
> 
> vars <- spca@variances
> obj <- AuerGervini(vars, dd=dim(spca@scores))
> summary(obj)
An 'AuerGervini' object that estimates the number of principal components to be 1.
> 
> plot(obj)
> 
> f <- makeAgCpmFun("Exponential")
> agfuns <- list(twice=agDimTwiceMean,
+                km=agDimKmeans, km3=agDimKmeans3, 
+                tt=agDimTtest, tt2=agDimTtest2,
+                cpt=agDimCPT, cpm=f)
> plot(obj, agfuns)
> 
> compareAgDimMethods(obj, agfuns)
twice    km   km3    tt   tt2   cpt   cpm 
    1     1     1     3     1     1     3 
> 
> proc.time()
   user  system elapsed 
   2.20    0.32    2.51