Function for estimating Generalized Pareto Distribution parameters in R not working-CodePudding

In this paper https://www.jstor.org/stable/27867255?seq=1#metadata_info_tab_contents they provide a method of estimation for the GPD, based on the profile likelihood function and the empirical Bayesian method. The code they provide is the following:

# x is the sample data from the GPD
f <- function(x) {
  n <- length(x)
  x <- sort(x)
  lx <- function(b,x) {
    k <- -mean(log(1-b*x))
    if (b==0) {
      k-1-log(mean(x))
    } else {
      k-1 log(b/k)
    }
    }
  p <- (3:9)/10
  xp <- x[round(n*(1-p) .5)]
  m <- 20 round(n^.5)
  xq <- x[round(n*(1-p*p) .5)]
  k <-  log(xq/xp-1,p)
  a <- k*xp/(1-p^k)
  a[k==0] <- (-xp/log(p))[k==0]
  k <- -1
  b <- w <- L <- (n-1)/(n 1)/n[n]-(1-((1:m-.5)/m)^k)/k/median(a)/2
  for (i in 1:m) L[i] <- n*lx(b[i],x)
  for (i in 1:m) w[i] <- 1/sum(exp(L-L[i]))
  b <- sum(b*w)
  k <- -mean(log(1-b*x))
  list(sigma=k/b,k=k)
}

I simulated a vector of x values with size n = 100 from a GPD with shape parameter equals to 1 and scale parameter equals to -1. Then I get the error:

Error in if (condition) { : missing value where TRUE/FALSE needed

According to this question Error in if/while (condition) {: missing Value where TRUE/FALSE needed the error can be attributed to the NA values in the b vector, but I don't understand why this happens.

My x vector:

x
  [1]   0.24264350   0.71418670   0.90929131   1.31011008   1.26467953   0.66491141  40.16132902
  [8]   0.24815930   0.83783287   1.18879459   1.87167421   0.00396890   1.58125086   0.32042532
 [15]   0.01330135   0.80800732   0.67832942   0.32623323  10.87447448   0.64779973   0.16550198
 [22]   1.03595024   0.08463867   6.17820208  15.07862358   0.19188437   0.14925543   2.20031552
 [29]   1.14093185   2.92184174   0.79279984   0.14736789   0.32954571  24.41728358   5.64925376
 [36]   1.31141081   1.62897827   0.01491289   0.86641233   0.82289179   0.12105522   1.66701079
 [43]   4.65817711  53.08764134   8.87696704   0.61560327   0.77817268   0.65148331   2.33976096
 [50]   2.96629395   0.43783850   3.11005777   0.44923740   0.27073261   2.98854135   1.85962571
 [57]   0.67065561   1.09661475   1.55934896   1.65683579   0.62373160  27.04125805   0.16137916
 [64]   0.32925148   0.40487288   1.41714764   0.48798225   0.07624247   2.32993211   2.13228723
 [71]   6.74938943   2.10121108   0.56772145   2.84504482   1.67119601   1.68640938  23.65633183
 [78]   7.78195561   0.21317910   0.33639542   1.14508402   6.61466064   2.03818446   2.28166528
 [85]   0.28095778  30.93008603   0.16512528   1.21974281   0.32121843 103.90277963   0.09672460
 [92]   1.04882224   1.68444513  29.93956683   0.43961534   2.60623811   0.25039076   0.76208631
 [99] 568.46662349   1.66133511

CodePudding user response：

I think you have some transcription errors in your function (to be fair, the code is not well formatted in the paper). I think this should be what you're looking for:

f <- function(x) {
  n <- length(x) 
  x <- sort(x)
 lx <- function(b, x) {
   k <- -mean(log(1 - b * x))
   if (b == 0) k - 1 - log(mean(x)) else k - 1   log(b/k) 
 }
 p <- (3:9)/10 
 xp <- x[round(n*(1 - p)  0.5)] 
 m <- 20   round(n^0.5)
 xq <- x[round(n*(1 - p * p)   .5)]
 k <- log(xq / xp - 1, p)
 a <- k*xp/(1 - p^k)
 a[k == 0] <- (-xp/log(p))[k == 0]
 k <- -1
 b <- w <- L <- (n-1)/(n 1)/x[n] -(1-((1:m-.5)/m)^k)/k/median(a)/2
 for (i in 1:m) L[i] <- n*lx(b[i],x) 
 for (i in 1:m) w[i] <- 1/sum(exp(L-L[i])) 
 b <- sum(b * w) 
 k <- -mean(log(1 - b * x))
 list(sigma = k / b, k = k)
}

Which gives, on your data:

f(x)
#> $sigma
#> [1] 1.050348
#> 
#> $k
#> [1] -1.070265

data

x <- c(0.2426435, 0.7141867, 0.90929131, 1.31011008, 1.26467953, 0.66491141, 
40.16132902, 0.2481593, 0.83783287, 1.18879459, 1.87167421, 0.0039689, 
1.58125086, 0.32042532, 0.01330135, 0.80800732, 0.67832942, 0.32623323, 
10.87447448, 0.64779973, 0.16550198, 1.03595024, 0.08463867, 
6.17820208, 15.07862358, 0.19188437, 0.14925543, 2.20031552, 
1.14093185, 2.92184174, 0.79279984, 0.14736789, 0.32954571, 24.41728358, 
5.64925376, 1.31141081, 1.62897827, 0.01491289, 0.86641233, 0.82289179, 
0.12105522, 1.66701079, 4.65817711, 53.08764134, 8.87696704, 
0.61560327, 0.77817268, 0.65148331, 2.33976096, 2.96629395, 0.4378385, 
3.11005777, 0.4492374, 0.27073261, 2.98854135, 1.85962571, 0.67065561, 
1.09661475, 1.55934896, 1.65683579, 0.6237316, 27.04125805, 0.16137916, 
0.32925148, 0.40487288, 1.41714764, 0.48798225, 0.07624247, 2.32993211, 
2.13228723, 6.74938943, 2.10121108, 0.56772145, 2.84504482, 1.67119601, 
1.68640938, 23.65633183, 7.78195561, 0.2131791, 0.33639542, 1.14508402, 
6.61466064, 2.03818446, 2.28166528, 0.28095778, 30.93008603, 
0.16512528, 1.21974281, 0.32121843, 103.90277963, 0.0967246, 
1.04882224, 1.68444513, 29.93956683, 0.43961534, 2.60623811, 
0.25039076, 0.76208631, 568.46662349, 1.66133511)