# Bayesian Update for two Normal Random Variables following one observation of difference between them

Let's say we have a Normal Random Variable with unknown mean and variance. We want to estimate the mean and variance given a normal choice for a prior and some data. Let's call it X.

Also there's another Normal Random Variable, also with unknown mean and variance It has it's own prior, which can be similarly parameterized. Let's call it Y.

Then we observe the following: d = x - y where x and y are occurrences of X and Y respectively. We don't actually know what x and y are; all we see is d. We wish to update the hyperparameters of X and Y based on observing this data.

My questions are:

1. Is this a known problem with a solution worked out?
2. If not, does anyone have any ideas of how to attack it?

If $X \sim \mathcal{N}(\mu_X, \sigma_X^2)$ and $Y \sim \mathcal{N}(\mu_Y, \sigma_Y^2)$, and $X$ and $Y$ are independent, then

$$X-Y \sim \mathcal{N}(\mu_X-\mu_Y, \sigma_X^2+\sigma_Y^2)$$

This means that your question simplifies to asking how many real numbers $a,b$ such that $a-b = c$ exist, and how many real numbers $d,e$ such that $d+e=f$ exist? In both cases you have unlimited number of possibilities. So the model you describe is not identifiable (there is no unique solution for it) and it cannot be estimated.

Notice: As Glen_b noted in the comments (see comments under this answer), theoretically it may be possible to estimate such model when using informative priors, however I do not see how it could be implemented.

### Comment to code posted by adamwlev in answer to his question

Answers are not meant to comment other answers, but there is too much to be said to use comments for it, so I'm going to refer to the code in my answer. Moreover, while OP considers his solution to answer the question, it seems that commenting it also answers his misconceptions about mentioned model.

First, notice that your code doesn't work. To convince yourself, try using realistic example for testing the code that in fact is a difference of two normals and check if your code correctly finds the appropriate parameters. For example, try the following data and check if the parameters are recovered

ds <- c(rnorm(200,7,10) - rnorm(200,2,10)) ##simulated data


...they will not. But, as Glen_b noticed, Bayesian estimation can work in here by employing prior, out-of-data knowledge into the model to facilitate identifiability. Let's make one more change in your code and assume perfect knowledge about the searched parameters, where your priors are equal to the true parameters, so the model should not have any problems with finding them (as they are known in advance). So, let's change the following lines:

p_x <- dnorm(x,7,10)/sum(dnorm(x,7,10)) ##prior probabilities for X
p_y <- dnorm(y,2,10)/sum(dnorm(y,2,10)) ##prior probabilities for Y


It will diverge from the initial correct values assumed in priors. In fact, it doesn't work exactly for the reason outlined in my answer: such model is not identifable.

The presented algorithm does not estimate the model described in the question. The model that was described was

$$X_i \sim \mathrm{Normal}(\mu_X, \sigma_X^2) \\ Y_i \sim \mathrm{Normal}(\mu_Y, \sigma_Y^2) \\ D_i = X_i - Y_i \sim \mathrm{Normal}(\mu_X-\mu_Y, \sigma_X^2 + \sigma_Y^2)$$

while the model described in the algorithm is something closer to

$$X_i \sim \mathrm{Normal}(\mu_X, \sigma_X^2) \\ Y_i \sim \mathrm{Normal}(\mu_Y, \sigma_Y^2) \\ D_i \sim \mathrm{Normal}(X_i - Y_i, \sigma_D^2)$$

So $D$ in here is not a random variable that is a difference of two other random variables $X,Y$, but instead you have $n$ random variables $D_i$ where each of them has different mean equal to $X_i-Y_i$, where you have only one observation per each random variable $D_i$. So in fact from single datapoint you want to estimate two parameters, $X_i,Y_i$ -- you have not identifable model and insufficient data for this. Moreover, calling exactly the same likelihood twice

p_x1 <- rowSums(outer(p_x,p_y)*dnorm(d,outer(x,y,"-"),liklihood_shape))
p_x1 <- p_x1/sum(p_x1)
p_y1 <- colSums(outer(p_x,p_y)*dnorm(d,outer(x,y,"-"),liklihood_shape))
p_y1 <- p_y1/sum(p_y1)


won't lead you to calculate different parameters...

