威布尔分布的更新函数
Renewal Function for Weibull Distribution
威布尔分布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。我的迭代过程是否正确？谢谢。
我认为这行不通。首先，将 Γ(2k+1) 从 m(t) 的分母移动到 Ak。因此，Ak 将大致表现为 1/k!.
在 m(t) 项的提名者中有 t2k，所以粗略地说，你是在用项
计算总和
100k/k!
来自斯特林公式
k！ ~kk, 做条件
(100/k)k
所以是的，他们会开始减少并收敛到某个东西，但是在第 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!
这实际上是 well-known 总和，它等于带有负参数的 exp()（好吧，你必须为 k=0 添加项）
m(t) ~= 1 - exp(-t2)
结论
近似值为正。毕竟可能会保持正数，Ak/Γ(2k+1) 与 1/k!.
有点不同
我们说的是 1 - exp(-100)，即 1-3.72*10-44！我们正在尝试计算它精确地求和和减去 10100 或更高数量级的值。即使有 MPFR，我也不认为这是可能的。
需要另一种方法
好的，所以我最终在这方面走了一条完全不同的路。我实现了定义更新函数的积分方程的简单离散化：
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 的近似值（蓝色）得到的图。尝试了一些不同步长的例子，我看到随着步长变小，值趋于增加一点，所以我认为红线可能有点低，蓝线可能更接近正确。
这是我的 Maxima 代码：
/* 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);
威布尔分布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。我的迭代过程是否正确？谢谢。
我认为这行不通。首先，将 Γ(2k+1) 从 m(t) 的分母移动到 Ak。因此，Ak 将大致表现为 1/k!.
在 m(t) 项的提名者中有 t2k，所以粗略地说，你是在用项
计算总和100k/k!
来自斯特林公式
k！ ~kk, 做条件
(100/k)k
所以是的，他们会开始减少并收敛到某个东西，但是在第 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!
这实际上是 well-known 总和，它等于带有负参数的 exp()（好吧，你必须为 k=0 添加项）
m(t) ~= 1 - exp(-t2)
结论
近似值为正。毕竟可能会保持正数，Ak/Γ(2k+1) 与 1/k!.
有点不同我们说的是 1 - exp(-100)，即 1-3.72*10-44！我们正在尝试计算它精确地求和和减去 10100 或更高数量级的值。即使有 MPFR，我也不认为这是可能的。
需要另一种方法
好的，所以我最终在这方面走了一条完全不同的路。我实现了定义更新函数的积分方程的简单离散化：
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 的近似值（蓝色）得到的图。尝试了一些不同步长的例子，我看到随着步长变小，值趋于增加一点，所以我认为红线可能有点低，蓝线可能更接近正确。
这是我的 Maxima 代码：
/* 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);