[Vinecopula-commits] r122 - pkg/R

[noreply at r-forge.r-project.org]

Author: tnagler
Date: 2015-08-24 12:50:15 +0200 (Mon, 24 Aug 2015)
New Revision: 122

Modified:
   pkg/R/BiCopEst.r
   pkg/R/BiCopPar2Tau.r
   pkg/R/BiCopTau2Par.r
Log:
* faster implementation of par<->tau conversion for Frank copula
* check input whether tau \in (-1,1) (BiCopTau2Par)

Modified: pkg/R/BiCopEst.r
===================================================================
--- pkg/R/BiCopEst.r	2015-08-21 11:40:31 UTC (rev 121)
+++ pkg/R/BiCopEst.r	2015-08-24 10:50:15 UTC (rev 122)
@@ -365,17 +365,10 @@
         a <- -1
         tau <- -tau
     }
-    f <- function(x) {
-        x/(exp(x) - 1)
-    }
-    tauF <- function(x) 1 - 4/x + 4/x^2 * integrate(f,
-                                                    lower = 0 + .Machine$double.eps^0.5,
-                                                    upper = x)$value
-    v <- uniroot(function(x) tau - tauF(x),
-                 lower = 0, 
-                 upper = 500, 
+    v <- uniroot(function(x) tau - (1 - 4/x + 4/x * debye1(x)), 
+                 lower = 0 + .Machine$double.eps^0.5, upper = 5e5, 
                  tol = .Machine$double.eps^0.5)$root
-    return(a * v)
+    return(a*v)
 }
 
 

Modified: pkg/R/BiCopPar2Tau.r
===================================================================
--- pkg/R/BiCopPar2Tau.r	2015-08-21 11:40:31 UTC (rev 121)
+++ pkg/R/BiCopPar2Tau.r	2015-08-24 10:50:15 UTC (rev 122)
@@ -58,14 +58,7 @@
     } else if (family == 4 || family == 14) {
         tau <- 1 - 1/par
     } else if (family == 5) {
-        f <- function(x) x/(exp(x) - 1)
-        fu <- function(x) integrate(f, lower = 0, upper = x)$value
-        fl <- function(x) integrate(f, lower = x, upper = 0)$value
-        if (any(par > 0)) {
-            tau <- 1 - 4/par + 4/par^2 * sapply(par, fu)
-        } else {
-            tau <- 1 - 4/par - 4/par^2 * sapply(par, fl)
-        }
+        tau <- 1 - 4/par + 4/par * debye1(par)
     } else if (family == 6 || family == 16) {
         # tau = 1 + 4/par^2 * integrate(function(x) log(x)*x*(1-x)^(2*(1-par)/par), 0,
         # 1)$value

Modified: pkg/R/BiCopTau2Par.r
===================================================================
--- pkg/R/BiCopTau2Par.r	2015-08-21 11:40:31 UTC (rev 121)
+++ pkg/R/BiCopTau2Par.r	2015-08-24 10:50:15 UTC (rev 122)
@@ -1,4 +1,7 @@
 BiCopTau2Par <- function(family, tau) {
+    ## sanity check
+    if (any(abs(tau) > 0.99999))
+        stop("some tau is too close to -1 or 1")
     
     ## adjust length for input vectors; stop if not matching
     n <- max(length(family), length(tau))