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:
