--- title: "Bivariate Chen distributions with BvCh" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Bivariate Chen distributions with BvCh} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5) set.seed(2024) ``` ## Introduction `BvCh` provides the probability density, survival and hazard functions, random generation, graphical displays and maximum likelihood estimation for three bivariate Chen type lifetime models: | Model | Construction | Singular part? | |:------|:-------------|:---------------| | **BvCh** | Marshall-Olkin type, `min` of three Chen variables | yes, on the diagonal | | **BEC** | Marshall-Olkin type with extended Chen variables, adding `lambda` | yes, on the diagonal | | **IFGM-Chen** | Inverse FGM copula with Chen marginals | no | All univariate Chen building blocks are taken from the [reliaR](https://cran.r-project.org/package=reliaR) package. ```{r load} library(BvCh) ``` ## The BvCh model With $U_i \sim \mathrm{Chen}(\alpha_i,\beta)$ independent and $Z_1=\min(U_1,U_3)$, $Z_2=\min(U_2,U_3)$, the joint survival function is $$S(z_1,z_2)=e^{\alpha_1(1-e^{z_1^\beta})}e^{\alpha_2(1-e^{z_2^\beta})}e^{\alpha_3(1-e^{w^\beta})}, \qquad w=\max(z_1,z_2),$$ which is a product of three univariate Chen survival functions. The density is piecewise, with a branch on each side of the diagonal and a singular branch on the diagonal itself. ```{r bvch} a1 <- 1.5; a2 <- 2.5; a3 <- 0.5; b <- 1.8 z1 <- c(0.4, 0.9, 0.7, 1.3) z2 <- c(0.6, 0.4, 0.7, 1.1) cbind(z1, z2, pdf = dbvch(z1, z2, a1, a2, a3, b), surv = sbvch(z1, z2, a1, a2, a3, b), haz = hbvch(z1, z2, a1, a2, a3, b)) ``` The survival function can be checked against the product form, and the hazard is exactly $f/S$. ```{r check} S <- sbvch(z1, z2, a1, a2, a3, b) Sref <- schen(z1, b, a1) * schen(z2, b, a2) * schen(pmax(z1, z2), b, a3) c(max_abs_diff_surv = max(abs(S - Sref)), max_abs_diff_haz = max(abs(hbvch(z1, z2, a1, a2, a3, b) - dbvch(z1, z2, a1, a2, a3, b) / S))) ``` ### The mixture representation Proposition 1 writes the density as a mixture of an absolutely continuous part and a singular part on the diagonal. Both components are available directly and the mixture identity holds to machine precision. ```{r mixture} K <- a1 + a2 + a3 ac <- dbvch_ac(z1, z2, a1, a2, a3, b) sg <- dbvch_sing(z1, z2, a1, a2, a3, b) mix <- (a1 + a2) / K * ac + a3 / K * sg data.frame(ac = ac, sing = sg, mixture = mix, dbvch = dbvch(z1, z2, a1, a2, a3, b)) ``` ### Mass carried by each region Numerically integrating the density over the three regions recovers $\alpha_i/K$ exactly, confirming that the density is properly normalised. ```{r mass} m1 <- integrate(function(z) vapply(z, function(v) integrate(function(v2) dbvch(v, v2, a1, a2, a3, b), v, Inf, rel.tol = 1e-10, subdivisions = 2000L)$value, numeric(1)), 0, Inf, rel.tol = 1e-8, subdivisions = 2000L)$value m2 <- integrate(function(z) vapply(z, function(v) if (v > 0) integrate(function(v2) dbvch(v, v2, a1, a2, a3, b), 0, v, rel.tol = 1e-10, subdivisions = 2000L)$value else 0, numeric(1)), 0, Inf, rel.tol = 1e-8, subdivisions = 2000L)$value u <- seq(1e-8, 6, length.out = 200001L) m0 <- a3 / K * sum(reliaR::dchen(u, b, K)) * (u[2] - u[1]) rbind(f1 = m1, f2 = m2, diagonal = m0, total = m1 + m2 + m0, expected = c(a1/K, a2/K, a3/K, 1)) ``` ### Simulation Random generation draws three independent Chen variables and applies the Marshall-Olkin construction. ```{r sim} sim <- rbvch(100000, a1, a2, a3, b) q <- seq(0.4, 1.6, length.out = 6) data.frame(q = q, theo = 1 - schen(q, b, a1) * schen(q, b, a3), emp = sapply(q, function(v) mean(sim[, 1] <= v))) ``` The frequency of ties also matches $\alpha_3/K$: ## The BEC model The Bivariate Extended Chen model replaces the Chen building block with the extended Chen distribution, whose survival function is $[1+\lambda(e^{x^\beta}-1)]^{-\alpha}$. ```{r bec} x1 <- c(0.4, 0.9, 0.7, 1.3); x2 <- c(0.6, 0.4, 0.7, 1.1) cbind(x1, x2, pdf = dbec(x1, x2, a1, a2, a3, b, lambda = 2), surv = sbec(x1, x2, a1, a2, a3, b, lambda = 2), haz = hbec(x1, x2, a1, a2, a3, b, lambda = 2)) ``` Off the diagonal the joint hazard factorises into the componentwise hazard gradients: ```{r becgrad} off <- x1 != x2 h12 <- hbec1(x1, x2, a1, a2, a3, b, 2) * hbec2(x1, x2, a1, a2, a3, b, 2) data.frame(hbec = hbec(x1, x2, a1, a2, a3, b, 2), hbec1_times_hbec2 = h12, off_diagonal = off) ``` ## The IFGM-Chen model This model pairs the Inverse FGM copula with Chen marginals and is absolutely continuous. ```{r ifgm} t1 <- 1.2; b1 <- 1.3; t2 <- 1.5; b2 <- 1.1; g <- 0.6; om <- 0.4 cbind(x = 0.5, y = 0.9, reliability = relifgmchen(0.5, 0.9, t1, b1, t2, b2, g, om), pdf = difgmchen(0.5, 0.9, t1, b1, t2, b2, g, om), haz = hifgmchen(0.5, 0.9, t1, b1, t2, b2, g, om)) ``` ```{r copula} data.frame(C = c(cifgm(0.3, 0.7, g, om), cifgm(0.5, 0.5, g, om)), tau = tau_ifgm(g, om), rho = rho_ifgm(g, om)) ``` ```{r tauscan} data.frame(gamma = c(0, 0.3, 0.6, 0.9, 1), tau = sapply(c(0, 0.3, 0.6, 0.9, 1), tau_ifgm, omega = 0.4), rho = sapply(c(0, 0.3, 0.6, 0.9, 1), rho_ifgm, omega = 0.4)) ``` > With the copula formula exactly as given in the source document, $\tau$ and > $\rho$ are non-positive for $\gamma,\omega \geq 0$. The source quotes > $\tau \approx 0.156>0$ for the same parameters, which is inconsistent with > the copula formula it also quotes. See `?cifgm` for details. ## Estimation `fitbvch`, `fitbec` and `fitifgmchen` estimate the models by maximum likelihood. For the Marshall-Olkin type models the observed-data log-likelihood is exact, because the contribution of an observation is fully determined by which side of the diagonal it falls on, so direct numerical maximisation reaches the same point as the EM algorithm. ```{r fit} d <- rbvch(1000, a1, a2, a3, b) fit <- fitbvch(d) fit ``` Standard errors come from the numerical Hessian and are reported on the natural scale. ```{r se} coef(fit) fit$se ``` ### Goodness of fit Because $Z_1=\min(U_1,U_3)$, the marginal of $Z_1$ is $\mathrm{Chen}(\alpha_1+\alpha_3,\beta)$; likewise for $Z_2$, and $\min(Z_1,Z_2)=\min(U_1,U_2,U_3)$ is $\mathrm{Chen}(\alpha_1+\alpha_2+\alpha_3,\beta)$. ```{r ks} ksfit_bvch(fit, d) ``` ### Model comparison ```{r compare} f1 <- fitbvch(rbvch(500, a1, a2, a3, b)) f2 <- fitbec(rbec(500, a1, a2, a3, b, 2)) f3 <- fitifgmchen(rifgmchen(500, t1, b1, t2, b2, g, om)) compare_models(f1, f2, f3) ``` > The IFGM-Chen copula parameters `gamma` and `omega` are only weakly > identified: the marginal parameters are estimated precisely but the > likelihood is nearly flat in the copula direction. See `?fitifgmchen`. ## Graphics The absolutely continuous part of the BvCh density is plotted for a range of $\beta$ values; as $\beta$ increases the surface becomes more peaked near the origin. ```{r surface, fig.height = 6} plot_surface_bvch(alpha1 = a1, alpha2 = a2, alpha3 = a3, beta.vec = c(0.8, 1.8, 5, 10)) ``` Contour plot of the absolutely continuous part, with the singular support marked. ```{r contour} plot_contour_bvch(alpha1 = a1, alpha2 = a2, alpha3 = a3, beta = b) ``` Reliability surface of the IFGM-Chen model. ```{r relsurface} plot_reliability_ifgmchen(theta1 = t1, beta1 = b1, theta2 = t2, beta2 = b2, gamma = g, omega = om) ``` ## References * Chen, Z. (2000). A new two-parameter lifetime distribution with bathtub shape or increasing failure rate function. *Statistics & Probability Letters*, 49, 155-161. * Nelsen, R. B. *An Introduction to Copulas*, 2nd edition, Springer, 2006. * Klein, M. and Schumacher, C. (2002). Bivariate survival and three dimensional Marshall-Olkin type copulas. *Journal of the Royal Statistical Society B*, 64, 721-734. ```{r ties} c(empirical = mean(sim[, 1] == sim[, 2]), expected = a3 / K) ```