data(PlackettPlackettABKGtest) # grab a solution table for a
# Plackett-Plackett composited copula from the copBasic package
# Then create an environment to house the "table".
PlackettPlackettABKG <- new.env()
assign("NeedToCreateForDemo", PlackettPlackettABKGtest, envir=PlackettPlackettABKG)
# Now that the table is assigned into the environment, the parameter estimation
# function can be used. In reality a much much larger solution set is needed.
# Assume one had the following six L-comoments, extract a possible solution.
PPcop <- lcomoms2.ABKGcop2parameter(solutionenvir=PlackettPlackettABKG,
T2.12=-0.5059, T2.21=-0.5110,
T3.12= 0.1500, T3.21= 0.1700,
T4.12=-0.0500, T4.21= 0.0329,
uset3err=TRUE, uset4err=TRUE)
# Now take that solution and setup a parameter object.
para <- list(cop1=PLACKETTcop, cop2=PLACKETTcop,
alpha=PPcop$alpha, beta=PPcop$beta, kappa=PPcop$kappa, gamma=PPcop$gamma,
para1=PPcop$Cop1Thetas, para2=PPcop$Cop2Thetas)
set.seed(57) # set for reproducability
# Example Plot Number 1
D <- simCOP(n=2000, cop=composite3COP, para=para, col=rgb(0,0,0,0.1), pch=16)
print(lmomco::lcomoms2(D, nmom=4)) # See the six extacted sample values for this seed.
T2.12 <- -0.4877171; T2.21 <- -0.4907403
T3.12 <- 0.1642508; T3.21 <- 0.1715944
T4.12 <- -0.0560310; T4.21 <- -0.0350028
PPcop <- lcomoms2.ABKGcop2parameter(solutionenvir=PlackettPlackettABKG,
T2.12=T2.12, T2.21=T2.21,
T3.12=T3.12, T3.21=T3.21,
T4.12=T4.12, T4.21=T4.21, uset4err=TRUE)
para <- list(cop1=PLACKETTcop, cop2=PLACKETTcop,
alpha=PPcop$alpha, beta=PPcop$beta, kappa=PPcop$kappa, gamma=PPcop$gamma,
para1=PPcop$Cop1Thetas, para2=PPcop$Cop2Thetas)
# Example Plot Number 2
D <- simCOP(n=1000, cop=composite3COP, para=para, col=rgb(0,0,0,0.1), pch=16)
level.curvesCOP(cop=composite3COP, para=para, delt=0.1, ploton=FALSE)
qua.regressCOP.draw(cop=composite3COP, para=para,
ploton=FALSE, f=c(seq(0.05, 0.95, by=0.05)))
qua.regressCOP.draw(cop=composite3COP, para=para, swap=TRUE,
ploton=FALSE, f=c(seq(0.05, 0.95, by=0.05)), col=c(3,2))
diag <- diagCOP(cop=composite3COP, para=para, ploton=FALSE, lwd=4)
# Compare plots 1 and 2, some generalized consistency between the two is evident.
# One can inspect alternative solutions like this
# S <- PPcop$solutions$solutions[,1:18]
# B <- S[abs(S$t2.12res) < 0.02 & abs(S$t2.21res) < 0.02 &
# abs(S$t3.12res) < 0.02 & abs(S$t3.21res) < 0.02, ]
#print(B)
Run the code above in your browser using DataLab