r - Weibull 分布的更新函数

Question

m(t)Weibull 分布的更新函数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。我的迭代程序正确吗？谢谢。

score 2 · Accepted Answer

我不认为它会起作用。首先，让 Γ(2k+1) 从 m(t) 的分母移到 A _k中。因此，A _k将大致表现为 1/k!。

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

100^千/千！

来自斯特林公式

！~ k ^k , 制定条款

(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)

更新

好的，我在你的问题中计算了 A _k，得到了相同的答案。我想从 m(t) 中估计项 A _k /Γ(2k+1)，我相信它几乎会被 1/k 支配！学期。为此，我制作了另一个数组 k!*A _k /Γ(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)

表明

我得到了和你一样的 A _k值。
B _k确实非常接近1

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

m(t) ~= - Sum(k=1, k=Infinity) (-1) ^k (t ² ) ^k / k！= 1 - Sum(k=0, k=Infinity) (-t ² ) ^k / k！

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

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

结论

近似值为正。毕竟可能会保持正数，A _k /Γ(2k+1) 与 1/k! 有点不同。
我们说的是 1 - exp(-100)，即 1-3.72*10 ^-44！^{我们正在尝试精确地计算它，将 10 100}甚至更高的值相加和相减。即使使用 MPFR，我也不认为这是可能的。

需要另一种方法

score 2 · Accepted Answer

好的，所以我最终走上了一条完全不同的道路。我已经实现了定义更新函数的积分方程的简单离散化：

m(t) = F(t) + integrate (m(t - s)*f(s), s, 0, t)

积分用矩形规则近似。逼近不同 t 值的积分给出了一个线性方程组。我编写了一个函数来生成方程并从中提取系数矩阵。在看了一些例子之后，我猜到了一个直接定义系数的规则，并用它来为一些例子生成解决方案。特别是我尝试了 shape = 2, t = 10，就像在 OP 的例子中一样，step = 0.1（所以 101 个方程）。

我发现该结果与我在一篇论文中找到的近似结果非常吻合（Baxter 等人，在代码中引用）。由于更新函数是事件的预期数量，因此对于大的 t，它大约等于 t/mu，其中 mu 是事件之间的平均时间；这是了解我们是否在附近任何地方的便捷方式。

