# The setup

I've tried to reduce the problem to a self-contained subset for this question, but it still ended up being pretty long. Sorry about that!

I have a set of observations of a set of objects where the position of each object is known up to an overall common shift $b$, but the amplitude $a_s$ of each object $s$ is unknown, leading to the data model

$${\bf d} = P\mathbf{a} + \mathbf{n}$$

where $P$ is the profile matrix of each object (including its offset $b$ which $P$ implicitly depends on) and $\mathbf{n} \leftarrow \mathcal{N}(0,N)$. For example, a possible shape of $P$ is $P_{is} = \mathrm{exp}(-\frac12 [p_s+b-x_i]^2)$ where $x_i$ is the location of sample $i$ and $p_s$ is the known position of object $s$, but the exact shape of the profiles is not very important here.

The marginalized probability for $b$ is

$$P(b|\mathbf{d}) = \int d\mathbf{a} P(b,\mathbf{a}|\mathbf{d}) = \int d\mathbf{a} P(\mathbf{d}|b,\mathbf{a})P(b,\mathbf{a})$$

If I assume a uniform prior over $b,\mathbf{a}$, we get

$$P(b|\mathbf{d}) \propto \int d\mathbf{a} P(\mathbf{d}|b,\mathbf{a}) = \frac{1}{\sqrt{|2\pi N|}} \int \mathrm{exp}(-\frac12 [\mathbf{d}-P\mathbf{a}]^TN^{-1}[\mathbf{d}-P\mathbf{a}])d\mathbf{a}$$

By completing the square in $\mathbf{a}$, we can rewrite the term in the exponent as

$$[\mathbf{d}-P\mathbf{a}]^TN^{-1}[\mathbf{d}-P\mathbf{a}] = [\mathbf{a}-\mathbf{\hat a}]^TA^{-1}[\mathbf{a}-\mathbf{\hat a}] + \mathbf{d}^TN^{-1}\mathbf{d} - \mathbf{\hat a}^TA^{-1}\mathbf{\hat a}$$

with $\mathbf{\hat a} = (P^TN^{-1}P)^{-1}P^TN^{-1}\mathbf{d}$ and $A = (P^TN^{-1}P)^{-1}$. With this, the integral is simply

$$\int \mathrm{exp}(-\frac12 [\mathbf{d}-P\mathbf{a}]^TN^{-1}[\mathbf{d}-P\mathbf{a}])d\mathbf{a} = \sqrt{|2\pi A|}\mathrm{exp}(\frac12\mathbf{\hat a}^TA^{-1}\mathbf{\hat a})\mathrm{exp}(-\frac12\mathbf{d}^TN^{-1}\mathbf{d}) \propto \sqrt{|2\pi A|}\mathrm{exp}(\frac12\mathbf{\hat a}^TA^{-1}\mathbf{\hat a})$$

So

$$P(b|\mathbf{d}) \propto \sqrt{|2\pi A|}\mathrm{exp}(\frac12\mathbf{\hat a}^TA^{-1}\mathbf{\hat a})$$

with the actual $b$ dependence hidden in $A$ and $\mathbf{\hat a}$, and where we can recover the normalization by requiring that $\int P(b|\mathbf{d})db=1$.

# The problem

So far so good. But now consider what happens if we have a large number of objects, but very low signal-to-noise. The number in the exponent will then be chisquare-distributed with $N_\textrm{obj}$ degrees of freedom and hence have a standard deviation of $\sqrt{\frac{N_\textrm{obj}}{2}}$ (including the $\frac12$ that was already there). Hence we can expect the probability to fluctuate as

$$P_\pm(b|\mathbf{d}) \propto \sqrt{|2\pi A_\pm|} \mathrm{exp}\left(\pm \sqrt{\frac{N_\textrm{obj}}{2}}\right)$$

Hence due to simple noise fluctuations, one must expect the marginal probability for $b$ to fluctuate by at least about

$$\frac{P(b|\mathbf{d})}{P(b_\textrm{max}|\mathbf{d})} \lesssim \textrm{exp}(-\sqrt{2N_\textrm{obj}})$$

(I've neglected the $\sqrt{|2\pi A|}$ term here, which does not qualitatively change the conclusion for evenly distributed data points).

For example, with $N_\textrm{obj}=100$ objects, one typically finds that large regions of the parameter space for $b$ are disfavored by a factor of about $\textrm{exp}(-\sqrt{2\cdot100})\approx\textrm{exp}(-14)$ compared to the maximum posterior value, even in noise-only simulations. Which exact regions are disfavored vary based on the noise realization, and on average no particular region is favored (for the case of even sample density, at least), but given a data set seems unreasonably certain about where the objects should be.

Here is a plot illustrating this, for 5 different noise-only realizations, where it really should be impossible to determine a significant value for $b$:

These symptoms are similar to those of overfitting, and reducing the number of amplitudes in the model reduces the problem. But in this case I am explicitly marginalizing over all the extra parameters, so I don't see how this could be overfitting.

The problem can also be reduced by imposing a tighter prior on the amplitudes $\mathbf{a}$. But there is a quite large physically reasonable parameter range for $\mathbf{a}$ to occupy, so I do not think this can be a real solution to the problem. I have not explored very far in this direction though.

# The question

So my question is: Why am I getting this excessive significance? It seems to follow inexorably from my model and assumptions, but they all seem reasonable to me. How can I modify my model/assumptions/etc. such that a one wouldn't end up claiming a strong detection in a signal-less simulation?

• What do you mean by a "noise-only realization"? Are you drawing from the model? – Tom Minka Nov 26 '14 at 12:20
• Tom Minka: For that figure I'm drawing from the model in the limit where there is only noise, i.e. $\mathbf{d}\approx \mathbf{n}$, where $\mathbf{n}\leftarrow \mathcal{N}(0,N)$. – amaurea Nov 26 '14 at 13:45
• You can expect Bayes to give you the best possible answer under the assumptions you've built into the model. Unless the full model nests the noise-only model within it such that posterior mass can concentrate on the noise-only possibility, Bayes will hallucinate, just as you've told it to. – Cyan Dec 2 '14 at 3:16
• To get the noise-only model onto the stage, I guess what you'd have to do is make $N_\rm{obj}$ a model parameter and do inference on it. You'll need proper priors over the object-specific parameters; otherwise, you'll run into the Jeffreys-Lindley thingy. – Cyan Dec 2 '14 at 3:27

As Cyan pointed out, this behaviour arises from the fact that the model does not put enough emphasis on the case $\mathbf{a}=0$. You say that you don't want to impose a 'tighter' prior on $\mathbf{a}$. I'm not sure what you mean by 'tighter', but increasing the prior probability of $\mathbf{a}=0$ is the only solution. This could be done by putting a point mass on $\mathbf{a}=0$, i.e. model averaging as suggested by Cyan, or using a prior with a sharp peak at 0 such as a T distribution with small $\nu$.