Home - News ( United States | United Kingdom | Italy | Germany ) - Football scores

Showing content from http://cran.rstudio.com/web/packages/Matrix/../copula/vignettes/logL_visualization.html below:

Log-Likelihood Visualization for Archimedean Copulas

n <- 200
d <- 100
tau <- 0.2
theta <- copJoe@iTau(tau)
cop <- onacopulaL("Joe", list(theta,1:d))
theta

## [1] 1.443824

Here, the three different methods work “the same”:

set.seed(1)
U1 <- rnacopula(n,cop)
enacopula(U1, cop, "mle") # 1.432885 --  fine

## [1] 1.432898

th4 <- 1 + (1:4)/4
mL.tr <- c(-3558.5, -3734.4, -3299.5, -2505.)
mLt1 <- sapply(th4, function(th) mLogL(th, cop@copula, U1, method="log.poly")) # default
mLt2 <- sapply(th4, function(th) mLogL(th, cop@copula, U1, method="log1p"))
mLt3 <- sapply(th4, function(th) mLogL(th, cop@copula, U1, method="poly"))
stopifnot(all.equal(mLt1, mL.tr, tolerance=5e-5),
          all.equal(mLt2, mL.tr, tolerance=5e-5),
          all.equal(mLt3, mL.tr, tolerance=5e-5))
system.time(r1l  <- curveLogL(cop, U1, c(1, 2.5), X=list(method="log.poly")))

##    user  system elapsed 
##   0.499   0.007   0.498

mtext("all three polyJ() methods on top of each other")
system.time({
    r1J <- curveLogL(cop, U1, c(1, 2.5), X=list(method="poly"),
                     add=TRUE, col=adjustcolor("red", .4))
    r1m  <- curveLogL(cop, U1, c(1, 2.5), X=list(method="log1p"),
                      add=TRUE, col=adjustcolor("blue",.5))
})

##    user  system elapsed 
##   0.759   0.000   0.762

U2 <- rnacopula(n,cop)
summary(dCopula(U2, cop)) # => density for the *correct* parameter looks okay

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
##  0.000e+00  4.900e+01  6.430e+02 2.777e+175  1.932e+04 5.553e+177

## hmm: max = 5.5e177
if(doExtras)
    system.time(r2 <- curveLogL(cop, U2, c(1, 2.5)))
stopifnot(all.equal(enacopula(U2, cop, "mle"), 1.43992755, tolerance=1e-5),
          all.equal(mLogL(1.8, cop@copula, U2), -4070.1953,tolerance=1e-5)) # (was -Inf)

U3 <- rnacopula(n,cop)
(th. <- enacopula(U3, cop, "mle")) # 1.4495

## [1] 1.449569

system.time(r3 <- curveLogL(cop, U3, c(1, 2.5)))

##    user  system elapsed 
##   0.355   0.000   0.356

axis(1, at = th., label = quote(hat(theta)))

U4 <- rnacopula(n,cop)
enacopula(U4, cop, "mle") # 1.4519  (prev. was 2.351 : "completely wrong")

## [1] 1.451916

summary(dCopula(U4, cop)) # ok (had one Inf)

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
##  0.000e+00  7.500e+01  9.080e+02 1.981e+259  1.434e+04 3.961e+261

if(doExtras)
    system.time(r4 <- curveLogL(cop, U4, c(1, 2.5)))
mLogL(2.2351, cop@copula, U4)

## [1] -1789.59

mLogL(1.5,    cop@copula, U4)

## [1] -3882.819

mLogL(1.2,    cop@copula, U4)

## [1] -3517.366

if(doExtras) # each curve takes almost 2 sec
    system.time({
        curveLogL(cop, U4, c(1, 1.01))
        curveLogL(cop, U4, c(1, 1.0001))
        curveLogL(cop, U4, c(1, 1.000001))
    })
## --> limit goes *VERY* steeply up to  0
## --> theta 1.164 is about the boundary:
stopifnot(identical(setTheta(cop, 1.164), onacopula(cop@copula, C(1.164, 1:100))),
      all.equal(600.59577,
            cop@copula@dacopula(U4[118,,drop=FALSE],
                    theta=1.164, log = TRUE), tolerance=1e-5)) # was "Inf"

n <- 200
d <- 150
tau <- 0.3
(theta <- copJoe@iTau(tau))

## [1] 1.772108

cop <- onacopulaL("Joe",list(theta,1:d))

set.seed(47)
U. <- rnacopula(n,cop)
enacopula(U., cop, "mle") # 1.784578

## [1] 1.78459

system.time(r. <- curveLogL(cop, U., c(1.1, 3)))

##    user  system elapsed 
##   0.571   0.000   0.574

## => still looks very good

d <- 180
tau <- 0.4
(theta <- copJoe@iTau(tau))

## [1] 2.219066

cop <- onacopulaL("Joe",list(theta,1:d))

U. <- rnacopula(n,cop)
enacopula(U., cop, "mle") # 2.217582

## [1] 2.217593

if(doExtras)
system.time(r. <- curveLogL(cop, U., c(1.1, 4)))
## => still looks very good

RetroSearch is an open source project built by @garambo | Open a GitHub Issue

Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo

HTML: 3.2 | Encoding: UTF-8 | Version: 0.7.4