# ============================================================ # Age-replacement for discrete PH (Eryilmaz, 2021) — General R code # Inputs: # pi : 1 x m row vector (initial phase probabilities) # Q : m x m transient transition matrix (spectral radius < 1) # c1 : failure replacement cost (c1 > c2) # c2 : preventive replacement cost # Nmax: upper cap for search (default 500) # Outputs: # list(N_star, C_star, details = list(S, p, r, L, lambdas, omegas)) # ============================================================ ph_age_replacement <- function(pi, Q, c1, c2, Nmax = 500, tol = 1e-10) { # ---------- basic dimensions and helpers ---------- pi <- as.numeric(pi) m <- length(pi) e <- rep(1, m) I <- diag(1, m) u <- as.numeric((I - Q) %*% e) # exit probabilities (row-sum deficits) # sanity checks if (abs(sum(pi) - 1) > 1e-8) warning("sum(pi) != 1") if (any(u < -1e-10)) warning("Some exit probs negative? Check Q.") if (c1 <= c2) stop("Need c1 > c2 by assumption.") # ---------- E[T] ---------- E_T <- as.numeric(pi %*% solve(I - Q) %*% e) # ---------- p(n), S(n), r(n) up to Nmax ---------- # definitions: # p(n) = pi Q^{n-1} u', n >= 1 # S(n) = P{T >= n} = pi Q^{n-1} e', n >= 1 ; and S(0)=1 # r(n) = p(n)/S(n), n >= 1 p <- numeric(Nmax) S <- numeric(Nmax + 1) r <- numeric(Nmax) S[1] <- 1.0 # S(0) = 1 Qn <- diag(1, m) for (n in 1:Nmax) { Qn <- Qn %*% Q S[n + 1] <- as.numeric(pi %*% Qn %*% e) # S(n) p[n] <- as.numeric(pi %*% (Qn / Q) %*% u) # WRONG if written like this! } # Correct p(n): need Q^{n-1}; we already have Qn = Q^n now. So compute separately: Qn <- diag(1, m) for (n in 1:Nmax) { if (n == 1) { Qpow <- diag(1, m) # Q^{0} } else { Qpow <- Qpow %*% Q # Q^{n-1} } p[n] <- as.numeric(pi %*% Qpow %*% u) } # r(n) for n>=1 uses S(n)=S[n] where S[n] stores S(n-1) actually; we defined S[1]=S(0) # So S(n) is S[n+1] for (n in 1:Nmax) { denom <- S[n + 1] r[n] <- ifelse(denom > tol, p[n] / denom, NA_real_) } # ---------- eigen-decomposition (lambdas) ---------- # According to the algorithm text: roots of lambda^m + b1 lambda^{m-1} + ... + bm = 0 # are the eigenvalues of Q (under the mapping z=1/lambda). So: ev <- eigen(Q, only.values = TRUE)$values lambdas <- as.numeric(Re(ev)) + 1i*as.numeric(Im(ev)) # ---------- compute omegas via solving Vandermonde (if roots are real/distinct) ---------- omegas <- rep(NA_real_, m) roots_status <- "not_used" if (all(abs(Im(lambdas)) < 1e-12)) { # real lambdas lam <- as.numeric(Re(lambdas)) if (length(unique(round(lam, 12))) == m) { # distinct real roots: p(n) = sum_k omega_k (1-lam_k) lam_k^{n-1} # Build Vandermonde for n=1..m: p(n) = sum_k A_k lam_k^n # Relation between A_k and omega_k: A_k = omega_k (1 - lam_k) lam_k^{-1} # So for n=1..m: p(n) = sum_k [omega_k (1 - lam_k) lam_k^{n-1}] V <- outer(1:m, lam, function(n, lk) lk^(n-1)) # powers 0..m-1 rhs <- p[1:m] # Solve V * b = rhs where b_k = omega_k (1 - lam_k) b <- as.numeric(solve(V, rhs)) omegas <- b / (1 - lam) # omega_k roots_status <- "real_distinct" } else { roots_status <- "real_but_repeated" } } else { roots_status <- "complex" } # ---------- monotonicity (increasing hazard) check ---------- # We check numerically that r(n+1) >= r(n) for n=1..Ncheck where both finite Ncheck <- min(Nmax - 1, 200L) mono_ok <- TRUE for (n in 1:Ncheck) { if (is.finite(r[n]) && is.finite(r[n+1])) { if (r[n+1] + 1e-12 < r[n]) { mono_ok <- FALSE; break } } } # ---------- limit hazard ---------- # For real-distinct roots, r(n) -> 1 - max(lambdas) # In general (spectral radius rho), asymptotics governed by rho. rho <- max(Mod(lambdas)) r_limit <- if (roots_status == "real_distinct") { 1 - max(Re(lambdas)) } else { 1 - rho # a conservative estimate } # ---------- threshold in (9): (c1/(c1 - c2)) * (1/E[T]) ---------- thresh <- (c1 / (c1 - c2)) * (1 / E_T) # ---------- L(n) and C(n) ---------- # L(n) = [ (pi Q^n u')/(pi Q^n e') ] * [ pi (I - Q^n) (I - Q)^{-1} e' ] + pi Q^n e' - 1 # C(n) = c1 + (c2 - c1) * [ pi Q^n e' ] / [ pi (I - Q^n) (I - Q)^{-1} e' ] L <- numeric(Nmax) C <- numeric(Nmax) invIQ <- solve(I - Q) Qn <- diag(1, m) for (n in 1:Nmax) { Qn <- Qn %*% Q num1 <- as.numeric(pi %*% Qn %*% u) # pi Q^n u' den1 <- as.numeric(pi %*% Qn %*% e) # pi Q^n e' den2 <- as.numeric(pi %*% (I - Qn) %*% invIQ %*% e) # pi (I - Q^n) (I - Q)^{-1} e' # L(n) L[n] <- if (abs(den1) > tol) (num1/den1) * den2 + den1 - 1 else NA_real_ # C(n) C[n] <- if (abs(den2) > tol) c1 + (c2 - c1) * (den1 / den2) else NA_real_ } # ---------- decide N* ---------- # Case A: if mono hazard (increasing) AND r_limit > thresh -> finite unique N* # Find smallest n with L(n) >= c2/(c1 - c2) target <- c2 / (c1 - c2) N_star <- Inf if (mono_ok && (r_limit > thresh)) { idx <- which(L >= target) if (length(idx) > 0) { N_star <- min(idx) } } # Case B: otherwise N* = Inf if (is.infinite(N_star)) { C_star <- c1 / E_T } else { C_star <- C[N_star] } return(list( N_star = N_star, C_star = C_star, details = list( E_T = E_T, thresh = thresh, r_limit = r_limit, monotone_hazard = mono_ok, roots_status = roots_status, lambdas = lambdas, omegas = omegas, p = p, S = S[-1], r = r, L = L, C = C ) )) } 
An Error occurred while handling another error:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response->sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/ErrorHandler.php(136): yii\web\Response->send()
#2 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/ErrorHandler.php(135): yii\web\ErrorHandler->renderException()
#3 [internal function]: yii\base\ErrorHandler->handleException()
#4 {main}
Previous exception:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response->sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(656): yii\web\Response->send()
#2 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/models/Files.php(696): yii\base\Application->end()
#3 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/controllers/FilesController.php(484): faravaghi\filemanager\models\Files->getFile()
#4 [internal function]: faravaghi\filemanager\controllers\FilesController->actionGetFile()
#5 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/InlineAction.php(57): call_user_func_array()
#6 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Controller.php(180): yii\base\InlineAction->runWithParams()
#7 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Module.php(528): yii\base\Controller->runAction()
#8 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Application.php(103): yii\base\Module->runAction()
#9 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(386): yii\web\Application->handleRequest()
#10 /var/www/html/prof-homepages/frontend/web/index.php(18): yii\base\Application->run()
#11 {main}