Estimate the probability that a crack length exceeds a threshold value after N cycles

Question

The length $a$ of a crack after $N$ fatigue cycles $N$ is

$$a(N;C,m) = \left(a_0^{\left(1-\frac m 2\right)}+ C\left(1-\frac m 2\right)B^mN\right)^{\frac{2}{2-m}}$$

where $a_0$ is the initial crack length, $B$ is a value which depends on the range of stresses applied during a fatigue cycle, and $C$ and $m$ are unknown model parameters. We have a dataset $D$ of iid measurements $D=(N_i, a_i)_{i=1}^{n}$ ($n$ is the number of measurements, not to be confused with the number of cycles $N$), where, for fixed $a_0$ and $B$, each time we perform a different number of fatigue cycles $N_i$ and we measure the corresponding final crack length $a_i$.

The goal of the analysis is to estimate the number of cycles $N_c$ such that the probability $Pr\left(a(N;C,m)>a_c\right)$ is equal to a threshold value $p_c$. We work in a Bayesian framework, so we assume Gaussian priors $\ln{C}\sim\mathcal{N}(\mu_C, \sigma_C)$, $m\sim\mathcal{N}(\mu_m, \sigma_m)$ for the model parameters, and measurements errors $\epsilon\sim\mathcal{N}(0, \sigma_{\epsilon})$. In other words, we assume that the likelihood $\mathcal{L}(D|C,m)$ of the data $D$ under our crack growth model is

$$\mathcal{L}(D|C,m)=\frac{1}{(2\pi\sigma_{\epsilon}^2)^{n/2}}e^{-SS/2\sigma_{\epsilon}^2}$$

where the residual sum of squares is

$$SS=\sum_{i=1}^n(a_i-a(N_i;C,m))^2 $$

How can I estimate $N_c$?

