I am working on an example in which I have obtained parameter estimates using Stan.

In a real life scenario, I would receive more data every week. The data are covariates of a product which are used as predictors to estimate the failure-time of the product. I would also like to put more weighting on the local data since it is a better indicator of the current use-rate/use-environment of the product than the past data.

I would like to update my parameters every week to make inferences about the failure-time for the product. I have read that I have two options:

(1) Approximate the MCMC posterior in a mathematical form and then use this as the prior for the next update.

(2) Run the MCMC code with all of the data $\{y_0, y_1\}$ and use a power prior, where $y_0$ is the past data and $y_1$ is the data obtained one week from now. The likelihood can be written as $L(\theta \mid y_0)^aL(\theta \mid y_1)$, where $a$ is a parameter to be estimated.

Option (1) would be very difficult with the number of correlated parameters in the model.

Option (2) seems like an okay option. However, suppose I now have the following data $\{y_0, y_1,\dots, y_n\}$. I would, once again, like the most recent observations to have more weighting in the likelihood function, but I do not want to give the same weight to all new observations as $n$ gets sufficiently large. Instead, I am imagining the likelihood function would be

\begin{equation} L(\theta \mid y_0)^{a_1}L(\theta \mid y_1)^{a_2}\dots L(\theta \mid y_{n-1})^{a_n}L(\theta \mid y_n). \end{equation}

This way I would be able to assign weights (where the weights are parameters to be estimated) to the observations as they arrive each week. But, this would mean that I have to constantly change my code and add more parameters each week.

I could instead write the likelihood as

\begin{equation} L(\theta \mid y_{old})^{a}L(\theta \mid y_{new}). \end{equation}

Fixing the number of parameters in the model. But how would I decide when a "new" observation becomes old?

Can anyone recommend an approach that they would take?

Is option 1 feasible with approximation methods?

Edit: The likelihood function is given by

\begin{equation} L(\theta_T, \theta_X \mid \text{Data}) = L(\theta_T \mid \text{Failure-time Data, Covariate History}) \times L(\theta_X \mid \text{Covariate History}), \end{equation}

where $\theta_T = (\mu_0, \sigma_0, \beta)$ and $\theta_X = (\eta, \sigma_1, \sigma_2, \rho, \sigma)$.

The expression of the failure-time data conditional on the covariate history is

\begin{equation} \begin{aligned} L&(\theta_T \mid \text{Failure-time Data, Covariate History}) \\ = &\prod_{i=1}^n\{\exp(\beta x_i(t_i)]f_0(u[t_i;\beta, x_i(t_i)], \theta_0)\}^{\delta_i} \times \{1 - F_0(u[t_i;\beta, x_i(t_i)], \theta_0)\}^{1-\delta_i}, \end{aligned} \end{equation}

where $f_0$ and $F_0$ are the pdf and cdf of a Weibull distribution with shape $1/\sigma_0$ and scale $\exp(\mu_0)$.

The likelihood function of the covariate history is

\begin{equation} L(\theta_X \mid \text{Covariate History}) = \prod_{i=1}^n\int_{w_i}\bigg\{\prod_{t_{ij} \leq t_i} f_{\text{NOR}}[x_i(t_{ij}) - \eta - Z_i(t_{ij})w_i; \sigma^2]\bigg\} \times f_{\text{BVN}}(w_i; \Sigma_{w})dw_i, \end{equation}

where, $f_{\text{NOR}}( \cdot ; \sigma^2)$ is the pdf of a univariate normal distribution with mean 0 and variance $\sigma^2$, and $f_{\text{BVN}}( \cdot; \Sigma_{w})$ is the pdf of a bivariate normal distribution with mean 0 and variance-covariance matrix $\Sigma_{w}$.

The following model for $X_i(t_{ij})$ is used

\begin{equation} X_i(t_{ij}) = \eta + Z_i(t_{ij})w_i + \epsilon_{ij}, \end{equation}

