2

作为一项任务，我必须使用 PMF 为给定的几何分布开发和算法并生成样本

在此处输入图像描述

使用逆变换方法，我想出了以下用于生成值的表达式：

在此处输入图像描述

其中 U 代表一个值或 n 个值，具体取决于样本的大小，取自 Unif(0,1) 分布，p 为 0.3，如上述 PMF 所述。

我有算法，R 中的实现，我已经生成了 QQ 图，以直观地评估经验值对理论值的调整（用 R 生成），即，如果生成的样本确实遵循几何分布。

现在我想将生成的样本提交给拟合优度测试，即卡方，但我在 R 中遇到了麻烦。

4

3 回答 3

2

[尽管您对 whuber 的问题做出了回应，但我认为这有点仓促，因为我认为在解决“我如何在 R 中编写此算法”问题之前，处理“您是什么”可能更重要做不是解决问题的最佳方法（当然属于您发布它的地方）。既然它在这里，我将处理“在 R 中做”方面，但我会敦促你回过头来询问第二个问题（作为新帖子）。]

首先，卡方检验略有不同，具体取决于您是否测试

H0：数据来自参数为p的几何分布

或者

H0：数据来自参数为 0.3 的几何分布

如果你想要第二个，这很简单。首先，对于几何，如果要对检验统计量的分布使用卡方近似，则需要对尾部的相邻单元格进行分组。“通常”的规则 - 过于保守 - 表明您需要每个 bin 中的预期计数至少为 5。

我假设你有一个很好的大样本量。在这种情况下，您将拥有许多具有大量预期计数的垃圾箱，您无需担心保持如此高的数量，但您仍然需要选择如何将尾部装箱（无论您是否只选择一个例如，所有值都分组的截止值）。

我会继续说 n 是 1000 （尽管如果你正在测试你的几何随机数生成，那是相当低的）。

首先，计算您的预期计数：

 dgeom(0:20,.3)*1000
 [1] 300.0000000 210.0000000 147.0000000 102.9000000  72.0300000  50.4210000
 [7]  35.2947000  24.7062900  17.2944030  12.1060821   8.4742575   5.9319802
[13]   4.1523862   2.9066703   2.0346692   1.4242685   0.9969879   0.6978915
[19]   0.4885241   0.3419669   0.2393768

警告，dgeom朋友们从 x=0 开始，而不是 x=1；虽然您可以将输入和输出转移到 R 函数，但如果您从所有几何值中减去 1 并对其进行测试，则会容易得多。我将继续进行，就好像您的样本已减去 1，以便它从 0 开始。

我将在第 15 个学期 (x=14) 将其切断，并将 15+ 分组到它自己的组中（在这种情况下为单个组）。如果你想遵循“大于五”的经验法则，你会在第 12 个学期 (x=11) 之后将其切断。在某些情况下（例如较小的 p），您可能希望将尾部拆分为多个 bin 而不是一个。

> expec <- dgeom(0:14,.3)*1000
> expec <- c(expec, 1000-sum(expec))
> expec
 [1] 300.000000 210.000000 147.000000 102.900000  72.030000  50.421000
 [7]  35.294700  24.706290  17.294403  12.106082   8.474257   5.931980
[13]   4.152386   2.906670   2.034669   4.747562

最后一个单元格是“15+”​​类别。我们还需要概率。

现在我们还没有样品；我将只生成一个：

y <- rgeom(1000,0.3)

但现在我们想要一个观察计数表：

 (x <- table(factor(y,levels=0:14),exclude=NULL))

   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14 <NA> 
 292  203  150   96   79   59   47   25   16   10    6    7    0    2    5    3

现在您可以直接计算卡方，然后计算 p 值：

> (chisqstat <- sum((x-expec)^2/expec))
[1] 17.76835
(pval <- pchisq(chisqstat,15,lower.tail=FALSE))
[1] 0.2750401

但你也可以让 R 来做：

> chisq.test(x,p=expec/1000)

        Chi-squared test for given probabilities

data:  x 
X-squared = 17.7683, df = 15, p-value = 0.275

