如何在R中绘制阶乘函数

Question

我试图绘制:

使用以下R代码但没有成功:

N= seq(from=150, to=2000)
P=((factorial(60) / factorial(50))*(factorial(N-60) /factorial(N-150))) /(factorial(N) /factorial(N-100))
plot(N,P)

Answer 1

几乎总是,涉及阶乘的概率表达式是"N选择K"计算的一些结果:

但是通过阶乘来计算它是非常低效的,最重要的是,它在数值上并不稳定.看看你的代码使用factorial():你得到了NaN.

在R中,该choose(N, K)函数快速稳定地计算"N选择K".

现在,仔细检查您给定的配方表明它相当于:

choose(N-100, 50) / choose(N, 60)

所以,你可以这样做:

P <- choose(N-100, 50) / choose(N, 60)
plot(N, P, type = "l")

跟进

嗨,这是一个非常有效的功能.但是这个情节的平均值,模式和中位数与我在同一情节的课程材料中的含义不一致？平均值应为727,Mode = 600,中位数= 679 !! 如何从建议的情节中获取这些描述？

我对你的课程资料想要做的事感到困惑.你给出的概率是条件概率P(D | N),即随机变量的概率D.我们画草图P时N.因此,上图不是概率质量函数!然后,我们如何使用它来计算随机变量N???的均值,模式和中位数等统计数据？

好吧无论如何,既然你要求并坚持得到一个答案,让我们假装这是随机变量的概率质量函数N.但是因为它不是真的,sum(P)不是甚至接近1.我们实际上有sum(P) = 3.843678e-12.因此,要将其用作适当的概率质量函数,我们需要首先对其进行标准化.

P <- P / sum(P)

现P总结为1.

为了计算平均值,我们这样做

sum(N * P)
# [1] 726.978

为了计算模式,我们这样做

N[which.max(P)]
# 599

为了计算中位数,我们这样做

N[which(cumsum(P) > 0.5)[1]]
# 679