我正在使用 Maxima ( http://maxima.sourceforge.net )，它对数值的东西效率不高，但它可以很容易地尝试不同的方面。在这一点上，将最终的数字内容移植到另一种语言（如 Python）将很简单。

感谢 OP 提出问题，感谢 S. Pappadeux 进行有见地的讨论。这是我将离散近似值（红色）与大 t 的近似值（蓝色）进行比较的图。尝试了一些不同步长的示例，我发现随着步长变小，值往往会增加一点，所以我认为红线可能有点低，蓝线可能更接近正确。

这是我的千里马代码：

/* discretize weibull renewal function and formulate system of linear equations
 * copyright 2020 by Robert Dodier
 * I release this work under terms of the GNU General Public License
 *
 * This is a program for Maxima, a computer algebra system.
 * http://maxima.sourceforge.net/
 */

"Definition of the renewal function m(t):" $

renewal_eq: m(t) = F(t) + 'integrate (m(t - s)*f(s), s, 0, t);

"Approximate integral equation with rectangle rule:" $

discretize_renewal (delta_t, k) :=
  if equal(k, 0)
    then m(0) = F(0)
    else m(k*delta_t) =   F(k*delta_t)
                        + m(k*delta_t)*f(0)*(delta_t / 2)
                        + sum (m((k - j)*delta_t)*f(j*delta_t)*delta_t, j, 1, k - 1)
                        + m(0)*f(k*delta_t)*(delta_t / 2);

make_eqs (n, delta_t) :=
  makelist (discretize_renewal (delta_t, k), k, 0, n);

make_vars (n, delta_t) :=
  makelist (m(k*delta_t), k, 0, n);

"Discretized integral equation and variables for n = 4, delta_t = 1/2:" $

make_eqs (4, 1/2);
make_vars (4, 1/2);

make_eqs_vars (n, delta_t) :=
  [make_eqs (n, delta_t), make_vars (n, delta_t)];

load (distrib);
subst_pdf_cdf (shape, scale, e) :=
  subst ([f = lambda ([x], pdf_weibull (x, shape, scale)), F = lambda ([x], cdf_weibull (x, shape, scale))], e);

matrix_from (eqs, vars) :=
 (augcoefmatrix (eqs, vars),
  [submatrix (%%, length(%%) + 1), - col (%%, length(%%) + 1)]);

"Subsitute Weibull pdf and cdf for shape = 2 into discretized equation:" $

apply (matrix_from, make_eqs_vars (4, 1/2));
subst_pdf_cdf (2, 1, %);

"Just the  right-hand side matrix:" $

rhs_matrix_from (eqs, vars) :=
 (map (rhs, eqs),
  augcoefmatrix (%%, vars),
  [submatrix (%%, length(%%) + 1), col (%%, length(%%) + 1)]);

"Generate the right-hand side matrix, instead of extracting it from equations:" $

generate_rhs_matrix (n, delta_t) :=
  [delta_t * genmatrix (lambda ([i, j], if i = 1 and j = 1 then 0
                                        elseif j > i then 0
                                        elseif j = i then f(0)/2
                                        elseif j = 1 then f(delta_t*(i - 1))/2
                                        else f(delta_t*(i - j))), n + 1, n + 1),
   transpose (makelist (F(k*delta_t), k, 0, n))];

"Generate numerical right-hand side matrix, skipping over formulas:" $

generate_rhs_matrix_numerical (shape, scale, n, delta_t) :=
  block ([f, F, numer: true], local (f, F),
         f: lambda ([x], pdf_weibull (x, shape, scale)),
         F: lambda ([x], cdf_weibull (x, shape, scale)),
         [genmatrix (lambda ([i, j], delta_t * if i = 1 and j = 1 then 0
                                               elseif j > i then 0
                                               elseif j = i then f(0)/2
                                               elseif j = 1 then f(delta_t*(i - 1))/2
                                               else f(delta_t*(i - j))), n + 1, n + 1),
          transpose (makelist (F(k*delta_t), k, 0, n))]);

"Solve approximate integral equation (shape = 3, t = 1) via LU decomposition:" $

fpprintprec: 4 $
n: 20 $
t: 1;
[AA, bb]: generate_rhs_matrix_numerical (3, 1, n, t/n);
xx_by_lu: linsolve_by_lu (ident(n + 1) - AA, bb, floatfield);

"Iterative solution of approximate integral equation (shape = 3, t = 1):" $

xx: bb;
for i thru 10 do xx: AA . xx + bb;
xx - (AA.xx + bb);
xx_iterative: xx;

"Should find iterative and LU give same result:" $

xx_diff: xx_iterative - xx_by_lu[1];
sqrt (transpose(xx_diff) . xx_diff);

"Try shape = 2, t = 10:" $

n: 100 $
t: 10 $
[AA, bb]: generate_rhs_matrix_numerical (2, 1, n, t/n);
xx_by_lu: linsolve_by_lu (ident(n + 1) - AA, bb, floatfield);

"Baxter, et al., Eq. 3 (for large values of t) compared to discretization:" $
/* L.A. Baxter, E.M. Scheuer, D.J. McConalogue, W.R. Blischke.
 * "On the Tabulation of the Renewal Function,"
 * Econometrics, vol. 24, no. 2 (May 1982).
 * H(t) is their notation for the renewal function.
 */
H(t) := t/mu + sigma^2/(2*mu^2) - 1/2;

tx_points: makelist ([float (k/n*t), xx_by_lu[1][k, 1]], k, 1, n);

plot2d ([H(u), [discrete, tx_points]], [u, 0, t]), mu = mean_weibull(2, 1), sigma = std_weibull(2, 1);

r - Weibull 分布的更新函数

2 回答 2

Related

Reference