Parameter estimation in a drifting normal distribution Suppose we observe $n$ pairs of points $(a_1,b_1),~(a_2,b_3)~..~(a_n,b_n)$. The underlying data generating process is known to be as follows:

*

*$u_i \sim N(0,1)$ and $a_i \sim N(u_i,1)$

*$b_i \sim N(u_i + k*a_i, 1)$
independently for each $i$ and $k$ is a constant. What is a good unbiased estimator for $k$?

I am going to start to discuss @bleh's idea in order to explain why my solution seems more robust.

*

*@bleh's idea:

Assume:

*

*$N_{i,j} \sim N(\mu_j,\sigma^2_j)$ are independent random variables $\forall i,j$

*$U_i=N_{i,1}$

*$A_i\sim N(U_i,\sigma_2^2)\Rightarrow A_i=N_{i,1}+N_{i,2}$

*$B_i \sim N(U_i+kA_i,\sigma_3^2)\Rightarrow B_i=(1+k)N_{i,1}+kN_{i,2}+N_{i,3}$
As a consequence, $B_i\sim N((1+k)\mu_1+k\mu_2+\mu_3,k^2(\sigma^2_1+\sigma^2_2)+2k\sigma^2_1+\sigma^3_1)$
It is tempting to compute an estimator $\widehat{\sigma^2_B}$ of the variance of $B$ based on the squared differences from the means $\sum_{i=1}^n\left(b_i-\mathbb{E}(B_i)\right)^2$, to equate it to $\mathbb{V}ar(B)$, and to solve for $k$.
This is a problem for 2 main reasons:

*

*To solve this polynomial, the value of the discriminant $\Delta$ will depend on the parameters $\sigma^2_j$ as well as on the estimator $\widehat{\sigma^2_B}$. However, in some unlikely cases (ex: very small $\widehat{\sigma^2_B}$ being computed despite $\sigma^2_3$ being large), the calculation will fail (due to $\Delta<0$) because the observed values $b_i$ were not dispersed enough compared to the expected variance $\sigma^2_3$.

*Even if this process worked, determining an estimator $\widehat{\sigma^2_B}$ that would allow us to find an unbiased estimator for $k$ after all this steps, including the calculations of squared deviations to a mean, and a square root, would be very difficult.

In addition, we would end up loosing some information because we would use our observations $b_i$, but not $a_i$.
Hence, I propose my solution:


*Use bayesian statistics and the normal conjugate prior.

You can see here that:
If $\mu \sim N(\mu_0,\sigma^2_0)$ and $X_i\sim N(\mu,\sigma^2)$ (the expected value of $X_i$ depends on a random value $\mu$), then:

*

*$\mu\vert x_i, i=1,...,n \sim N(\mu_0\prime, \sigma^{2}_0\prime)$, with $\mu_0\prime = \frac{1}{\frac{1}{\sigma^2_0}+\frac{n}{\sigma^2}}\left(\frac{\mu_0}{\sigma^2_0}+\frac{\sum_{i=1}^{n}x_i}{\sigma^2}\right)$ and $\sigma_0^2\prime=\frac{1}{\frac{1}{\sigma_0^2}+\frac{n}{\sigma^2}}$ (more observations of $x_i$ allow to estimate a more accurate distribution for $\mu$)

*$X_{i+1} \sim N(\mu_0\prime, \sigma^{2}_0\prime + \sigma^2)$ (updated knowledge about $\mu$ allow to estimate a more accurate distribution for future observations of $X$)

Apply this with your example:

*

*$U_i\sim N(\mu_U,\sigma_U^2)$

*$A_i \sim N(U_i, \sigma_A^2)$

*$B_i \sim N(U_i+kA_i, \sigma_B^2)$
So this is a "tower" of normal variables using the previous variables as conjugate prior.
The observations we have are the couples $\{a_i,b_i\}$.
Let's start by the end:
$B_i\vert A_i \sim N(U_i+kA_i, \sigma_B^2)\vert A_i$
In this "equation", we know $A_i$, but not $U_i$. But based on our bayesian knowledge, we can determine the distribution of $U_i \vert A_i$:
Based on the equations above:
$$U_i \lvert A_i \sim N\left(\tilde{\mu_i}=\frac{1}{\frac{1}{\sigma_U^2}+\frac{1}{\sigma_A^2}}\left(\frac{\mu_U}{\sigma_U^2}+\frac{A_i}{\sigma_A^2}\right),\tilde{\sigma^2_{U,i}}=\frac{1}{\frac{1}{\sigma_U^2}+\frac{1}{\sigma^2_A}}\right)$$
Hence, you will have:
$$B_i\vert A_i \sim N(kA_i+\tilde{\mu_i},\sigma^2_B+\tilde{\sigma^2_{U,i}})$$
Hence, $\mathbb{E}(B_i\vert A_i)=kA_i+\tilde{\mu_i}$, and we can propose a range of estimators $\hat{k}_i=\frac{b_i-\tilde{\mu_i}}{a_i}$. They are unbiased, and each has a variance which can be computed with the parameters above.
Since you have $i$ estimators $\hat{k}_i$ with variance $\sigma^2_{k, i}=\frac{\sigma^2_B+\tilde{\sigma^2_{U,i}}}{a_i^2}$, you can determine a linear combination of these estimators minimising the variance of the final estimator:
$\min_w w' \Sigma w$, st $w'\mathbb{1}=1$ where $\Sigma$ is the covariance matrix of the estimators $\hat{k}_i$.
Finally, the unbiased estimator with the smallest variance will be:
$$w^*=\frac{\Sigma^{-1}\mathbb{1}}{\mathbb{1}'\Sigma^{-1}\mathbb{1}}$$
