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$?

  • $\begingroup$ We have $u = n_0$, $a = n_0 + n_1$, $b = n_0 + k*(n_0 + n_1)+ n_2 = (1+k)*n_0 + k*n_1 + n_2$ where $n_i$'s are all standard normals. Variance of $b$ then is $(1+k)^2+k^2+1$. Equating this with empirical estimate of variance and solving gives you a pretty good estimate for $k$. But it is not unbiased. $\endgroup$
    – user284807
    Jul 25, 2022 at 23:13
  • $\begingroup$ @bleh So simple, but so efficient... I developped formulas based on quadratic optimisations of bayesian estimators... $\endgroup$
    – FP0
    Jul 25, 2022 at 23:19
  • 1
    $\begingroup$ Haha. It would be interesting to see how to estimate the drift in general. Say if we have $b \sim N(u+f(a), 1)$ for smooth function $f$ and we need to estimate it. Maybe your methods will come handy then? $\endgroup$
    – user284807
    Jul 25, 2022 at 23:25
  • $\begingroup$ I tried to write a very comprehensive answer. @bleh, user363943, may I kindly ask you to help me to check that what I wrote is correct and contains no typo? $\endgroup$
    – FP0
    Jul 26, 2022 at 1:31

1 Answer 1


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

  1. @bleh's idea:


  • $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:

  1. 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$.
  2. 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:

  1. 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:



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.