# Uncertainties on Bayes' theorem

• A
kelly0303
Hello! I have a setup consisting of some charged particles each of which is produced at a different position, ##(x_i,y_i,z_i)##. I don't know the exact position, but I know that each of the 3 variables is normally distributed with mean zero and standard deviation of 3 mm. What I measure in the experiment is the position of the particles on a 2D screen perpendicular to the trajectory of motion of the particles, call it x' and y' and the time from the moment they are produced to when they hit the 2D screen, call it t. So what i measure is ##(x',y',t)## with the associated uncertainties, given by the detector, ##(dx',dy',dt)##. What I want to do is, for each ion, from these measurements, to get the most probable starting point (together with the associated uncertainty). I am not sure what is the best way to proceed. What I want to get basically is:

$$p(x_i,y_i,z_i|x_f,y_f,t)$$
and I will write this just as ##p(x_i|x_f)## for simplicity. For example, assuming we had no uncertainty on our measurements, I could use Bayes theorem and have:

$$p(x_i|x_f) = \frac{p(x_f|x_i)p(x_i)}{p(x_f)}$$

I know ##p(x_i)## (which is Gaussian) and for a given event, ##p(x_f)## is constant (I don't know what it is, but it is constant). For ##p(x_f|x_i)## I can get a value by running some simulations. Assuming everything would be ideal, this would be 1 for just one value and 0 for everything else. Given that we can't model the electrodes perfectly, and I will run the simulation several times with slightly different potentials, this will probably have a distribution, too (say Gaussian for simplicity). From here I can get ##p(x_i|x_f)## and the associated standard deviation, which is exactly what I need.

However, given that I have noise on ##(x_f,y_f,t)##, I can't just apply the Bayes theorem directly as above. How can I account for these errors? Thank you!

EDIT: I attached before a diagram of the trajectory of the particle (red line) from the source to the 2D detector.

Last edited:

I don't know the exact position, but I know that each of the 3 variables is normally distributed with mean zero and standard deviation of 3 mm.

If that is the model, then the final position of a particle is distributed 3-D normally about the point (0,0,0). Is that the case?

Or perhaps you mean that the ##i##th particle is produced at some unknown location ##(x_i,y_i,z_i)## and travels to position ##(x_i+\triangle x, y_i + \triangle y, z_i + \triangle z)## where each of ##\triangle x, \triangle y , \triangle z## is distributed normally with mean zero and standard deviation 3 mm.

kelly0303
If that is the model, then the final position of a particle is distributed 3-D normally about the point (0,0,0). Is that the case?

Or perhaps you mean that the ##i##th particle is produced at some unknown location ##(x_i,y_i,z_i)## and travels to position ##(x_i+\triangle x, y_i + \triangle y, z_i + \triangle z)## where each of ##\triangle x, \triangle y , \triangle z## is distributed normally with mean zero and standard deviation 3 mm.
Hello! The initial position is distributed 3-D normally about the point (0,0,0). The final position is the one I measure.

Staff Emeritus
Gold Member
I have a couple questions

1.) Is there any uncertainty about the angle of motion, or is it 100% guaranteed to be in the z direction?

2.) Is the uncertainty in the detector also guassian?

3.) What kind of uncertainty do you have on the speed of the particle?

DaveE
kelly0303
I have a couple questions

1.) Is there any uncertainty about the angle of motion, or is it 100% guaranteed to be in the z direction?

2.) Is the uncertainty in the detector also guassian?

3.) What kind of uncertainty do you have on the speed of the particle?
1) The motion is not in the z direction. As I mentioned the x' and y' of the detector are in a different plane compared to the starting point (I will add a simulation picture in the original post)

2) Yes

3) The uncertainty on the energy is also gaussian, but for simplicity we can assume that the energy is constant for all particles.

Staff Emeritus
Gold Member
Got it. Does the uncertainty in the detector measurement include uncertainty about how literally the particle follows the red line that you've drawn? I'm a little nervous that if the tube is long, very small deflections in the path of the particle could result in it striking a place on the screen at the end that is arbitrarily far from where you would "expect" it to hit.

