威 bool 分布的更新函数

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

我不认为它会起作用。首先，让 Γ(2k+1) 从 m(t) 的分母移到 Ak 中。因此，Ak 的表现大致为 1/k!。
在 m(t) 项的提名者中，有 t2k，所以粗略地说，您正在用项计算总和
10 万/千!
从斯特林公式
克! ~ kk, 制定条件
(100/千)千
所以是的，他们会开始减少并收敛到某个东西，但在第 100 个学期之后
无论如何，这是代码，您可以尝试改进它，但它在 k~70 处中断

N <- 20
A <- rep(0, N)

# compute A_k/gamma(2k+1) terms
ps <- 0.0 # previous sum
A[1] = 1.0
for(k in 2:N) {
    ps <- ps + A[k-1]*gamma(2*(k-1) + 1)/factorial(k-1)
    A[k] <- 1.0/factorial(k) - ps/gamma(2*k+1)
}

print(A)

t <- 10.0
t2 <- t*t

r <- 0.0
for(k in 1:N){
    r <- r + (-t2)^k*A[k]
}

print(-r)

更新
好的，我按照你的问题计算了 Ak，得到了相同的答案。我想从 m(t) 估计 Ak/Γ(2k+1) 项，我相信它几乎会被 1/k 所支配!学期。为此，我制作了另一个数组 k!*Ak/Γ(2k+1)，它应该接近于 1。
代码

N <- 20
A <- rep(0.0, N)

psum <- function( pA, k ) {
    ps <- 0.0
    if (k >= 2) {
        jmax <- k - 1
        for(j in 1:jmax) {
            ps <- ps + (gamma(2*j+1)/factorial(j))*pA[k-j]
        }
    }
    ps
}

# compute A_k/gamma(2k+1) terms
A[1] = gamma(3)
for(k in 2:N) {
    A[k] <- gamma(2*k+1)/factorial(k) - psum(A, k)
}

print(A)

B <- rep(0.0, N)
for(k in 1:N) {
    B[k] <- (A[k]/gamma(2*k+1))*factorial(k)
}

print(B)

表明

我得到了和你一样的 Ak 值。

Bk 确实非常接近 1

这意味着术语 Ak/Γ(2k+1) 可以替换为 1/k!快速估计我们可能会得到什么(替换)
m(t) ~= - Sum(k=1, k=Infinity) (-1)k (t2)k/k! = 1 - Sum(k=0, k=Infinity) (-t2)k/k!
这实际上是众所周知的总和，它等于带有负参数的 exp() (好吧，您必须为 k=0 添加项)
m(t) ~= 1 - exp(-t2)
结论

近似值为正。毕竟可能会保持正数，Ak/Γ(2k+1) 与 1/k! 有点不同。

我们说的是 1 - exp(-100)，也就是 1-3.72*10-44!我们正试图计算它精确地求和和减去 10100 或更高数量级的值。即使使用 MPFR，我也不认为这是可能的。

需要另一种方法