# Bias from a simulation

I have some model $\vec{d}=A\vec{s}+\vec{n}$

$\vec{s}=(s_1,s_2)^T$ is a vector of two signals with a gaussian pdf prior with a zero mean and a diagonal covariance matrix A is a matrix the mixes the signals $\vec{n}$ is a noise 3 component vector with a gaussian pdf with zero mean and unit matrix covariance d is the data modeled in this way

My aim is given the data, find an estimation for the signal $s_1$ by finding a posterior distribution for it.

So, I used bayes theorem (G is a gaussian with zero mean, C covariance matrix for the components and N for the noise)

$p(\vec{s}|{\vec{d}}) \propto G(\vec{s}, C)G(\vec{n},N)$ Or

$p(\vec{s}|{\vec{d}}) \propto G(\vec{s}, C)G(\vec{d}-A\vec{s},N)$

Then I marginalized it in the component $s_2$ and I found another final pdf for the component $s_1$, the marginalized posterior, say $p(s_1)$ (in this case it is a Gaussian)

From $p(s_1)$ I found the maximum, finding $\hat{s}_1$ that maximeses it.

This $\hat{s}_1$ depends on the vector d and on the covariances.

Now I want to test it:

I generate a fixed value (the TRUE one) from the gaussian for $\vec{s}$ And I generate K noises From the noises I obtain my data as $\vec{d}=A\vec{s}+\vec{n}$

Then I use my estimator to see if this is equal to the true value.

In this plot the x axis is the 'measurement' number i On the y axis there are the estimations $\hat{s}_1$ in red dot and the black line is the true signal $s_1$.

I thought that an estimator should reproduce the true value. The problem is that all the dots are 'far' for the true signal. Also their mean is far. What am I doing wrong?

Have I some conceptual misunderstanding or the issue is another one?

And also:

Now the definition of the bias is that of the expectation value of the estimation - the true signal

But how can I find the expectation value of the estimation? Should I take a lot of data, make a lot of estimations and take the mean?

This is a typical graph

• First things first: do you know the true value of the signal?
– Tim
Sep 6, 2017 at 18:58
• Maybe I didn't understand the theory: but for me the true value,say S, of the signa is the one generated by the prior. Then I fix S, and I generated data in this way: generate a noise Ni (ex. gaussian noise with its prior) Di=S+Ni Then I insert Di in my posterior pdf and find the maximum of the posterior. This is my estimation. Why it isn't equal to the value S generated by the prior? It is because it is biased? Or there is something that I don't understand? Thank you! Sep 6, 2017 at 19:06
• Could you edit your question to describe your procedure in greater detail? How exactly do you conduct your simulation? How exactly do you estimate your posterior? Without those information it is hard to answer.
– Tim
Sep 6, 2017 at 19:17
• Updated! Tell me if I have to add anything else Sep 6, 2017 at 19:38
• I do not understand the "So, I used bayes theorem..." part: this does not seem to be a Bayes theorem. You seem to multiply the distribution of noise times the distribution of signal, instead you should rather multiply distribution of your data times the distributions of your priors for the mean and covariance of it.
– Tim
Sep 6, 2017 at 20:24