• But it does contain information and we can update the information we have. Imagine I start with a prior on a pair of variables (W and Z say) with variance L and L respectively (for L very large) and I observe 1 observation on one of the components (W). I can update the distribution of the mean of that component. Now let W be X-Y and Z be X+Y and we have the essentials of the present situation -- we have obtained an observation about a component of a transformed version of the problem. It should not be impossible to incorporate what that information tells us – Glen_b -Reinstate Monica Oct 3 '16 at 7:53
• I'm brought to mind of a paper by Kohn and Ansley (I think) on updating a state space model with initial state that has infinite variances - it does pretty much the sort of updating that I'm thinking of ... I'll see if I can recall the title. It updates the diffuse components as information comes in – Glen_b -Reinstate Monica Oct 3 '16 at 7:55
• Craig F. Ansley and Robert Kohn (1985), $\:$ "Estimation, Filtering, and Smoothing in State Space Models with Incompletely Specified Initial Conditions", $\:$ Ann. Statist., Vol. 13, No. 4, pp 1286-1316. $\:$ projecteuclid.org/euclid.aos/1176349739 $\quad$ (There's also a somewhat later paper by Piet de Jong on a similar topic) – Glen_b -Reinstate Monica Oct 3 '16 at 8:34
• @Glen_b I see your point, that's interesting. When I find some spare time I'll take a look at the papers. Maybe you would be willing to provide answer to this Q? – Tim Oct 3 '16 at 8:46
• @Glen_b I have done basically what you suggested. My bare bones version is just to define a PMF of each variable and then update each PMF after each observation of d observed. In the simple case of just two variables, the posterior means end up being opposite of one another. I can post my R code if anyone would like. – adamwlev Oct 4 '16 at 23:52

I have coded up a start towards answering the question. I assume independence of X and Y (which happens to suit my problem).

My approach was to discretize the two random variables X and Y and then update the probabilities upon observing data. I think this is like Approximate Bayes Computation. I needed a likelihood function which dictates how likely it is to observe a value given a difference in X and Y. I specified this to be a normal RV.

R Code:

require(colorRamps)

x <- seq(-300,300)/25 ##possible values of now discrete RV X
p_x <- dnorm(x,0,10)/sum(dnorm(x,0,10)) ##prior probabilities for X
y <- seq(-300,300)/25 ##possible values of Y
p_y <- dnorm(y,0,10)/sum(dnorm(y,0,10)) ##prior probabilities for Y

## this is the standard deviation of the liklihood function.
## it specifies how much information should be gleaned from
## any single observation of d.
liklihood_shape <- 60

## updates with bayes' formula. calculates posterior for each
## discrete value of X and Y given the liklidhood of the data
## under the joint distribution X and Y (assuming independence).
update <- function(x,p_x,y,p_y,d){
p_x1 <- rowSums(outer(p_x,p_y)*dnorm(d,outer(x,y,"-"),liklihood_shape))
p_x1 <- p_x1/sum(p_x1)
p_y1 <- colSums(outer(p_x,p_y)*dnorm(d,outer(x,y,"-"),liklihood_shape))
p_y1 <- p_y1/sum(p_y1)
return(list(p_x1,p_y1))
}

ds <- c(rnorm(200,-2,7)) ##simulated data
colorsX = blue2red(length(ds)+1)
colorsY = magenta2green(length(ds)+1)
plot(x,p_x,type='l',col=colorsX[1],ylim=c(0,.007)) ##plot the prior of X
plot(y,p_y,type='l',col=colorsY[1],ylim=c(0,.007)) ##plot the prior of Y
for (d in 1:length(ds)){ ##update the data and plot posteriors for each datum
res = update(x,p_x,y,p_y,ds[d])
p_x = res[[1]]
p_y = res[[2]]
lines(x,p_x,col=colorsX[d+1])
lines(y,p_y,col=colorsY[d+1])
}

• The main comment is that it doesn't work. Try defining X and Y as normal variables with equal sd's and different means, say -3 and 7. Then set prior means to the same values. You will notice that the posterior would "converge" to symmetric means that are far from the true/prior means no matter that you started with the correct values. – Tim Oct 5 '16 at 19:06
• @Tim hmm, I don't quite understand what you mean. So set the prior on X and Y to be centered around -3 and 7 respectively. Then do what? – adamwlev Oct 5 '16 at 23:52
• Your code is wrong. I provided extended comment to it in my answer. – Tim Oct 6 '16 at 8:14