In pderdeyn/Stats230VineCopulaClustering: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(echo = TRUE)

remotes::install_github("oezgesahin/vineclust")
library(vineclust)

install.packages("../",repos=NULL,type="source")
#devtools::install_github('https://github.com/pderdeyn/Stats230pieter')
library(Stats230VineCopulaClustering)

Simulate data

# Simulation setup given in Section 5.2 of the paper at https://arxiv.org/pdf/2102.03257.pdf
dims <- 3
obs <- c(500,500) 
RVMs <- list()
M1<-matrix(c(1,3,2,0,3,2,0,0,2),dims,dims)
F1<-matrix(c(0,3,4,0,0,14,0,0,0),dims,dims)
cop1<-matrix(c(0,0.8571429,2.5,0,0,5,0,0,0),dims,dims)

RVMs[[1]] <- VineCopula::RVineMatrix(Matrix=M1,
                        family=F1,
                        par=cop1,
                        par2=matrix(sample(0, dims*dims, replace=TRUE),dims,dims)) 
M2<-matrix(c(1,3,2,0,3,2,0,0,2), dims,dims)
F2<-matrix(c(0,6,5,0,0,13,0,0,0), dims,dims)
cop2<-matrix(c(0,1.443813,11.43621,0,0,2,0,0,0),dims,dims)


RVMs[[2]] <- VineCopula::RVineMatrix(Matrix=M2,
                        family=F2,
                        par=cop2,
                        par2=matrix(sample(0, dims*dims, replace=TRUE),dims,dims))
margin <- matrix(c('Skew Student-t', 'Normal', 'Logistic', 'Skew Student-t', 'Normal', 'Logistic'), 3, 2) 

margin_pars <- array(0, dim=c(4, 3, 2))
margin_pars[,1,1] <- c(1, 2, 7, 1.5)
margin_pars[,2,2] <- c(1.5, 0.4, 0, 0)
margin_pars[,2,1] <- c(1, 0.2, 0, 0)
margin_pars[,1,2] <- c(2, 2, 9, 1)
margin_pars[,3,1] <- c(0, 1, 0, 0)
margin_pars[,3,2] <- c(0.5, 0.5, 2.5, 0)
x_data <- rvcmm(dims, obs, margin, margin_pars, RVMs)

library(GGally)
cols <- character(nrow(x_data))
cols[1:obs[1]] <- "vine 1"
cols[(obs[1]+1):(obs[1]+obs[2])] <- "vine 2"
ggpairs(x_data,columns=1:3, aes(color = cols,alpha=0.35),upper="blank")

mix<-c(0.5,0.5)
vcmm_result<-test_vineclust(x_data, margin_pars, M1, M2, F1, F2, cop1, cop2, mix, iter=5)

library(KScorrect)
lower<-0.5
upper<-2
noise <- array(rlunif(prod(dim(margin_pars)),lower,upper),dim(margin_pars))
margin_parsp <- noise * margin_pars

noise <- array(rlunif(prod(dim(cop1)),lower,upper),dim(cop1))
cop1p <- noise * cop1

noise <- array(rlunif(prod(dim(cop2)),lower,upper),dim(cop2))
cop2p <- noise * cop2

noise <- rlunif(length(mix),lower,upper)
mixp <- noise * mix

vcmm_result<-test_vineclust(x_data, margin_parsp, M1, M2, F1, F2, cop1p, cop2p, mixp,iter=5)

r_perturbed <- fit_points_to_vcmm(x_data,margin_parsp,M1,M2,F1,F2,cop1p,cop2p,mixp)

labels_p <- apply(r_perturbed,1,which.max)
labels_p <- paste(labels_p)
ggpairs(x_data,columns=1:3, aes(color = labels_p,alpha=0.35),upper="blank")

plot(vcmm_result[[1]],ylab="log likelihood")

gam1fit<-vcmm_result[[2]]
gam2fit<-vcmm_result[[3]]
margin_parsfit<-array(c(gam1fit,gam2fit),dim(margin_pars))
cop1fit<-vcmm_result[[4]]
cop2fit<-vcmm_result[[5]]
mixfit<-vcmm_result[[6]]

r_fit <- fit_points_to_vcmm(x_data,margin_parsfit,M1,M2,F1,F2,cop1fit,cop2fit,mixfit)

labels_p <- apply(r_fit,1,which.max)
labels_p <- paste(labels_p)
ggpairs(x_data,columns=1:3, aes(color = labels_p,alpha=0.35),upper="blank")

```

pderdeyn/Stats230VineCopulaClustering documentation built on March 21, 2022, 12:56 a.m.