Your notation is confusing for me, so allow me to change it slightly in my answer. If your model is (for simplicity let's forget about $A$)

$$D = S + N$$

where $S,N,D$ are multivariate normal random variables. Here your notation becomes confusing, since you seem to use $\vec{s},\vec{n}$ as means of normally distributed random variables and $\vec{d}$ for a random variable itself. This leads to another confusion, you named $\vec{n}$ as "noise", but in fact there is a random noise around means of both random variables, so it could be re-written as

$$D = \mu_s + \varepsilon_s + \mu_n + \varepsilon_n$$

where $\mu_s,\mu_n$ are point values and $\varepsilon_s \sim \mathcal{N}(\boldsymbol{0}, \Sigma_s)$ and $\varepsilon_n \sim \mathcal{N}(\boldsymbol{0}, \Sigma_n)$, what is equivalent to saying that $S \sim \mathcal{N}(\mu_s, \Sigma_s)$ and $N \sim \mathcal{N}(\mu_n, \Sigma_n)$. So in fact you have two sources of noise (randomness) if you assume $S,N$ to be random variables.

Your general mistake is that you are not using Bayes theorem at all in your equations. What you write is

$$p(S|D) \propto p(S) \,p(N)$$

and this is not a Bayes theorem, moreover it is wrong. Distribution of a sum of random variables $S+N$ is not a product of their distributions. If $S$ and $N$ were independent, then their joint distribution would be a product of their marginal distributions. Joint distribution is not the same as a distribution of sum of variables. Moreover, for this to be a Bayes theorem, on the right hand side of the equation we would need to see $p(D|S)\,p(S)$.

Moreover, the idea behind this equation is wrong. Basically, what you are saying is that you want to learn about the distribution of $S$ using the distribution of $S$. If you know it, then what is the point of using it to estimate itself?

The proper use of Bayes theorem would be something like below. First, we would define the likelihood function, i.e. the function of your data conditional on parameters,

$$D \sim \mathcal{N}(\mu_d, \Sigma_d)$$

In your case the parameters are $\mu_d = \mu_s + \mu_n$ and $\Sigma_d = \Sigma_s + \Sigma_n$. Since you want to conduct a simulation to assess your estimator, you would assume that you don't know $\mu_s,\mu_n,\Sigma_s,\Sigma_n$ (otherwise there is nothing to estimate, as already described). The problem in here is that you cannot separate $S$ and $N$ knowing only their sum. Fortunately, if you know that $\mu_n = \boldsymbol{0}$ and $\Sigma_n = \boldsymbol{I}$, then you also know that $\mu_d = \mu_s + \boldsymbol{0} = \mu_s$, what simplifies your problem to estimating $\mu_d$. To estimate $\mu_d,\Sigma_d$ you would need to assume some kind of priors for them and use Bayes theorem

$$p(\mu_d,\Sigma_d \mid D) \propto p(D \mid \mu_d,\Sigma_d) \; p(\mu_d) \; p(\Sigma_d)$$

where $p(\mu_d), p(\Sigma_d)$ are priors for your parameters. To estimate it, you could use optimization (if you need only the MAP estimate), conjugate priors, or MCMC.

Finally, bias of the estimator would be

$$\mathrm{bias}(\hat \mu_d) = \mu_d - E(\hat \mu_d)$$

where $\hat \mu_d$ is your estimate of $\mu_d$. Yes, to estimate $E(\hat \mu_d)$ and the bias you could run a large number of simulations, and then take the empirical mean of the estimated $\hat \mu_d$'s as estimate of $E(\hat \mu_d)$.

• Thank you very much. But I want please to understand what it is the true signal. Can I theoretically set it? Sep 6, 2017 at 23:39
• @Saladino True signal is the point value of the original signal that is not distorted by noise, or any random interference. It depends on your problem. In statistics it is usually an abstract, unobservable quantity unless you conduct the simulation, so it's value is chosen by you and so, known to you.
– Tim
Sep 7, 2017 at 7:31
• @Saladino I edited for clarity.
– Tim
Sep 7, 2017 at 7:43
• So, if I have a prior can I generate the true signal from it and then make an estimation from the posterior? Or is this not correct? Thank you! Sep 7, 2017 at 10:43
• @Saladino please read my answer carefully once again, you can also check some handbook on Bayesian statistics. If you generate data $X$ from distribution $f$ (what you call a "prior"), then what do you want to estimate? You already know that the distribution of $X$ is $f$, so there is nothing to estimate. It does not make sense to estimate $f$ using $f$ -- how doing anything with $f$ would make your estimate better then taking $f$ unchanged?
– Tim
Sep 7, 2017 at 10:52