Weibull分布的更新函数

Question

给出了威布尔分布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

好吧，所以我在这件事上走了一条完全不同的路。我实现了积分方程的一个简单离散化，它定义了更新函数：

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

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

我发现这个结果与我在一篇论文中发现的近似结果非常吻合(Baxter等人，代码中引用)。由于更新函数是预期的事件数，对于大t，它大约等于t/亩，其中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);