kelly0303
Got it. Does the uncertainty in the detector measurement include uncertainty about how literally the particle follows the red line that you've drawn? I'm a little nervous that if the tube is long, very small deflections in the path of the particle could result in it striking a place on the screen at the end that is arbitrarily far from where you would "expect" it to hit.
In practice we have a very good control on the voltage on the potentials (I can check the actual relative uncertainty). But of course there will be small variations. This is why I mentioned that in practice I would slightly change the potential in the simulations to account for that. So basically, for a given initial point (that we don't know but it's fixed), the final point on the detector will have a distribution given by this variation in the voltage, even if the detector would have infinite resolution. In practice, on top of that, we also have the resolution of the detector which gives this ##(dx',dy',dt)##.

Gold Member
Am I right to think that you're kicking a Coulomb crystal out of a Paul trap? Sounds like you've got an MCP detector with a phosphor plate for imaging? I think I can help, but my internet connection is on the rocks right now thanks to black friday I'll write a post with the maths once I'm sure it won't suddenly disappear into the void thanks to a sudden disconnect

Actually, you can propagate the uncertainties on the voltages without needing to use a simulation, if your trap geometry is nice and simple. Can we assume a standard quadrupole trap with rod-shaped electrodes? And can you quote for us the fractional uncertainty on the pulsed voltages?

Last thing, I'm going to assume that the relative timing of your pulsed voltages is spot-on. If it isn't that can be another source of uncertainty.

Gold Member
Ok, sorry, I read the OP a little more carefully and I think understand better what you're asking for. My bad.

If you have uncertainties on your voltage, for example, you can still use Bayes theorem if you add the voltage in as a third variable. So instead of ##p(x_i|x_f,t)##, you'll want to use the simulation to get ##p(x_i|x_f,t,V)##. Then integrate: $$p(x_i|x_f,t) = \int p(x_i|x_f,t,V) p(V) dV$$
So if V is normal distributed with uncertainty ##\sigma_V## around a mean value ##\langle V \rangle##, then
$$p(x_i|x_f,t)= \int p(x_i|x_f,t,V) \frac{1}{\sigma_V \sqrt{2\pi}} \exp\left(-\frac{(V - \langle V \rangle)^2}{2\sigma_V^2} \right) dV$$

kelly0303
Ok, sorry, I read the OP a little more carefully and I think understand better what you're asking for. My bad.

If you have uncertainties on your voltage, for example, you can still use Bayes theorem if you add the voltage in as a third variable. So instead of ##p(x_i|x_f,t)##, you'll want to use the simulation to get ##p(x_i|x_f,t,V)##. Then integrate: $$p(x_i|x_f,t) = \int p(x_i|x_f,t,V) p(V) dV$$
So if V is normal distributed with uncertainty ##\sigma_V## around a mean value ##\langle V \rangle##, then
$$p(x_i|x_f,t)= \int p(x_i|x_f,t,V) \frac{1}{\sigma_V \sqrt{2\pi}} \exp\left(-\frac{(V - \langle V \rangle)^2}{2\sigma_V^2} \right) dV$$
Thank you! But I am still not sure how to use this for my end goal. For now, you can even assume there is no uncertainty on the voltages. The only uncertainty is on ##x_f## and ##t##. How do I account for these uncertainties when using Bayes theorem?

Gold Member
Sorry! Looks like I missed the point again

You need to measure ##p(x_f,t)## somehow. Equivalently, since you assume it's gaussian you just need to know ##\sigma_{x_f}## and ##\sigma_t##.

If your detector is a MCP and phosphor plate plus a CCD/CMOS camera, then ##\sigma_{x_f}## is bounded below by the distance between pores in the MCP (times the magnification of your camera lens, if you haven't factored that out already).

Question for you. Is there a source of uncertainty on the arrival time ##t## besides the uncertainty in the initial z-position of the ion? I'm choosing the z-axis here to be the direction the ions get kicked. If not, then I would just look at the MCP voltage vs time for a single ion and set ##\sigma_t## equal to the width of that single ion pulse.

Edit: I just re-read the OP and realized you wrote ##(x_f,y_f,t)## and not ##(x_f,y_f,z_f,t)##. Ignore the crossed-out part. Sorry, but I don't have a general method for you to find ##\sigma_t## in this case because it strongly depends on your experiment. We'd need more information: how are the ions being trapped (is it a standard type of ion trap)? how cold are they? how dense are they? and what's the trap frequency and mass of the ions? Ballpark numbers are fine!

