这是一个非常简化的例子。我产生了1,000个计数观察值,平均值为1.3。然后,如果只观察到两个或更高的那个,我将原始分布与我得到的分布进行比较。
由此代码生成:
<- rpois(1000, 1.3)\n\n# truncated version of data:\nb <- a[ a > 1]\n\n# graphic:\ndata_frame(value = c(a, b),\n variable = (c(\"Original data\", \"Truncated so only observations of 2 or more show up\"), c(length(a), length(b)))) %>%\n ggplot(aes(x = value)) +\n (binwidth = 1, colour = \"white\") +\n facet_wrap(~variable, ncol = 1) +\n ggtitle(\"Comparing full and truncated datasets from a Poisson distribution\") +\n labs(y = \"Number of observations\")\n\n# fitting a model to original works well:\nmean(a)\n (a, \"Poisson\")\n\n# but obviously not if naively done to the truncated version:\nmean(b)\nfitdistr(b, \"Poisson\")","classes":[],"lang":""}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet"> # 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提供数据的方式:
<- list(\n x = b,\n lower_limit = 2,\n n = length(),\n lambda_start_sigma = 1\n)\n\nfit <- stan(\"0120-trunc.stan\", data = data, cores = 4)\n\n\nplot(fit) + \n labs(y = \"Estimated parameters\") +\n theme_minimal(base_family = \"myfont\")","classes":[],"lang":""}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet">#-------------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。可信区间的图像:
联系客服