The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
require(copula)
require(grid)
require(lattice)
source(system.file("Rsource", "dnac.R", package="copula"))
set.seed(271)
Note that nacLL()
will use package
partitions
and polynom
.
Consider the following setup.
<- 250
n <- "Gumbel"
family <- c(0.2, 0.4, 0.6)
tau <- getAcop(family)@iTau(tau)
th <- onacopula(family, C(th[1], , list(C(th[2], 1:2), C(th[3], 3:5)))) G3
Sample and compute the log-likelihood.
<- rnacopula(n, G3)
U nacLL(G3, u=U) # log-likelihood at correct parameters
## [1] 450.654
Consider the following setup.
<- 250
n <- "Gumbel"
family <- getAcop(family)
cop. <- c(0.2, 0.4, 0.6)
tau <- cop.@iTau(tau)
th <- onacopula(family, C(th[1], c(1,4),
cop list(C(th[2], 2:3), C(th[3], 5:7))))
Sample and compute the log-likelihood.
<- rnacopula(n, cop)
U nacLL(cop, u=U) # log-likelihood at correct parameters
## [1] 455.4679
Consider the following setup.
<- 250
n <- "Gumbel"
family <- c(0.25, 0.5)
tau <- getAcop(family)@iTau(tau)
th <- onacopula(family, C(th[1], 1, C(th[2], 2:3))) copTrue
Sample and compute the log-likelihood.
<- rnacopula(n, copTrue)
U nacLL(copTrue, u=U) # log-likelihood at correct parameters
## [1] 136.6138
We consider a (1, (2,..,\(d\)))-structure structure here (we choose \(d=3\) here but larger \(d\) are of course possible; plotting can be done as long as we consider two parameters only).
<- "Gumbel" # choose "Clayton" or "Gumbel"
family <- 1 # non-sectorial indices; *Tr stands for the true (nesting structure/model)
compTr <- 2:3 # sectorial indices (for plotting, need 2:d)
scompTr stopifnot(compTr==1) # otherwise, sub (for plotting, see below) is wrong
We first define the negative log-likelihood.
##' Negative Log Likelihood for the two-parameter case
##' C_0({u_j}, C_1({u_k})) where j in 'comp'; k in 'scomp'
<- function(th, u, family, comp, scomp)
nLL2
{stopifnot(length(th) == 2)
if(th[1] > th[2]) # sufficient nesting condition not fulfilled
return(Inf) # for minimization
<- onacopulaL(family, list(th[1], comp, list(list(th[2], scomp))))
cop -nacLL(cop, u=u)
}
<- 100
n <- getAcop(family)
cop. <- c(0.25, 0.5)
tau <- cop.@iTau(tau)) (thTr
## [1] 1.333333 2.000000
<- onacopula(family, C(thTr[1], compTr, C(thTr[2], scompTr))) # copula
cop <- rnacopula(n, cop) # sample
U <- 0.2 # delta{tau} for defining a range of theta's
h <- cop.@iTau(c(tau[1]-h, tau[1]+h))) (th0
## [1] 1.052632 1.818182
<- cop.@iTau(c(tau[2]-h, tau[2]+h))) (th1
## [1] 1.428571 3.333333
<- 20 # number of grid points
m <- expand.grid(th0= seq(th0[1], th0[2], length.out=m),
grid th1= seq(th1[1], th1[2], length.out=m))
<- apply(grid, 1, nLL2,
val.grid u=U, family=family, comp=compTr, scomp=scompTr)
First we determine some plot supplements including the optimum of the negative log-likelihood on the grid.
<- c(th0=thTr[1], th1=thTr[2],
true.val nLL=nLL2(thTr, u=U, family=family, comp=compTr, scomp=scompTr)) # true value
<- which.min(val.grid)
ind <- c(grid[ind,], nLL=val.grid[ind]) # optimum on the grid
opt.grd <- rbind(true.val, opt.grd) # points to add to wireframe plot
pts <- paste("-log-likelihood of a nested", family, "copula") # title
title <- { if(length(scompTr)==2) bquote(italic(u[3]))
mysec else substitute(list(...,italic(u[j])), list(j=max(scompTr))) }
<- substitute(italic(C(bolditalic(u)))==italic(C[0](u[1],C[1](u[2],MSEC))) ~~~~~~
sub italic(n)==N ~~~~~~ tau(theta[0])==TAU0 ~~~~~~ tau(theta[1])==TAU1,
list(MSEC=mysec, N=n, TAU0=tau[1], TAU1=tau[2]))
<- as.expression(sub) # lattice "bug" (only needed by lattice)
sub <- expression(italic(theta[0]))
xlab <- expression(italic(theta[1]))
ylab <- list(as.expression(-log~L *
zlab group("(",italic(theta[0])*"," ~ italic(theta[1])*";"~bolditalic(u),")")),
rot = 90)
<- list(c(expression(group("(",list(theta[0],theta[1]),")")^T),
sTit expression(group("(",list(hat(theta)["0,n"],hat(theta)["1,n"]),")")^T)))
Now consider a wireframe and a level plot.
wireframe(val.grid~grid[,1]*grid[,2], aspect=1, zoom=1.02, xlim=th0, ylim=th1,
zlim= range(val.grid, as.numeric(pts[,3]), finite=TRUE),
xlab=xlab, ylab=ylab, zlab=zlab, main=title, sub=sub, pts=pts,
par.settings=list(standard.theme(color=FALSE),
layout.heights=list(sub=2.4), background=list(col="#ffffff00"),
axis.line=list(col="transparent"), clip=list(panel="off")),
alpha.regions=0.5, scales=list(col=1, arrows=FALSE),
## add wire/points
panel.3d.wireframe = function(x, y, z, xlim, ylim, zlim, xlim.scaled,
ylim.scaled, zlim.scaled, pts, ...) {panel.3dwire(x=x, y=y, z=z, xlim=xlim, ylim=ylim, zlim=zlim,
xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled,
zlim.scaled=zlim.scaled, ...)
panel.3dscatter(x=as.numeric(pts[,1]), y=as.numeric(pts[,2]),
z=as.numeric(pts[,3]), xlim=xlim, ylim=ylim, zlim=zlim,
xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled,
zlim.scaled=zlim.scaled, type="p", col=1,
pch=c(3,4), lex=2, cex=1.4, .scale=TRUE, ...)
},key =
list(x=-0.01, y=1, points=list(pch=c(3,4), col=1, lwd=2, cex=1.4),
text = sTit, padding.text=3, cex=1, align=TRUE, transparent=TRUE))
levelplot(val.grid~grid[,1]*grid[,2], aspect=1, xlab=xlab, ylab=ylab,
par.settings=list(layout.heights=list(main=3, sub=2),
regions=list(col=gray(140:400/400))),
xlim= extendrange(grid[,1], f = 0.04),
ylim= extendrange(grid[,2], f = 0.04),
main=title, sub=sub, pts=pts,
scales=list(alternating=c(1,1), tck=c(1,0)), contour=TRUE,
panel=function(x, y, z, pts, ...){
panel.levelplot(x=x, y=y, z=z, ...)
grid.points(x=pts[1,1], y=pts[1,2], pch=3,
gp=gpar(lwd=2, col="black")) # + true value
grid.points(x=pts[2,1], y=pts[2,2], pch=4,
gp=gpar(lwd=2, col="black")) # x optimum
},key =
list(x=0.18, y=1.09, points=list(pch=c(3,4), col=1, lwd=2, cex=1.4),
columns = 2, text = sTit, align=TRUE, transparent=TRUE))
For illustration purposes, we start not too closely to the optimum.
<- optim(c(1, 3), nLL2,
ropt u=U, family=family, comp=compTr, scomp=scompTr)
Now compare the different results (true values, optimum on the grid,
optimum via optim()
).
rbind(pts, optim=c(ropt$par, ropt$value))
## th0 th1 nLL
## true.val 1.333333 2 -82.54895
## opt.grd 1.374969 1.929825 -82.81026
## optim 1.361803 1.915033 -82.82763
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.