A possible solution (if I'm not forgetting something important) could be to:

1. use MCMC to compute the posterior distribution of random parameters $C, m$ $$p(C,m|D)=\frac{\mathcal{L}(D|C,m)p(C)p(m)}{\int\mathcal{L}(D|C,m)p(C)p(m)dCdm}$$
2. sample $C_i, m_i$ from $p(C,m|D)$ (actually, MCMC generates samples from $p(C,m|D)$, so I already have such samples from step 1)
3. substitute the samples $C_i,m_i$ into $a(N;C,m)$, for an array of values $N_j=0,\dots,N_{max}$. For each $N_j$, I thus have a set of samples $a_{ij}=a(N_j;C_i,m_i)$
4. for each value $N_j$ I estimate $Pr\left(a(N;C,m)>a_c\right)$ as $\frac{\#(a_{ij}>a_c)}{\#a_{ij}}$, i.e., the number of samples $a_{ij}>a_c$ divided by the total number of samples $a_{ij}$ (apologies if the notation is awkward/nonstandard: feel free to suggest improvements).
5. I estimate $N_c$ as the smallest $N_j$ such that $\frac{\#(a_{ij}>a_c)}{\#a_{ij}}>p_c$

However, I have two doubts:

• This isn't really much of an estimate. Ideally, given a dataset $D$, a statistical estimation procedure of $N_c$ should return an estimated values as well as some measure of the estimation uncertainty. Here I just got a value for $N_c$. What am I missing?
• for a series of reasons, I'd need to implement this in Python, and to have the whole procedure run really fast (ideally, a few seconds). It looks like MCMC can be slow even for these simple problems, so I was wondering if there could be an alternative way to sample from $p(C,m|D)$. I could compute the normalization factor $Z=\int\mathcal{L}(D|C,m)p(C)p(m)dCdm$ using the tensor product of two 1D Gauss-Hermite quadrature rules, similar to what it's illustrated here. However, how would I then sample from the posterior?

Answer 1

As I understand the formula at the top of the OP gives the expected length which depends on the parameters $C$ and $m$, whereas $a_0$ and $B$ are known quantities. I will denote by $A$ the random crack length and by $\mu_A$ or $\mu_A(N \,\vert\,C,\,m)$ the regression mean. Note that it is essential that $m < 2$. If so, $\mu_A$ is a strictly increasing function of $N$ for given $C$ and $m$, with the reciprocal function given by

\begin{equation} \tag{1} N = \frac{\mu_A^{1 - m/2} - a_0^{1 - m /2}}{1 - m/2} \, \frac{1}{C B ^m} =: \nu(\mu_A \, \vert \,C,\,m). \end{equation}

Maybe the limit for $m \approx 2$ can be used for $m=2$.

From the form of the likelihood, the model is actually a Gaussian non-linear regression with $2$ parameters $C$, $m$ and the error variance $\sigma^2_\epsilon$ seems to be considered as known.

If a crack length $a^\star$ and a probability $p^\star$ are given one can, for any given $C$ and $m$, find the value $N^\star(C,\,m)$ such that $F_{A}(a^\star,\,N^\star \vert\, C,\, m) = 1 - p^\star$ where $F_A(a,\,N\, \vert\, C,\, m) := \text{Pr}\{A \leqslant a \vert\, C,\, m\}$ which depends on $N$. Indeed the corresponding expectation of $A$ is then

$$ \mu_A(N^\star \, \vert \, C,\, m) = a^\star - \sigma_\epsilon \,q_Z(1- p^\star) $$ where $q_Z$ is the quantile function of the standard normal. Hence

\begin{equation} \tag{2} N^\star(C,\,m) = \nu\left\{a^\star - \sigma_\epsilon \,q_Z(1- p^\star) \, \vert \, C,\,m\right\} \end{equation} where the function $\nu$ is given by (1).

As an interesting feature of the Bayesian framework, the so-called predictive (posterior) distribution of the crack length $A$ after $N$ cycles is obtained by "marginalizing out" $C$ and $m$

\begin{equation} \tag{3} F_A(a;\,N) = \iint F_A(a;\, N \, \vert \, C,\, m)\,p(C,\,m \, \vert \, D) \,\text{d}C \, \text{d}m. \end{equation}

While $F_A(a;\, N \, \vert \, C,\, m)$ is a Gaussian distribution with mean $\mu_A(N \, \vert \, C,\, m)$ and variance $\sigma_\epsilon^2$, the distribution $F_A(a;\,N)$ is an over-dispersed version of the Gaussian. It may be heavy-tailed.

Since $N^\star$ is a known function of the parameter $[C,\, m]$ its posterior distribution can be retrieved from that of the parameter. One can give a credible interval on $N^\star$. Yet it makes sense to find instead the value $N^\circ$ such that $F_{A}(a^\star,\,N^\circ) = 1 - p^\star$, hence using the predictive distribution. This value is computed in the OP. It better assesses the risk related to a large crack since it takes into account the uncertainty on the parameters, the predictive distribution being more over-dispersed when the posterior is looser.

To see how the two notions differ, suppose that $K$ nearly independent MCMC iterates $[C^{[k]}, \, m^{[k]}]$ are available.

• One can approximate the posterior distribution of $N^\star(C,\,m)$ by the empirical distribution of the values $N^\star(C^{[k]}, m^{(k]})$ as computed by (2).

• By contrast, to find $N^\circ$ one has to use an approximation of (3) $$ \widetilde{F}_A(a;\,N) := \frac{1}{K} \, \sum_{k} F_{\texttt{norm}}(a; \mu_A^{[k]}, \, \sigma_\epsilon), \qquad \mu_A^{[k]} := \mu_A(N,\, C^{[k]}, \, m^{[k]}). $$ One can use this approximation for $a := a^\star$ and each $N$ in a grid of possible values, then solve approximately $\widetilde{F}_{A}(a^\star,\,N^\circ) = 1 - p^\star$.

Note that in the first computation one does not need to use empirical quantiles since the distribution of $A$ conditional on $C$ and $m$ is known. The second "predictive" approach is the one used in the OP. It does not need a credible interval which is the "first doubt". Yet the estimation uncertainty is accounted for in the result.

The cheapest alternative to MCMC would rely on the computation of the posterior mode, say $[\widehat{C}, \, \widehat{m}]$ and use Laplace's approximation. This should easily be implemented in a language like Python. Gauss-Hermite quadrature rules can indeed do better. Since the parameter is only bi-dimensional, a quadrature can also be used to get the joint posterior density $p(C,\,m \,\vert \,D)$ on a grid of say 1e4 nodes. For instance the bi-dimensional Simpson's rule should work well to get the normalizing constant. Then, for the first computation a simple weighted mean of $N^\star(C,\, m)$ given by (2) can be used to get the posterior mean. For the second (predictive) approach $F_A(a^\star, \,N)$ can be evaluated by using the approximation of (3) obtained by quadrature on a grid of candidate values for $N$. In both cases one does not need to sample from the posterior. Yet in the first computation, getting a credible interval on $N^\star$ would require more work. Of course the log-sum exponential trick must be used. The limitation is that using an extra unknown parameter such as $\sigma_{\epsilon}$ would be a pain.

Anyway, an efficient MCMC sampler such as Stan should be used to investigate the posterior before choosing a suitable approximation. Non-Gaussian priors should be used for $C$ and $m$ to cope with the constraints $C >0$ and $m<2$.

EDIT The condition $m< 2$ is actually not required, because $\mu_A(N)$ turns out to be an increasing function of $N$ even when $m > 2$ with the reciprocal function (1) above. However the quantity for which a non-integer exponent $2/(2-m)$ is taken should be positive so $$ a_0 + C \,\frac{2 -m}{2}\,\left[ \sqrt{a_0} B \right]^m N > 0 $$ must hold for the range of values of interest for $N$. This may impact the choice of the prior for $C$ and $m$.