Edit: Forgot to add some things

Once you have ##p(x_f,t)##, you can find the probability distribution on the initial position as follows: $$p(x_i) = \int \int p(x_i | x_f,t) p(x_f,t) dx_f dt$$

Last edited:
kelly0303
Sorry! Looks like I missed the point again

You need to measure ##p(x_f,t)## somehow. Equivalently, since you assume it's gaussian you just need to know ##\sigma_{x_f}## and ##\sigma_t##.

If your detector is a MCP and phosphor plate plus a CCD/CMOS camera, then ##\sigma_{x_f}## is bounded below by the distance between pores in the MCP (times the magnification of your camera lens, if you haven't factored that out already).

Question for you. Is there a source of uncertainty on the arrival time ##t## besides the uncertainty in the initial z-position of the ion? I'm choosing the z-axis here to be the direction the ions get kicked. If not, then I would just look at the MCP voltage vs time for a single ion and set ##\sigma_t## equal to the width of that single ion pulse.

Edit: I just re-read the OP and realized you wrote ##(x_f,y_f,t)## and not ##(x_f,y_f,z_f,t)##. Ignore the crossed-out part. Sorry, but I don't have a general method for you to find ##\sigma_t## in this case because it strongly depends on your experiment. We'd need more information: how are the ions being trapped (is it a standard type of ion trap)? how cold are they? how dense are they? and what's the trap frequency and mass of the ions? Ballpark numbers are fine!

Edit: Forgot to add some things

Once you have ##p(x_f,t)##, you can find the probability distribution on the initial position as follows: $$p(x_i) = \int \int p(x_i | x_f,t) p(x_f,t) dx_f dt$$
Hello! So I actually don't have an ion trap. I have some atoms traveling and I ionize them (perpendicularly) at some point, ##(x_i,y_i,z_i)##. This initial positions are basically given by the distribution of the ion beam cross section and the laser beam (which is Gaussian with diameter of 3mm). At this ionization point, I collect the electrons (for now one can assume I do this one atom at a time), which give me ##t_i## and ##dt_i##. The ion travels to the MCP detector, and when it hits, I get the 2D position on the detector (together with the associated uncertainty) and the time ##t_f## with ##dt_f##. In practice the uncertainty in the ##t_i## is given by the direction in which the electrons are emitted, and from simulations it is about 2ns. For ##t_f## the uncertainty is the MCP resolution which is 100ps, so basically one can assume the uncertainty in the time traveled is a gaussian with a standard deviation of 2ns. The MCP position resolution is 100 micros. I also know the average energy of the atoms (10keV) and for the energy spread we can assume 10eV. So knowing the initial position to the accuracy mentioned above, the final position, the time it takes to travel (and the uncertainties, ignore the voltages for now) and also the distribution of the kinetic energies, what I really want to do is to extract the energy of each individual atom, which hopefully would be smaller than 10eV (this is a lot more detail than the original post, I tried to simplify it in the beginning).

Gold Member
So I actually don't have an ion trap. I have some atoms traveling and I ionize them (perpendicularly) at some point, (xi,yi,zi).
Haha I struck out again. I should just keep my mouth shut more

Ok! With ##dt_f = 2\mathrm{ ns}## and ##dx_f = 100\mathrm{ \mu m}## from the MCP resolution, you have everything you need! I'm going to call these ##\sigma_{t_f}## and ##\sigma_{x_f}## for clarity. You can now construct the probability distribution on ##(x_f,t_f)##: $$p(x_f,t_f) = \frac{1}{\sigma_{t_f} \sqrt{2\pi}} \exp\left(-\frac{(t_f - \langle t_f \rangle)^2}{2\sigma_{t_f}^2} \right) \times \frac{1}{\sigma_{x_f} \sqrt{2\pi}} \exp\left(-\frac{(x_f - \langle x_f \rangle)^2}{2\sigma_{x_f}^2} \right)$$ This is of course assuming that ##x_f## and ##t_f## are uncorrelated, which I believe is a fair assumption for the transverse spatial coordinates.

Now you can run the simulations to find ##p(x_i | x_f,t)##, and finally you can plug into $$p(x_i) = \int \int p(x_i | x_f,t) p(x_f,t) dx_f dt$$ By the way, I never really explained this formula. If it isn't clear, just think of it as doing ##p(A) = \sum_{B} p(A \cap B) = \sum_B p(A|B) p(B)##.

There's a caveat here, and that's because the uncertainty in the final position might be larger than the 100 micron resolution of the MCP detector. The most annoying contributor would be static charges on surfaces in the time-of-flight (TOF for short) path, pushing the electrons around in the transverse direction. There's also Coulomb repulsion between ions at large densities. This actually broadens the ions in both space (spreading out) and time (getting "longer"). (Incidentally, this Coulomb time-broadening effect was one of the major hurdles in the development of ultrafast electron microscopes (UEM), which got a Nobel Prize in the late '90s, but anyways I'm getting distracted.) You'd have to characterize these spreading effects independently. Ideally, you'd want to create ions at the sample place every time and see where they end up, but I realize you have very little control over where the atoms get ionized.

kelly0303
Haha I struck out again. I should just keep my mouth shut more

Ok! With ##dt_f = 2\mathrm{ ns}## and ##dx_f = 100\mathrm{ \mu m}## from the MCP resolution, you have everything you need! I'm going to call these ##\sigma_{t_f}## and ##\sigma_{x_f}## for clarity. You can now construct the probability distribution on ##(x_f,t_f)##: $$p(x_f,t_f) = \frac{1}{\sigma_{t_f} \sqrt{2\pi}} \exp\left(-\frac{(t_f - \langle t_f \rangle)^2}{2\sigma_{t_f}^2} \right) \times \frac{1}{\sigma_{x_f} \sqrt{2\pi}} \exp\left(-\frac{(x_f - \langle x_f \rangle)^2}{2\sigma_{x_f}^2} \right)$$ This is of course assuming that ##x_f## and ##t_f## are uncorrelated, which I believe is a fair assumption for the transverse spatial coordinates.

Now you can run the simulations to find ##p(x_i | x_f,t)##, and finally you can plug into $$p(x_i) = \int \int p(x_i | x_f,t) p(x_f,t) dx_f dt$$ By the way, I never really explained this formula. If it isn't clear, just think of it as doing ##p(A) = \sum_{B} p(A \cap B) = \sum_B p(A|B) p(B)##.

There's a caveat here, and that's because the uncertainty in the final position might be larger than the 100 micron resolution of the MCP detector. The most annoying contributor would be static charges on surfaces in the time-of-flight (TOF for short) path, pushing the electrons around in the transverse direction. There's also Coulomb repulsion between ions at large densities. This actually broadens the ions in both space (spreading out) and time (getting "longer"). (Incidentally, this Coulomb time-broadening effect was one of the major hurdles in the development of ultrafast electron microscopes (UEM), which got a Nobel Prize in the late '90s, but anyways I'm getting distracted.) You'd have to characterize these spreading effects independently. Ideally, you'd want to create ions at the sample place every time and see where they end up, but I realize you have very little control over where the atoms get ionized.
Thanks a lot for this! I understand the formulas, but I am still not totally sure this is the whole story for my situation (I am sorry if I am missing something). I am not sure how would I simulate ##p(x_i|x_f,t)##. In my simulations (for reference I am using SIMION), what I can simulate is ##p(x_f,t|x_i)## i.e. I can set the starting point and then record the distribution of the ions on the final screen. But ##p(x_i|x_f,t)## is some sort of inverse problem, which unless I am missing something, is not something one can simulate easily.

Also, I am not sure I understand that formula for ##p(x_i)##. I do understand it mathematically, but I am not sure why in my situation I would need to sum over all ##x_f##? What I need is ##p(x_i)##, given the ##x_f## that I measured. I think what I need is ##p(x_i|x_f,t)##, but again, I don't think I can simulate that easily.

Gold Member
I am not sure how would I simulate p(xi|xf,t). In my simulations (for reference I am using SIMION), what I can simulate is p(xf,t|xi) i.e. I can set the starting point and then record the distribution of the ions on the final screen.
Ah got it. Yes you're right there's an inverse problem that needs to be solved. I think I have a solution, it just might not be the cleanest way to implement it.

For now, ignoring voltage fluctuations, you can simulate the deterministic ion trajectories in SIMION. This should give you a table of results: outputs ##x_f## vs inputs ##(x_i,t)##. Once you have enough data points in this table, you can do least-squares regression to a model of the form $$x_i = F(x_f,t)$$ where ##F(x_f,t)## can be any model function you like (piece-wise linear and polynomial should be applicable). Then you can write $$p(x_i | x_f,t) = \delta (x_i - F(x_f,t) )$$

For every fluctuating quantity you take into account (like electrode voltages), add another column to the table. For example, a table with outputs ##x_f## vs inputs ##(x_i, t, V)## would be fitted to a model like ##x_i = F(x_f,t,V)##, and the conditional probability distribution would be ##p(x_i | x_f,t,V) = \delta (x_i - F(x_f,t,V) )##. Then you would integrate $$p(x_i|x_f,t) = \int \frac{1}{\sigma_V \sqrt{2\pi}} \exp\left(-\frac{(V - \langle V \rangle)^2}{2\sigma_V^2} \right) \delta (x_i - F(x_f,t,V)) dV$$

I'm sure that mathematically, this holds water, but it might not be very pretty or easy to implement. Especially that last integral...

sysprog
kelly0303
Ah got it. Yes you're right there's an inverse problem that needs to be solved. I think I have a solution, it just might not be the cleanest way to implement it.

For now, ignoring voltage fluctuations, you can simulate the deterministic ion trajectories in SIMION. This should give you a table of results: outputs ##x_f## vs inputs ##(x_i,t)##. Once you have enough data points in this table, you can do least-squares regression to a model of the form $$x_i = F(x_f,t)$$ where ##F(x_f,t)## can be any model function you like (piece-wise linear and polynomial should be applicable). Then you can write $$p(x_i | x_f,t) = \delta (x_i - F(x_f,t) )$$

For every fluctuating quantity you take into account (like electrode voltages), add another column to the table. For example, a table with outputs ##x_f## vs inputs ##(x_i, t, V)## would be fitted to a model like ##x_i = F(x_f,t,V)##, and the conditional probability distribution would be ##p(x_i | x_f,t,V) = \delta (x_i - F(x_f,t,V) )##. Then you would integrate $$p(x_i|x_f,t) = \int \frac{1}{\sigma_V \sqrt{2\pi}} \exp\left(-\frac{(V - \langle V \rangle)^2}{2\sigma_V^2} \right) \delta (x_i - F(x_f,t,V)) dV$$

I'm sure that mathematically, this holds water, but it might not be very pretty or easy to implement. Especially that last integral...
Thanks for this idea! About fitting that function analytically, that function is a function (in the simplest case) of 3 variables (##(x_f,y_f,t)##), with 3 outputs (##(x_i,y_i,z_i)##). Also the variables are not really independent I think (for example if I change only ##x_i##, it won't be just ##x_f## that gets affected, but all variables). Is that easy to model that function (I might try to use a NN for that, tho, and not have to worry about the analytical form maybe)?

Also, where exactly does the position and time uncertainties come into play in this equation? I can't really integrate over ##x_f## as I do for the voltages, as I have a fixed (measured) ##x_f##, but shouldn't I account for the uncertainty on that value somehow?

sysprog
Gold Member
Is that easy to model that function (I might try to use a NN for that, tho, and not have to worry about the analytical form maybe)?
That's why I suggested piecewise linear or polynomial. You will need more data points, but you'll get a good empirical model. I have a funny feeling that the polynomial fit would be quite good. Just a hunch. I think an NN would be overkill.

Also, where exactly does the position and time uncertainties come into play in this equation? I can't really integrate over xf as I do for the voltages, as I have a fixed (measured) xf, but shouldn't I account for the uncertainty on that value somehow?
You don't account for the uncertainty in ##x_f## in ##p(x_i|x_f,t)##. Rather, you account for the uncertainty in ##p(x_f)##. That's what I was suggesting in the first and second large equations in post #13 (post numbers are in the top right of each post). Did that part not make sense?

I'm thinking that using conditional probability distributions may also be overkill for this problem. If you fit an empirical model ##(x_i,y_i,z_i) = F(x_f,y_f,t)##, then you can just do linear error propagation on ##x_f##, ##y_f##, and ##t##, where F is a vector-valued function (##F = (F_x,F_y,F_z)##). The errors on ##x_f##, ##y_f##, and ##t## are all uncorrelated, right? If that's true, then all you need is just the standard error propagation formula: $$\mathrm{Var}(x_i) = \left(\frac{\partial F_x}{\partial x_f} \right)^2 \mathrm{Var}(x_f) + \left(\frac{\partial F_x}{\partial y_f} \right)^2 \mathrm{Var}(y_f) + \left(\frac{\partial F_x}{\partial t} \right)^2 \mathrm{Var}(t)$$ etc.

kelly0303 and sysprog
kelly0303
That's why I suggested piecewise linear or polynomial. You will need more data points, but you'll get a good empirical model. I have a funny feeling that the polynomial fit would be quite good. Just a hunch. I think an NN would be overkill.

You don't account for the uncertainty in ##x_f## in ##p(x_i|x_f,t)##. Rather, you account for the uncertainty in ##p(x_f)##. That's what I was suggesting in the first and second large equations in post #13 (post numbers are in the top right of each post). Did that part not make sense?

I'm thinking that using conditional probability distributions may also be overkill for this problem. If you fit an empirical model ##(x_i,y_i,z_i) = F(x_f,y_f,t)##, then you can just do linear error propagation on ##x_f##, ##y_f##, and ##t##, where F is a vector-valued function (##F = (F_x,F_y,F_z)##). The errors on ##x_f##, ##y_f##, and ##t## are all uncorrelated, right? If that's true, then all you need is just the standard error propagation formula: $$\mathrm{Var}(x_i) = \left(\frac{\partial F_x}{\partial x_f} \right)^2 \mathrm{Var}(x_f) + \left(\frac{\partial F_x}{\partial y_f} \right)^2 \mathrm{Var}(y_f) + \left(\frac{\partial F_x}{\partial t} \right)^2 \mathrm{Var}(t)$$ etc.
Hello @Twigg, sorry for the late reply. I looked a bit at this and there is still something I am confused about. For now let's ignore any uncertainties (assume we just run a simulation). I realized that what I need in practice is the initial energy of the ion (this is what allows me to perform the Doppler correction). If I understand what you said above, once I have enough simulated data I would fit a function of the form:

$$E_i = F_1(x_f,y_f,t_f)$$

However, this seems to completely ignore any information about the ##z_i## (assume the particle travels along the z axis). Basically a particle with less energy produced further away from the target can end up having the same time of flight as a particle with more energy produced closer to the target. In the approach above (where ##F_1## would be obtained from some polynomial fitting) ##E_i## might not be uniquely determined. I guess my confusion is in treating the initial variables independently (i.e. fit a different function for each one of them), if I understand your argument well. Shouldn't we somehow account for the fact that ##E_i## and ##z_i## are correlated somehow?

Gold Member
You don't have enough information to untangle the energy and initial z position from the time-of-flight. The system is algebraically underconstrained.

If you know about the variance of the initial z position, then you can deduce how much it impacts your certainty on E. The math is the same old error propagation as before. This reflects your loss of precision due to E being correlated with the random variable z.

Set your coordinate system so ##z_i## is zero at the position of the MCP detector. Then, ##v = z_i / t##. This gives you $$\sigma_v^2 = \frac{1}{t^2} \left[\sigma_{z_i}^2 + \frac{z_i^2}{t^2} \sigma_t^2 \right]$$
Depending on the conditions of your experiment, you can estimate ##\sigma_{z_i}##. Say you have a trap frequency of ##\omega_z## and an initial temperature of ##T_i##. Then you can say $$\frac{1}{2} k_B T_i = \frac{1}{2} m \omega_z^2 \sigma_{z_i}^2$$ And thus, $$\sigma_v^2 = \frac{1}{t^2} \left[\frac{k_B T_i}{m \omega_z^2} + \frac{z_i^2}{t^2} \sigma_t^2 \right]$$
Since ##E = \frac{1}{2} mv^2##, you can say $$\sigma_E^2 = \frac{2mE}{t^2} \left[\frac{k_B T_i}{m \omega_z^2} + \frac{z_i^2}{t^2} \sigma_t^2 \right]$$ The main point is that your estimated uncertainty on ##z_i## increases your total uncertainty on ##E## beyond the contribution from your measured uncertainty from the time-of-flight (##t##).