威布尔分布的更新函数

Question

Weibull 分布的更新函数m(t)如下t = 10所示。

我想找到 的值m(t)。我写了下面的r代码来计算m(t)

last_term = NULL
gamma_k = NULL
n = 50
for(k in 1:n){
  gamma_k[k] = gamma(2*k + 1)/factorial(k)
}

for(j in 1: (n-1)){
  prev = gamma_k[n-j]
  last_term[j] = gamma(2*j + 1)/factorial(j)*prev
}

final_term = NULL
find_value = function(n){
  for(i in 2:n){
  final_term[i] = gamma_k[i] - sum(last_term[1:(i-1)])
  }
  return(final_term)
}
all_k = find_value(n)

af_sum = NULL
m_t = function(t){
for(k in 1:n){
af_sum[k] = (-1)^(k-1) * all_k[k] * t^(2*k)/gamma(2*k + 1)
}
  return(sum(na.omit(af_sum)))
}
m_t(20)

输出是m(t) = 2.670408e+93. 我的迭代程序正确吗？谢谢。

Answer 1

我认为这不会起作用。首先，将 \xce\x93(2k+1) 从 m(t) 的分母移至 A _k。因此，A _{k 的}表现大致为 1/k!。

\n

在 m(t) 项的提名者中有 t ^2k，所以粗略地说，您正在计算项的总和

\n

10万^/ k！

\n

由斯特林公式得

\n

k！~ k ^k，制定条件

\n

(100/k) ^k

\n

所以是的，它们将开始减少并收敛到某个值，但在第 100 个任期之后

\n

无论如何，这是代码，你可以尝试改进它，但它在 k~70 处中断

\n

N <- 20\nA <- rep(0, N)\n\n# compute A_k/gamma(2k+1) terms\nps <- 0.0 # previous sum\nA[1] = 1.0\nfor(k in 2:N) {\n    ps <- ps + A[k-1]*gamma(2*(k-1) + 1)/factorial(k-1)\n    A[k] <- 1.0/factorial(k) - ps/gamma(2*k+1)\n}\n\nprint(A)\n\nt <- 10.0\nt2 <- t*t\n\nr <- 0.0\nfor(k in 1:N){\n    r <- r + (-t2)^k*A[k]\n}\n\nprint(-r)\n

\n

更新

\n

好吧，我按照你的问题计算了 A _{k ，得到了相同的答案。}我想从 m(t) 估计项 A _k /\xce\x93(2k+1)，我相信它将几乎由 1/k 主导！学期。为此，我创建了另一个数组 k!*A _k /\xce\x93(2k+1)，它应该接近于 1。

\n

代码

\n

N <- 20\nA <- rep(0.0, N)\n\npsum <- function( pA, k ) {\n    ps <- 0.0\n    if (k >= 2) {\n        jmax <- k - 1\n        for(j in 1:jmax) {\n            ps <- ps + (gamma(2*j+1)/factorial(j))*pA[k-j]\n        }\n    }\n    ps\n}\n\n# compute A_k/gamma(2k+1) terms\nA[1] = gamma(3)\nfor(k in 2:N) {\n    A[k] <- gamma(2*k+1)/factorial(k) - psum(A, k)\n}\n\nprint(A)\n\nB <- rep(0.0, N)\nfor(k in 1:N) {\n    B[k] <- (A[k]/gamma(2*k+1))*factorial(k)\n}\n\nprint(B)\n

\n

表明

\n

\n
1. _{我得到的 A k}值与你相同。
\n
2. B _k确实非常接近1
\n

\n

这意味着项 A _k /\xce\x93(2k+1) 可以替换为 1/k！快速估计我们可能会得到什么（有替换）

\n

m(t) ~= - Sum(k=1, k=无穷大) (-1) ^k (t ² ) ^k / k! = 1 - Sum(k=0, k=无穷大) (-t ² ) ^k / k!

\n

这实际上是众所周知的和，它等于带负参数的 exp()（好吧，你必须为 k=0 添加项）

\n

m(t) ~= 1 - exp(-t ² )

\n

结论

\n

\n

1. 近似值为正。毕竟可能会保持正数，A _k /\xce\x93(2k+1) 与 1/k! 有点不同。

   \n
\n

2. 我们谈论的是 1 - exp(-100)，即 1-3.72*10 ^-44！^{我们尝试精确地计算 10 100}量级的求和和减去值^{甚至更高的值进行精确求和和减去来计算它。即使有了 MPFR，我也不认为这是可能的。}

   \n
\n

\n

需要另一种方法

\n