data(PlackettPlackettNP) # grab a solution table for a
# Plackett-Plackett composited copula from the copBasic package
# Then create an environment to house the "table."
PlackettPlackett <- new.env()
assign("NeedToCreateForDemo", PlackettPlackettNP, envir=PlackettPlackett)
# 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, but this effort is experimental.
# Now grab the closest Plackett-Plackett solution having the following six
# arbitrary L-comoments. Then simulate 1000 values and plot them to show
# the underlying bivariate distribution.
PPcop <- lcomoms2.ABcop2parameter(solutionenvir=PlackettPlackett,
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)
# A user in encouraged to inspect the contents of PPcop to "assess" the
# solution by a method of L-comoments, we will now proceed with showing the
# copula via a simulation of the fitted version.
para <- list(cop1=PLACKETTcop, cop2=PLACKETTcop, alpha=PPcop$alpha, beta=PPcop$beta,
para1=PPcop$Cop1Thetas, para2=PPcop$Cop2Thetas)
set.seed(1776)
D <- simCOP(n=5000, cop=composite2COP, para=para, col=rgb(0,0,0,0.1), pch=16)
# The sample L-comoments of the fitted Plackett-Plackett may be found by
lmomco::lcomoms2(D, nmom=4) # from the lmomco package, and six sample values shown
T2.12 <- -0.5151547; T2.21 <- -0.5139863
T3.12 <- 0.1502336; T3.21 <- 0.1721355
T4.12 <- -0.0326277; T4.21 <- 0.0233979
PPcop <- lcomoms2.ABcop2parameter(solutionenvir=PlackettPlackett,
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,
para1=PPcop$Cop1Thetas, para2=PPcop$Cop2Thetas)
D <- simCOP(n=5000, cop=composite2COP, para=para, col=rgb(0,0,0,0.1), pch=16)
level.curvesCOP(cop=composite2COP, para=para, delt=.1, ploton=FALSE)
qua.regressCOP.draw(cop=composite2COP, para=para,
ploton=FALSE, f=seq(0.05, 0.95, by=0.05))
qua.regressCOP.draw(cop=composite2COP, para=para, swap=TRUE,
ploton=FALSE, f=seq(0.05, 0.95, by=0.05), col=c(3,2))
diag <- diagCOP(cop=composite2COP, para=para, ploton=FALSE, lwd=4)
image(gridCOP(cop=composite2COP, para=para), col=terrain.colors(20))
# One can inspect alternative solutions like this
# S <- PPcop$solutions$solutions[,1:16]
# 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