where $\eta$ is the mean, $Z_i(t_{ij}) = [1, \log(t_{ij})]$, $w_i = (w_{0i}, w_{1i})' \sim N(0,\Sigma_{w})$, $\epsilon_{ij} \sim N(0, \sigma^2)$, and

\begin{equation*} \Sigma_{w} = \begin{pmatrix} \sigma^2_1& \rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma^2_2 \end{pmatrix}. \end{equation*}

Here, $w_i$ is the random effect for modelling variability in a unit’s covariate process over time.

Since the model for the covariate process is normal, I coded it using the built in normal log probability density in Stan (I am assuming this is correct but the likelihood function for the covariate history is not quite a product of normal densities, so I may be wrong?).

The code I used for the likelihood is given below:

target += weibull_lpdf(u_obs| 1/sigma0, exp(mu0));
target += Beta*x_DFD;
target += weibull_lccdf(u_cens|1/sigma0, exp(mu0));
target += normal_lpdf(y_tt|Mu, sig);

Here $u_{obs}$ are the values of $u$ for the products that failed, $x_{DFD}$ is the value of the covariate at the failure time, $u_{cens}$ are the values of $u$ for the products that are still working, and $y_{tt}$ is the value of the covariate $x_i(t_{ij})$.

I hope this makes sense. Further description can be found in https://www.tandfonline.com/doi/abs/10.1080/00401706.2013.765324

I also note and acknowledge that the model choice was made by the authors. I am using this model for future work.

  • 1
    $\begingroup$ Depending upon the likelihood function, it may be that $\alpha$ can be incorporated into the model cleanly (consider an exponential smoothing model, which effectively downweights past observations in a manner similar to what you are suggesting). Can you be more specific in this regard? $\endgroup$
    – jbowman
    Aug 2, 2019 at 18:29
  • $\begingroup$ @jbowman Thank you for your comment. I have edited by post. $\endgroup$
    – JLee
    Aug 2, 2019 at 22:04

1 Answer 1


From my understanding, option 1 or using a power prior are very similar in the sense that the power prior uses a weighted version of the previously obtained posterior as a new prior. Therefore, the power prior requires that your posterior as a nice mathematical form so that you can approximate it with a known distribution, but this is not always the case unless you are using conjugate priors.

I am not entirely convinced with option 2 since you plan to train the model on all your data and therefore you would lose the apparent benefit of using a power prior. Moreover, if we forget the power prior for a minute, you would expect your model to be valid whether you are using 10 or 20 observations, and by that I mean that you implicitly assume stationarity, i.e. the parameters shouldn't change. In this context, to me, $a$ is more a decision parameter and not something to be inferred.

If you want to use all the data, I would instead suggest that you explicitly model that the parameters are time-varying, then you could use a simple exponential smoothing and infer the corresponding smoothing factor.

For instance, let's consider the parameter $\theta(t)$ that you would define in the model block and it's smoothed version $\theta_s(t)$ that you would define in the transformed parameters block as: $$\theta_s(t) = \alpha \theta(t) + (1 - \alpha) \theta_s(t - 1)$$ Where $\alpha$ is the smoothing factor, to be inferred, which you can relate to the time constant $\tau$ of the learning process by $\alpha = 1 - e^{\frac{-\Delta T}{\tau}}$. You can model the evolution of $\theta(t)$ in the model block with: $$\theta(t + 1) \sim \mathcal{N}(\theta_s(t), \sigma^2)$$

Using this kind of approach, you will not only be able to infer when "a parameter becomes old" (cf. $\alpha$ or $\tau$) but also how much "a parameter is expected to change (cf. $\sigma$).

  • $\begingroup$ I am treating $a$ as an unknown parameter. This means that I need to evaluate a normalisation constant. This requires numeric integration. Do you have any suggestions about how to calculate this? The normalisation constant is described here: dx.doi.org/10.1002/sim.3722 $\endgroup$
    – JLee
    Aug 3, 2019 at 18:45
  • $\begingroup$ Thanks for the link, I will edit my answer in light of this. $\endgroup$
    – Guillem
    Aug 5, 2019 at 8:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.