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:


*

*Is this a known problem with a solution worked out?

*If not, does anyone have any ideas of how to attack it?

 A: 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) $$
so difference of two normals is also normally distributed.
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...
A: 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])
}

