From 63975ff1f6306cdd9b6b723956731fec5a82ef15 Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Fri, 20 Sep 2024 14:37:13 +0100 Subject: [PATCH 1/6] Add Weibull + minor typo fixes --- vignettes/why-it-works.Rmd | 64 ++++++++++++++++++++++++++++++++++---- 1 file changed, 58 insertions(+), 6 deletions(-) diff --git a/vignettes/why-it-works.Rmd b/vignettes/why-it-works.Rmd index a02ff62..0c403c0 100644 --- a/vignettes/why-it-works.Rmd +++ b/vignettes/why-it-works.Rmd @@ -154,6 +154,7 @@ Applying a uniform primary event time distribution to equation \@ref(eq:survival $$ Q_{S_+}(t) = Q_T(t + w_P) + { 1 \over w_P} \int_0^{w_P} f_T(t+p) p~ dp. $$ + This is analytically solvable whenever the upper distribution function of $T$ is known and the mean of $T$ is analytically solvable from its integral definition. In each case considered below it is easier to change the integration variable: @@ -177,12 +178,12 @@ can be reduced to an analytic expression. **General Discrete censored delay distribution** -First we note that equation \@ref(eq:disccensprob) can be written using the difference operator: $f_n = -\Delta_1 Q_{S_+}(n-1)$. We can insert this expression into equation \@ref(eq:unifprim) to give the discrete censored delay distribution for uniformly distributed primary event times: +First, we note that equation \@ref(eq:disccensprob) can be written using the difference operator: $f_n = -\Delta_1 Q_{S_+}(n-1)$. We can insert this expression into equation \@ref(eq:unifprim) to give the discrete censored delay distribution for uniformly distributed primary event times: $$ \begin{aligned} -f_n &= \Delta_1\Big[(n-1) \Delta_1F_T(n-1)\Big] - \Delta_1Q_T(n) + \Delta_1\Big[ \int_{n-1}^n f_T(z) z ~dz \Big] \\ -&= (n+1)F_T(n+1) + (n-1)F_T(n-1) - 2nF_T(n) + \Delta_1\Big[ \int_{n-1}^n f_T(z) z ~dz \Big]. +f_n &= \Delta_1\Big[(n-1) \Delta_1F_T(n-1)\Big] - \Delta_1Q_T(n) - \Delta_1\Big[ \int_{n-1}^n f_T(z) z ~dz \Big] \\ +&= (n+1)F_T(n+1) + (n-1)F_T(n-1) - 2nF_T(n) - \Delta_1\Big[ \int_{n-1}^n f_T(z) z ~dz \Big]. \end{aligned} (\#eq:disccensunifprim) $$ @@ -219,7 +220,7 @@ $$ **Survival function of $S_{+}$ for Gamma distribution** -By substituting equation \@ref(eq:gammapartexp) into equation \@ref(eq:disccensunifprim) we can solve for both the survival function of $S_+$ in terms of analytically available functions: +By substituting equation \@ref(eq:gammapartexp) into equation \@ref(eq:disccensunifprim) we can solve for the survival function of $S_+$ in terms of analytically available functions: $$ \begin{aligned} @@ -270,7 +271,7 @@ $$ **Survival function of $S_{+}$ for Log-Normal distribution** -By substituting equation \@ref(eq:lognormpartexp) into equation \@ref(eq:disccensunifprim) we can solve for both the survival function of $S_+$ in terms of analytically available functions: +By substituting equation \@ref(eq:lognormpartexp) into equation \@ref(eq:disccensunifprim) we can solve for the survival function of $S_+$ in terms of analytically available functions: $$ \begin{aligned} @@ -281,7 +282,7 @@ $$ **Log-Normal discrete censored delay distribution** -By substituting \@ref(eq:lognormpartex) into \@ref(eq:disccensprob) we get the discrete censored delay distribution in terms of analytically available functions: +By substituting \@ref(eq:lognormpartexp) into \@ref(eq:disccensprob) we get the discrete censored delay distribution in terms of analytically available functions: $$ \begin{aligned} @@ -292,5 +293,56 @@ f_n &= (n+1) F_T(n+1; \mu, \sigma) + (n-1) F_T(n-1; \mu, \sigma) - 2 n F_T(n; \m \end{aligned} $$ +#### Weibull distribution + +The Weibull distribution has the density function: + +$$ +f_T(z;\lambda,k) = +\begin{cases} +\frac{k}{\lambda}\left(\frac{z}{\lambda}\right)^{k-1}e^{-(z/\lambda)^{k}}, & z\geq0 ,\\ +0, & z<0, +\end{cases} +$$ +And distribution function: + +$$ +F_T(z;\lambda,k))=\begin{cases}1 - e^{-(z/\lambda)^k}, & z\geq0,\\ 0, & z<0.\end{cases} +$$ +Where $\Phi$ is the standard normal distribution function. + +**Weibull partial expectation** + +We use an integration substitution $y = (z / \lambda)^k$. This gives: + +$$ +\begin{aligned} +\int_t^{t+w_P} z~ f_T(z; \lambda,k) dz &= \lambda\int_{(t / \lambda)^k}^{((t + w_P) / \lambda)^k} y^{1/k} e^{-y} dy\\ +&= \lambda \Delta_{w_P} g(t; \lambda,k) +\end{aligned} (\#eq:weibullpartexp) +$$ + +Where $g(t; \lambda, k) = \gamma(1 + 1/k, (t/\lambda)^k)$ is a reparametrisation of the lower incomplete gamma function. + +**Survival function of $S_{+}$ for Weibull distribution** + +By substituting equation \@ref(eq:weibullpartexp) into equation \@ref(eq:disccensunifprim) we can solve for the survival function of $S_+$ in terms of analytically available functions: + +$$ +Q_{S+}(t ;\lambda,k) = Q_T(t + w_P; \lambda,k) + { 1 \over w_P} \Big[ \lambda \Delta_{w_P} g(t; \lambda,k) - t\Delta_{w_P}F_T(t; \lambda,k)\Big]. +$$ + +**Weibull discrete censored delay distribution** + +By substituting \@ref(eq:weibullpartexp) into \@ref(eq:disccensprob) we get the discrete censored delay distribution in terms of analytically available functions: + +$$ +\begin{aligned} +f_n &= (n+1)F_T(n+1) + (n-1)F_T(n-1) - 2nF_T(n) - \Delta_1\Big[ \int_{n-1}^n f_T(z) z ~dz \Big] \\ +&= (n+1)F_T(n+1) + (n-1)F_T(n-1) - 2nF_T(n) - \lambda \Delta_1^{(2)} g(n-1; \lambda,k) \\ +&= (n+1)F_T(n+1) + (n-1)F_T(n-1) - 2nF_T(n) \\ +&+ \lambda [2 g(n; \lambda,k) - g(n+1; \lambda,k) - g(n-1; \lambda,k)] \qquad n = 0, 1, \dots. +\end{aligned} +$$ ## References From 74e715e5c46928604ea5315f92ca8ca08af147ad Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Fri, 20 Sep 2024 14:51:24 +0100 Subject: [PATCH 2/6] small rewrite and acknowledge cori et al --- vignettes/why-it-works.Rmd | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/vignettes/why-it-works.Rmd b/vignettes/why-it-works.Rmd index 0c403c0..a4bcc5c 100644 --- a/vignettes/why-it-works.Rmd +++ b/vignettes/why-it-works.Rmd @@ -322,7 +322,12 @@ $$ \end{aligned} (\#eq:weibullpartexp) $$ -Where $g(t; \lambda, k) = \gamma(1 + 1/k, (t/\lambda)^k)$ is a reparametrisation of the lower incomplete gamma function. +Where + +$$ +g(t; \lambda, k) = \gamma(1 + 1/k, (t/\lambda)^k) = \frac{1}{k}\gamma(1/k, (t/\lambda)^k) - \frac{t}{\lambda}e^{-(t/\lambda)^k}. +$$ +is a reparametrisation of the lower incomplete gamma function. **Survival function of $S_{+}$ for Weibull distribution** @@ -345,4 +350,6 @@ f_n &= (n+1)F_T(n+1) + (n-1)F_T(n-1) - 2nF_T(n) - \Delta_1\Big[ \int_{n-1}^n f_ \end{aligned} $$ +Which was also found by Cori _et al_ [@cori2013new]. + ## References From 533e40552a5b1db196ea9fb781dede198d5d88c9 Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Fri, 20 Sep 2024 14:57:09 +0100 Subject: [PATCH 3/6] Update WORDLIST --- inst/WORDLIST | 2 ++ 1 file changed, 2 insertions(+) diff --git a/inst/WORDLIST b/inst/WORDLIST index a1e4485..9d050fe 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -60,6 +60,7 @@ primaryeventdistributions propto pwindow qquad +reparametrisation rfloor rightarrow rprimary @@ -73,4 +74,5 @@ survivalfunc swindow u unifprim +weibullpartexp zsusswein From 310c332a787d32034bd8a93f5c0b7bd0f1038041 Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Fri, 20 Sep 2024 14:57:36 +0100 Subject: [PATCH 4/6] Update NEWS.md --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 626bbc2..c5ac224 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,7 @@ This is the current development version. * Add `{touchstone}` based benchmarks for benchmarking R utility functions, and fitting the `stan` and `fitdistplus` models. * Added a "How it works" vignette. * Added R infrastructure for analytical solutions via the `primary_censored_dist` S3 class. +* Added Weibull analytical solution to "How it works" vignette. # primarycensoreddist 0.4.0 From e28bb0148b2eb74496249ab651d0ac5d37ac08ef Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Fri, 20 Sep 2024 16:36:26 +0100 Subject: [PATCH 5/6] more steps and caught some typos --- vignettes/why-it-works.Rmd | 32 ++++++++++++++++++++++++++------ 1 file changed, 26 insertions(+), 6 deletions(-) diff --git a/vignettes/why-it-works.Rmd b/vignettes/why-it-works.Rmd index a4bcc5c..d17ae69 100644 --- a/vignettes/why-it-works.Rmd +++ b/vignettes/why-it-works.Rmd @@ -174,7 +174,9 @@ $$ \int_t^{t+w_P} f_T(z) z~ dz (\#eq:partexp) $$ -can be reduced to an analytic expression. +can be reduced to an analytic expression. + +The insight here is that this will be possible for any distribution where the average of the distribution can be calculated analytically, which includes commonly used non-negative distributions such as the Gamma, Log-Normal and Weibull distributions. **General Discrete censored delay distribution** @@ -210,6 +212,8 @@ Where $\gamma$ is the lower incomplete gamma function. **Gamma partial expectation** +We know that the full expectation of the Gamma distribution is $\mathbb{E}[T] = k\theta$, which can be calculated as a standard integral. Doing the same integral for the partial expectation gives: + $$ \begin{aligned} \int_t^{t+w_P} f_T(z) z~ dz &= {1 \over \Gamma(k) \theta^k} \int_t^{t+w_P}z z^{k-1} \exp(-z/\theta)~dz \\ @@ -258,12 +262,19 @@ Where $\Phi$ is the standard normal distribution function. **Log-Normal partial expectation** -We use an integration substitution $y = (\ln z - \mu) / \sigma$ and using the shift property of the standard normal distribution function. This gives: +We know that the full expectation of the Log-Normal distribution is $\mathbb{E}[T] = e^{\mu + \frac{1}{2} \sigma^2}$, which can be calculated by integration with the integration substitution $y = (\ln z - \mu) / \sigma$. This has transformation Jacobian: + +$$ +\frac{dz}{dy} = \sigma e^{\sigma y + \mu}. +$$ + +Doing the same integral for the partial expectation, and using the same integration substitution gives: $$ \begin{aligned} -\int_t^{t+w_P} z~ f_T(z; \mu, \sigma) dz &= {1 \over \sigma \sqrt{2\pi}} \int_{(\ln t - \mu)/\sigma}^{(\ln(t+w_P) - \mu)/\sigma} e^{\sigma y + \mu} e^{-y^2/2} dy\\ -&= e^{\mu + \frac{1}{2} \sigma^2} \int_{(\ln t - \mu)/\sigma}^{(\ln(t+w_P) - \mu)/\sigma} e^{-(y- \sigma)^2/2} dy \\ +\int_t^{t+w_P} z~ f_T(z; \mu, \sigma) dz &= {1 \over \sigma \sqrt{2\pi}} \int_t^{t+w_P} \exp\left( - {(\log(z) - \mu)^2 \over 2 \sigma^2} \right) dz \\ +&= {1 \over \sqrt{2\pi}} \int_{(\ln t - \mu)/\sigma}^{(\ln(t+w_P) - \mu)/\sigma} e^{\sigma y + \mu} e^{-y^2/2} dy\\ +&= {e^{\mu + \frac{1}{2} \sigma^2} \over \sqrt{2 \pi} } \int_{(\ln t - \mu)/\sigma}^{(\ln(t+w_P) - \mu)/\sigma} e^{-(y- \sigma)^2/2} dy \\ &= e^{\mu + \frac{1}{2} \sigma^2} \Big[\Phi\Big({\ln(t+w_P) - \mu \over \sigma} - \sigma\Big) - \Phi\Big({\ln(t) - \mu \over \sigma} - \sigma\Big) \Big]\\ &= e^{\mu + \frac{1}{2} \sigma^2} \Delta_{w_P}F_T(t; \mu + \sigma^2, \sigma). \end{aligned} (\#eq:lognormpartexp) @@ -313,11 +324,20 @@ Where $\Phi$ is the standard normal distribution function. **Weibull partial expectation** -We use an integration substitution $y = (z / \lambda)^k$. This gives: +We know that the full expectation of the Weibull distribution is $\mathbb{E}[T] = \lambda \Gamma(1 + 1/k)$, which can be calculated by integration using the integration substitution $y = (z / \lambda)^k$, which has transformation Jacobian: + +$$ +\frac{dz}{dy} = \frac{\lambda}{k}y^{1/k - 1}. +$$ + +Doing the same integral for the partial expectation, and using the same integration substitution gives: $$ \begin{aligned} -\int_t^{t+w_P} z~ f_T(z; \lambda,k) dz &= \lambda\int_{(t / \lambda)^k}^{((t + w_P) / \lambda)^k} y^{1/k} e^{-y} dy\\ +\int_t^{t+w_P} z~ f_T(z; \lambda,k) dz &= \int_t^{t+w_P} \frac{kz}{\lambda}\left(\frac{z}{\lambda}\right)^{k-1}e^{-(z/\lambda)^{k}} dz \\ +&= k\int_t^{t+w_P} \left(\frac{z}{\lambda}\right)^{k}e^{-(z/\lambda)^{k}} dz \\ +&= \lambda k \int_{(t / \lambda)^k}^{((t + w_P) / \lambda)^k} y y^{1/k - 1} e^{-y} dy \\ +&= \lambda\int_{(t / \lambda)^k}^{((t + w_P) / \lambda)^k} y^{1/k} e^{-y} dy\\ &= \lambda \Delta_{w_P} g(t; \lambda,k) \end{aligned} (\#eq:weibullpartexp) $$ From a4d5a210a43efaa13a7b68ff0728d8cdaccb5558 Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Fri, 20 Sep 2024 16:43:55 +0100 Subject: [PATCH 6/6] Update WORDLIST --- inst/WORDLIST | 1 + 1 file changed, 1 insertion(+) diff --git a/inst/WORDLIST b/inst/WORDLIST index 9d050fe..9068137 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -40,6 +40,7 @@ gammapartexp geq hexsticker int +kz leq lfloor lim