这可能没有一个硬性答案，但我想知道在尝试很好地估计广义帕累托分布的参数时是否需要收集更多数据？

我问的原因是因为我正在尝试使用贝叶斯估计来估计广义帕累托分布的参数，当我有大量数据（例如，1000 多个数据点）时，我的参数估计似乎非常好，但是当我放弃时数据大小说 100 那么估计可能很差。

例如，如果我有一个具有真实参数的广义帕累托分布$\mu=0$,$\sigma=1.2$， 和$\xi=0.8$, 我采样$N=1000$然后观察（运行我的贝叶斯算法）我得到估计$\hat\sigma=1.27$95% 置信区间：$(1.12, 1.46)$和$\hat\xi=0.83$95% 置信区间：$(0.72, 0.98)$.

但是，如果我将相同的大小放到$N=100$我明白了$\hat\sigma=0.87$95% 置信区间：$(0.55, 1.24)$和$\hat\xi=0.94$95% 置信区间：$(0.61, 1.38)$.

如果我继续减少$N$, 的点估计$\sigma$和$\xi$只会变得更糟。是否有关于极值分布需要多少数据的经验法则？在大多数情况下，100 个数据点足以对没有太极端值（例如，正态、指数、伽马等）的分布进行建模。在我的应用程序中，我将始终处理少于 100 个数据点，因此使用广义帕累托分布是一个坏主意吗？

以下是我试图解释的代码示例：

# log-likelihood
likelihood <- function(x, xi, sigma){
  
  llik <- -log(sigma) - (1 / xi + 1) * log(1 + xi * x / sigma)
  lik  <- sum(llik)
  
  return(lik)
}

# log(prior)
prior <- function(xi, sigma){
  
  prior1 <- dgamma(xi, .01, .01, log = TRUE)
  prior2 <- dgamma(sigma, .01, .01, log = TRUE)
  prior <- (prior1 + prior2)
  
  return(prior)
}

# log(posterior)
posterior <- function(x, xi, sigma){
  
  post <- likelihood(x, xi, sigma) + prior(xi, sigma)
  
  return(post)
}


##############################################################
### Function to simulate data from GPD
##############################################################

gpd <- function(n, mu, sigma, xi){
  
  u <- runif(n)
  x = mu + sigma * (u^-xi - 1) / xi
  
  return(x)
}

set.seed(4)
N = 1000 # Number of data points
x = gpd(N, 0, 1.2, .8)  # Here mu = 0, sigma = 1.2, and xi = 0.8



S <- 10000
xi <- rep(NA, S)
sigma <- rep(NA, S)
xi[1] <- 1
sigma[1] <- 1

for(i in 2:S){
  
  
  # MCMC for xi
  xi.star = xi[i-1] + rnorm(1,0)
  
  if(xi.star < 0){
    alpha = 0
  }else{
    ratio <- exp(posterior(x, xi.star, sigma[i-1]) - posterior(x, xi[i-1], sigma[i-1]))
    alpha <- min(1, ratio)
  }

  if(runif(1) < alpha){
    xi[i] <- xi.star
  }else{
    xi[i] <- xi[i - 1]
  }
  
  
  # MCMC for sigma
  sigma.star = xi[i-1] + rnorm(1,0)
  
  if(sigma.star < 0){
    alpha = 0
  }else{
    ratio <- exp(posterior(x, xi[i-1], sigma.star) - posterior(x, xi[i-1], sigma[i-1]))
    alpha <- min(1, ratio)
  }
  
  if(runif(1) < alpha){
    sigma[i] <- sigma.star
  }else{
    sigma[i] <- sigma[i - 1]
  }  
}


sigma <- sigma[5000:S]
xi <- xi[5000:S]

mean(sigma)
mean(xi)