Warning message:
In chisq.test(x, p = expec/1000) :
  Chi-squared approximation may be incorrect

现在未指定 p 的情况是相似的，但是（据我所知）你不能再chisq.test直接做，你必须用第一种方法做，但你必须从数据中估计参数（通过最大似然或最小卡方），然后按上面的方法进行测试，但估计参数的自由度要少一个。

请参阅此处为具有估计参数的泊松做卡方的示例；几何遵循与上述大致相同的方法，在链接处进行调整（处理未知参数，包括 1 个自由度的损失）。

于 2013-12-03T19:13:02.593 回答
1

让我们假设您在 vector 中获得了随机生成的变量x。您可以执行以下操作：

x <- rgeom(1000,0.2)

x_tbl <- table(x)
x_val <- as.numeric(names(x_tbl))
x_df <- data.frame(count=as.numeric(x_tbl), value=x_val)

# Expand to fill in "gaps" in the values caused by 0 counts
all_x_val <- data.frame(value = 0:max(x_val))
x_df <- merge(all_x_val, x_df, by="value", all.x=TRUE)
x_df$count[is.na(x_df$count)] <- 0

# Get theoretical probabilities 
x_df$eprob <- dgeom(x_df$val, 0.2)

# Chi-square test: once with asymptotic dist'n, 
# once with bootstrap evaluation of chi-sq test statistic
chisq.test(x=x_df$count, p=x_df$eprob, rescale.p=TRUE)
chisq.test(x=x_df$count, p=x_df$eprob, rescale.p=TRUE, 
   simulate.p.value=TRUE, B=10000)
于 2013-12-03T18:57:10.023 回答
1

在包“vcd”中有一个“goodfit”函数，描述为“离散数据的拟合优度测试”。

G.fit <- goodfit(x, type = "nbinomial", par = list(size = 1))

我打算使用您在之前的问题中发布的代码，但现在看来您已经删除了该代码。我觉得这很冒犯。您是否正在使用此论坛收集作业答案，然后对其进行破坏以删除证据？（我们这些有足够代表的人仍然可以看到已删除的问题，并且该界面会阻止删除带有赞成答案的问题，因此您应该无法删除此问题。）

生成用于测试几何分布样本的 QQ 图

- - 问题 - -

我有一个在 R 中生成的 n 个元素的样本

sim.geometric <- function(nvals)
{
    p <- 0.3
    u <- runif(nvals)
    ceiling(log(u)/log(1-p))
}

我想测试它的分布，特别是它是否确实遵循几何分布。我想生成一个 QQ PLot 但不知道如何生成。

--------重新发布的答案----------

与从具有相同概率参数的几何分布中抽取的“真实”样本相比，QQ 图应该是一条直线。一个给函数两个向量，这些向量本质上是在每个分位数比较它们的逆 ECDF。（你的尝试并不是特别成功：）

sim.res <- sim.geometric(100) sim.rgeom <- rgeom(100, 0.3) qqplot(sim.res, sim.rgeom)

在这里，我跟随 qqplot 帮助页面的作者的引导（这导致围绕身份线翻转上曲线）：

png("QQ.png")
qqplot(qgeom(ppoints(100),prob=0.3), sim.res,
       main = expression("Q-Q plot for" ~~ {G}[n == 100]))
dev.off()

---图片不包括---

您可以通过绘制一条穿过每个分布的第 25 个和第 75 个百分位点的线来添加“良好拟合线”。（我为此添加了一个抖动功能，以便更好地了解“概率质量”的位置：）

sim.res <- sim.geometric(500)
qqplot(jitter(qgeom(ppoints(500),prob=0.3)), jitter(sim.res),
       main = expression("Q-Q plot for" ~~ {G}[n == 100]), ylim=c(0,max( qgeom(ppoints(500),prob=0.3),sim.res )),
xlim=c(0,max( qgeom(ppoints(500),prob=0.3),sim.res )))
 qqline(sim.res, distribution = function(p) qgeom(p, 0.3),
       prob = c(0.25, 0.75), col = "red")
于 2013-12-04T07:20:09.500 回答