拓端tecdat|在R语言辅导和Stan中估计截断泊松分布

时间：2022-11-11 07:15:47

 

 数据

 

这是一个非常简化的例子。我产生了1,000个计数观察值，平均值为1.3。然后，如果只观察到两个或更高的那个，我将原始分布与我得到的分布进行比较。 

拓端tecdat|在R语言辅导和Stan中估计截断泊松分布


 

 由此代码生成：

# original data:
set.seed(321)
a <- rpois(1000, 1.3)

# truncated version of data:
b <- a[ a > 1]

# graphic:
data_frame(value = c(a, b),
           variable =  (c("Original data", "Truncated so only observations of 2 or more show up"), c(length(a), length(b)))) %>%
  ggplot(aes(x = value)) +
   (binwidth = 1, colour = "white") +
  facet_wrap(~variable, ncol = 1) +
  ggtitle("Comparing full and truncated datasets from a Poisson distribution") +
  labs(y = "Number of observations")

# fitting a model to original works well:
mean(a)
 (a, "Poisson")

# but obviously not if naively done to the truncated version:
mean(b)
fitdistr(b, "Poisson")


估计​​lambda​​​完整数据（​​a​​）的关键参数效果很好，估计值为1.347，刚好超过1.3的真实值的一个标准误差。

最大似然

 

需要​​dpois​​​和​​ppois​​​函数的截断版本并在​​fitdist​​其中使用这些版本。

#-------------MLE fitting in R-------------------
dtruncated_poisson <- function(x, lambda) {
}
ptruncated_poisson <- function(x, lambda) {
}

fitdist(b, "truncated_poisson", start = c(lambda = 0.5))

请注意，要执行此操作，我将下限阈值指定为1.5; 因为所有数据都是整数，这实际上意味着我们只观察2或更多的观察结果。我们还需要为估计值指定一个合理的起始值​​lambda​​- 这样做太远会导致错误。

 

贝叶斯

对于替代贝叶斯方法，Stan可以很容易地将数据和概率分布描述为截断的。除了我​​x​​​在这个程序中调用的原始数据之外，我们需要告诉它有多少观察（​​n​​​），​​lower_limit​​它被截断了，以及表征我们估计的参数的先验分布所需的任何东西。

以下程序的关键部分是：

  • 在​​data​​​块中，指定数据的​​x​​​下限为​​lower_limit​​
  • 在​​model​​​块中，指定​​x​​​通过截断的分布​​T[lower_limit, ]​​ 
data {
  int n;
  int lower_limit;
  int  x[n];
  real lambda_start_mu;
  real lambda_start_sigma;
}

parameters {
  reallambda;
}

model {
  lambda ~ normal(lambda_start_mu, lambda_start_sigma);

  for(i in 1:n){
    x[i] ~ poisson(lambda) T[lower_limit, ];
  }
}

以下是从R向Stan提供数据的方式：

#-------------Calling Stan from R--------------
data <- list(
x = b,
lower_limit = 2,
n = length(),
lambda_start_sigma = 1
)

fit <- stan("0120-trunc.stan", data = data, cores = 4)


plot(fit) + 
  labs(y = "Estimated parameters") +
  theme_minimal(base_family = "myfont")


这为我们提供​​lambda​​​了与该​​fitdistrplus​​方法匹配的后验分布：1.35，标准偏差为0.08。可信区间的图像：

拓端tecdat|在R语言辅导和Stan中估计截断泊松分布


标签：

相关文章