paulnorthrop / threshr
Showing 1 of 1 files from the diff.

@@ -42,7 +42,7 @@
Loading
42 42
43 43
# =========================== gp_obs_info ===========================
44 44
45 -
gp_obs_info <- function(gp_pars, y) {
45 +
gp_obs_info <- function(gp_pars, y, eps = 1e-5) {
46 46
  # Observed information for the generalized Pareto distribution
47 47
  #
48 48
  # Calculates the observed information matrix for a random sample \code{y}
@@ -51,18 +51,46 @@
Loading
51 51
  #
52 52
  # Args:
53 53
  #   gp_pars : A numeric vector. Parameters sigma and xi of the
54 -
  #   generalized Pareto distribution.
54 +
  #             generalized Pareto distribution.
55 55
  #   y       : A numeric vector. A sample of positive data values.
56 +
  #   eps     : A (small, positive) numeric scalar.  If abs(xi) is smaller than
57 +
  #             eps then we approximate the [2, 2] element of the information
58 +
  #             matrix using a Taylor series approximation.
56 59
  #
57 60
  # Returns:
58 61
  #   A 2 by 2 numeric matrix.  The observed information matrix.
59 62
  #
63 +
  if (eps <= 0) {
64 +
    stop("'eps' must be positive")
65 +
  }
66 +
  # sigma
60 67
  s <- gp_pars[1]
68 +
  # xi
61 69
  x <- gp_pars[2]
62 70
  i <- matrix(NA, 2, 2)
63 -
  i[1,1] <- -sum((1 - (1 + x) * y * (2 * s + x * y) / (s + x * y) ^ 2) / s ^ 2)
64 -
  i[1,2] <- i[2,1] <- -sum(y * (1 - y / s) / (1 + x * y / s) ^ 2 / s ^ 2)
65 -
  i[2,2] <- sum(2 * log(1 + x * y / s) / x ^ 3 - 2 * y / (s + x * y) / x ^ 2 -
66 -
                  (1 + 1 / x) * y ^ 2 / (s + x * y) ^ 2)
71 +
  i[1, 1] <- -sum((1 - (1 + x) * y * (2 * s + x * y) / (s + x * y) ^ 2) / s ^ 2)
72 +
  i[1, 2] <- i[2, 1] <- -sum(y * (1 - y / s) / (1 + x * y / s) ^ 2 / s ^ 2)
73 +
  # Direct calculation of i22 is unreliable for x close to zero.
74 +
  # If abs(x) < eps then we expand the problematic terms (all but t4 below)
75 +
  # in powers of z up to z ^ 2. The terms in 1/z and 1/z^2 cancel leaving
76 +
  # only a quadratic in z.
77 +
  z <- x / s
78 +
  t0 <- 1 + z * y
79 +
  if (any(t0 <= 0)) {
80 +
    stop("The log-likelihood is 0 for this combination of data and parameters")
81 +
  }
82 +
  if (abs(x) < eps) {
83 +
    s1 <- 12 * z ^ 2 * y ^ 2 / 5
84 +
    s2 <- 3 * z * y / 2
85 +
    s3 <- 2 / 3
86 +
    i[2, 2] <- sum(y ^ 3 * (s1 - s2 + s3) / s ^ 3 - (y / s) ^ 2 / t0 ^ 2)
87 +
  } else {
88 +
    t1 <- 2 * log(t0) / z ^ 3
89 +
    t2 <- 2 * y / (z ^ 2 * t0)
90 +
    t3 <- y ^ 2 / (z * t0 ^ 2)
91 +
    t4 <- y ^ 2 / t0 ^ 2
92 +
    i[2, 2] <- sum((t1 - t2 - t3) / s ^ 3 - t4 / s ^ 2)
93 +
  }
94 +
  dimnames(i) <- list(c("sigma[u]", "xi"), c("sigma[u]", "xi"))
67 95
  return(i)
68 96
}
Files Coverage
R 65.94%
Project Totals (14 files) 65.94%
cl5bf4voumh6d9jb

No yaml found.

Create your codecov.yml to customize your Codecov experience

Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading