0
$\begingroup$

I have a real world question but I think it'll be clearer if I phrase it like a "homework" question so here goes*:

A physicist is trying to use a laser measuring device to determine the length of a piece of string with high accuracy. Their device has a nominal resolution of $\pm 10 \mu m$, which the user manual goes on to explain refers to the edges of a triangular probability distribution of where the true value the device reads might lie. That is the readings it outputs are in 10 micron increments and a value of $1,350 \mu m$ actually represents a triangular probability density function centred around $1,350$ with width $\pm 10$ (see diagram below) diagram showing triangular distribution

The physicist also notices though that the device is inconsistent; when they try to calibrate it by measuring a known length of string, there is some seemingly random noise in the readings. They decide that by taking a few samples of known lengths and seeing how far off the device is, they can find a probability distribution for the typical error in the device and subtract this from the readings to calibrate it.

The experimental setup is simple. For a known length calibration piece of $1000\mu m$, several measurements are taken. These might be as below:

Measurement No Value ($\mu m$)
1 1050
2 1010
3 980
4 990
5 1010
6 1000
7 1040
8 980

... so clearly there is some rounding owing to the instrument resolution, and some uncertainty from other factors. Using this information, the physicist would like to model the error as a Normal distribution $Err\sim\cal N\left( {\mu,\sigma^2} \right)$ and calculate the most likely values for those parameters...

* The caveat for asking this in homework style is I may have missed out some important details required to give an answer, so if you need clarification or more details please ask

My approach

Ok so I have 2 approaches:

  1. At first I thought of fitting a normal distribution to the data. This is done by first plotting each measurement and its triangular resolution pdf. I then find $\mu$ (which is just $\bar{x}$, easy to calculate since the resolution error is symmetric). Finally I find the value of $\sigma$ which maximises the overlap between the measurements and the resultant distribution with integration - i.e. maximises the likelihood of observing those measurements weighted by their resolution distributions.
  2. However I remembered that samples from a normal population are t-distributed - appearing closer to the centre than expected and having a smaller standard deviation than the population would (making the first approach inaccurate). To solve this, I could try and fit a Student t-distribution to my samples, however I wasn't sure how many degrees of freedom to use since my measurements are pdfs, not individual values.

So I'm a bit stumped on how to fit to a t-distribution or what adjustments I need (given the small sample count and large relative resolution error, $\sigma$ I expect will be larger than just the sample standard deviation of the data itself).

$\endgroup$
4
  • 1
    $\begingroup$ So a few clarifying questions. (1) What does a histogram of the means of the calibration data look like? Triangle-like? Gaussian? Other? (2) What are the sample mean and sample variance? (And how do they compare with the expected calibration mean and variance of the triangle distribution) (3) How precisely do you want to characterize this uncertainty, keeping in mind "effort" scales pretty steeply with "precision", on a scale of "a few simple statistics that give me a ballpark idea" to "a full Bayesian analysis producing a posterior distribution for parameters in a model of the uncertainty"? $\endgroup$
    – Andrew
    Aug 2 '21 at 22:22
  • $\begingroup$ @Andrew good questions, thanks:) 1) The data looks Gaussian but a bit flatter when you account for the resolution error smudging things out 2) Sample mean has a constant offset from the calibration length. Sample standard deviation is a bit complicated to calculate with the resolution taken into account but I can have a think. 3) IRL this will be done numerically on a computer (my low resolution device to calibrate is a code timer API, not a laser) so full Bayesian theoretical precision is better since once I've coded it up I may as well be as comprehensive as possible. $\endgroup$
    – Greedo
    Aug 2 '21 at 22:49
  • $\begingroup$ OK sounds good. I can try to write up some ideas for how to set up a Bayesian analysis. For (2) though I literally meant the sample standard deviation of the means, not accounting for the spread. Eg -- is this comparable to the width of the triangle distribution, or larger? $\endgroup$
    – Andrew
    Aug 2 '21 at 23:30
  • $\begingroup$ @Andrew Ah I see. Yes it is supposed to be comparable - in the data I gave for example it is $26 \mu m$ (compared with the $10 \mu m$ spread from the resolution). $\endgroup$
    – Greedo
    Aug 3 '21 at 9:40
1
$\begingroup$

First things first, in the comments you mentioned there is an offset from the sample mean and the calibration length, implying there is a bias in your measurements. This is very important to understand -- for example, can it be corrected by subtracting a constant value, or by rescaling by a certain factor? My suggestion is to try different calibration lengths, and see how the bias depends on the length of the string. A subtlety here is that you need a way to know the length of your calibration string with at least as much precision as your measuring device, but without using the measuring device you're trying to calibrate; I trust you have a way to do this. Below, I'm assuming you have a way to correct the bias, and apply it to $d_i$. You can also extend all the analysis below to include a "calibration parameter", say by rescaling the data in the likelihood; I can add some comments if you want.

In terms of the error analysis, in the comments you mentioned a Bayesian analysis would be worth doing -- as I mentioned this is on the upper end of the "effort" scale (so this answer is pretty involved!), but the net result will be a flexible framework to get a posterior probability distribution for the length $\ell$, able to incorporate multiple kinds of uncertainty. Before getting stuck in, I'll just mention some simpler (but less rigorous) approaches would be to: (a) look at the standard deviation of the means of the measurements you've done to get a "calibration error" that you can add in quadrature to other uncertainties, (b) fit a Gaussian directly to your observed data with the mean and standard deviation as parameters (this is your option 1 above, and also the first Bayesian model I consider below, but you don't need an MCMC sampler to do this if that's all you need).

If the parameter space is not too large, it's not too much work to code up something that will evaluate the likelihood on a grid and compute/plot the posterior (as well as desired credible intervals), in your favorite programming language (I am partial to python). Alternatively, if you don't want to do that or if you have more than 2 or 3 parameters, it's better to use a sampler. There are several off-the-shelf samplers that are pretty good: emcee is a good MCMC sampler in python, and dynesty is a nested sampler. If you don't know much about samplers, it's probably better to stick with an MCMC sampler to start with, since the way they work is easier to understand and their output is easier to interpret (in my opinion), and also practically speaking it's possible to get intermediate results while the sampler is running (unlike nested sampling).

The main ingredients to a Bayesian analysis are the likelihood and prior.

The likelihood is the probability of the data given the model, $p(d_i|\Theta)$, where $d_i = \{d_1, d_2, \cdots, d_N\}$ are the data (in your case, a list of $N$ observed lengths) and $\Theta$ are the model parameters (we'll consider a few options). My philosophy is to start with a simple model and make it more complicated if it fails. A simple starting point would be to say the combination of the "resolution error" and other errors can be described with a Gaussian distribution: \begin{equation} p_{\rm G.L.}(d_i|\ell,\sigma) = \prod_{i=1}^N \frac{1}{\sqrt{2\pi \sigma^2}}e^{-(d_i-\ell)^2/2\sigma^2} \end{equation} where "G.L." stands for "Gaussian likelihood" (to be compared with a few other options below), and where the model parameters are

  • $\ell$: the length of the string -- this is what you want to fit to the data). I am assuming all the data are measurements of the same string, so note that $\ell$ does not depend on the sample number, $i$.
  • $\sigma$: the uncertainty -- in our simple model, all uncertainty (the "triangle" resolution, and all other sources of uncertainty) will be folded into this parameter. Again, I'm assuming the uncertainty for each measurement is the same. I'm also assuming the measurements are independent, so there is no covariance.

Regardless of whether or not you decide to make the model more complicated, I would start off by implementing this Gaussian model in your sampler of choice, because it is simple and you want to get your code working on a simple model before extending it.

Another simple model would be to say that the noise comes entirely from the resolution (triangle distribution) error. Again, I think this model is worth coding up even it it is not a perfect fit. In this case, the model is \begin{equation} p_{\rm T.L.}(d_i|\ell,\Delta) = \prod_{i=1}^N p_i \end{equation} where "T.L." means "Triangle Likelihood", and where the likelihood for a single measurement $p_i$ is \begin{equation} p_i = \begin{cases} \frac{4}{\ell \Delta }\left(d_i - \ell\right) , \ \ \ell - \frac{\Delta}{2}< d_i < \ell \\ \frac{4}{\ell \Delta }\left(\ell - d_i\right), \ \ \ell < d_i < \ell + \frac{\Delta}{2} \\ 0, \ \ {\rm otherwise} \end{cases} \end{equation} The parameters in this case are

  • $\ell$ (the true length of the string), and
  • $\Delta$, the width of the resolution, which you could either assume took the value given by the manufacturer (in which case you would not sample over this parameter but just fix it to its given value), or try to measure from the data by taking it as a free parameter.

Alternatively, you could try to make the model more complicated in various ways. One way would be to say that the data consists of the length plus the noise, $d=n+\ell$, and the noise $n$ consists of two contributions that are added together, $n=n_t + n_g$, where $n_t$ is drawn from the triangle distribution and $n_g$ is drawn from a gaussian distribution. Then your likelihood would look something like \begin{eqnarray} p_{\rm M.L.}(d_i|\ell,\sigma,\Delta) &=& \prod_{i=1}^N \int {\rm d} n_t \int {\rm d} n_g \delta(d_i-\ell-n_g-n_t) p_{\rm T.L.}(n_t|\ell,\Delta) p_{\rm G.L.}(n_t|\ell, \sigma) \\ &=& \prod_{i=1}^N \int_{-\Delta/2}^{\Delta/2} {\rm d} n_t p_{\rm T.L.}(n_t|\ell,\Delta) p_{\rm G.L.}(d_i-\ell-n_t|\ell, \sigma) \end{eqnarray} where "M.L." means "Mixture Likelihood." The logic of the first line is that we add all possible values of $\ell, n_t, n_g$, consistent with the idea that the data are a sum of these three quantities $d_i=\ell+n_t+n_g$ (this is the integral + delta function). For each possibility, we multiply the probability $p_{\rm T.L.}$ that the "triangle noise" will be $n_t$ with the probability $p_{\rm G.L.}$ that the "Gaussian noise" will be $n_g$. The second line is one way of evaluating the integral over the delta function, but seems natural to me.

Anyway, the point of this discussion is not so much to suggest any one particular model (I would use the first two as training and the third one as my first "real" case, based on what you said in the question and comments). But, to point out the kind of logic you can use to build models of your uncertainty and fit them into a Bayesian framework.

You also need to decide on priors. Generally, a good starting point is to place uninformative priors on your parameters. While there are lots of debates as to what this means, a simple and justifiable choice is to take uniform priors on every parameter.

  • For $\ell$, you want a wide enough region that the posterior does not "hit the rails" of the range you choose (sometimes this takes some trial and error, running the sampler multiple times). Based on the numbers in your table, a conservative choice would be to take a prior on $\ell$ that ranged from $0$ to twice the "true" value.
  • For $\sigma$, taking a range that went from, say, half the observed standard deviation of the sample means, to twice, should probably be sufficient.
  • For $\Delta$, as stated above, you could fix it to the manufacturer value. Alternatively, you could try to determine this parameter from the data. If you do the latter, I would choose a uniform prior centered on the manufacturer's value, that was not too broad -- unless you start seeing evidence in the posterior that a much different value is strongly preferred.

Once you code up the likelihood and prior in whatever format is needed by your chosen sampler, run the code and look at the results! The results will take the form of samples drawn from the posterior distribution $p(\Theta|d_i)$ of the parameters given the data, obtained using Bayes theorem \begin{equation} p(\Theta|d_i) = \frac{p(d_i|\Theta) p(\Theta)}{p(d_i)} \end{equation} where $p(d_i|\Theta)$ is the likelihood (which we discussed in detail), $p(\Theta)$ is the prior (I gave some advice above), and $p(d_i)$ is a normalization constant sometimes called the "evidence" (which you don't need to compute if you are only interested in the probability of different parameters). At the end of the day, I think you will want to construct the "1-dimensional marginalized posterior" for $\ell$. In other words, assuming we work in the G.L. model, you would want to plot \begin{equation} p(\ell|d_i) = \int_{\sigma_{\rm min}}^{\sigma_{\rm max}} {\rm d} \sigma p(\ell,\sigma|d_i) \end{equation} where $\sigma_{\rm min}$ and $\sigma_{\rm max}$ are the bounds of the prior on $\sigma$. With the output of the MCMC sampler (samples drawn from the posterior distribution), this procedure can be done by simply histogramming the samples by the value of $\ell$.

Some practical advice on interpreting the output of MCMC samplers (and in particular checking for convergence and other problems that can arise) can be found here:

$\endgroup$
6
  • 1
    $\begingroup$ "You need a way to know the length of your calibration string with at least as much precision as your measuring device" The lowest tolerance of sets of engineering "gauge blocks" are 0.25micrometers, (and the highest grade are 0.05 micrometers), so that is not a big deal. $\endgroup$
    – alephzero
    Aug 3 '21 at 1:29
  • $\begingroup$ This is really thorough and well explained, a single upvote doesn't really do it justice! One thing I don't see mentioned is the t-distribution stuff. This may be wrapped up in your analysis but for example the simplest $p_{\rm G.L.}$ is just the product of likelihoods of each data mean in a distribution. If I have a very small number of measurements, then this will be maximised in the MCMC sampler when $\sigma$ is very small because those measurements will initially cluster close to the mean. IIUC the t-distribution fixes this by flattening the distribution based on number of data points DOF $\endgroup$
    – Greedo
    Aug 3 '21 at 7:57
  • $\begingroup$ Also I'm a bit confused - your assumption was correct that all measurements refer to the same $l$. Why then do I need to have it as a free variable in the MCMC sampler? (FWIW, in my actual problem the calibration "length" is 0 by definition, which is why in this problem I'm fixing the calibration to 1000 with infinite resolution) $\endgroup$
    – Greedo
    Aug 3 '21 at 8:08
  • $\begingroup$ @Greedo In the Bayesian approach, you can derive the t-distribution by marginalizing over $\sigma$. So... backing up, if you have a small number of points, what will happen is not that $\sigma$ will be very small, but that the uncertainty of $\sigma$ will be very large (so you don't know what $\sigma$ is). If the data are completely uninformative, then the distribution for $\sigma$ will be flat. If you marginalize (integrate) over $\sigma$, then the distribution for $\ell$ will be a student t distribution. reference $\endgroup$
    – Andrew
    Aug 3 '21 at 13:20
  • $\begingroup$ (Actually strictly speaking the "uninformative" distribution you need for the above argument to produce a student's t distribution is a "log uniform" prior $\sim 1/\sigma^2$, but I'm speaking more conceptually to get across the idea of how the Bayesian formalism builds in the student t distribution). Having said that, a reasonable approach would be to replace $p_{\rm G.L.}$ with a student t distribution if you think a better model of the noise is that it has wider tails than a Gaussian. $\endgroup$
    – Andrew
    Aug 3 '21 at 13:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.