Parameter Updates in Stan 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. 
 A: 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$).
