Bayes Questions

I'll walk through how I'd analyze this problem, let me know if I haven't answered your questions by the end.

First, problem structure. You have three unknowns, which you want to estimate from the data: a shape parameter, a scale parameter, and an indicator of which model is correct.

"Which model is correct" is probably easiest to start with. I'm not completely sure that I've followed what your spreadsheet is doing, but if I understand correctly, then that's probably an overly-complicated way to tackle the problem. You want $P\left(model|data\right)$, which will be determined via Bayes' Rule by $P\left(data|model\right)$ and your prior probability for each model being correct. The prior is unlikely to matter much unless your dataset is tiny, so $P\left(data|model\right)$ is the important part. That's an integral:

$P\left(data|model\right)={\int }_{\mu ,\sigma }{\prod }_{i}P(x_i|\mu ,\sigma ,model)dP(\mu ,\sigma )$

In your case, you're approximating that integral with a grid over $\mu$ and $\sigma$ ($dP(\mu ,\sigma )$ is a shorthand for $\rho (\mu ,\sigma )d\mu d\sigma$ here). Rather than whatever you're doing with timesteps, you can probably just take the product of $P(x_i|\mu ,\sigma ,model)$, where $x_i$ is the lifetime of the $i^{th}$ component in your dataset, then sum over the grid. (If you are dealing with online data streaming in, then you would need to do the timestep thing.)

That takes care of the model part. Once that's done, you'll probably find that one model is like a gazillion times more likely than the other two, and you can just throw out the other two.

On to the 95% CI for the worst part in a million.

The distribution you're interested in here is $P(x_{worst}|model,data)$. $x_{worst}$ is an order statistic. Its CDF is basically the CDF for any old point $x$ raised to the power of 1000000; read up on it to see exactly what expression to use. So if we wanted to do this analytically, we'd first compute $P(x|model,data)$ via Bayes' Rule:

$P(x|model,data)=\frac{P(data\cup x|model)}{P\left(data|model\right)}$

... where both pieces on the right would involve our integral from earlier. Basically, you imagine adding one more point $x$ to the dataset and see what that would do to $P\left(data|model\right)$. If we had a closed-form expression for that distribution, then we could just raise the CDF to the millionth power, we'd get a closed-form expression for the millionth order statistic, and from there we'd get a 95% CI in the usual way.

In practice, that's probably difficult, so let's talk about how to approximate it numerically.

First, the order statistic part. As long as we can sample from the posterior distribution $P(x|model,data)$, that part's easy: generate a million samples of $x$, take the worst, and you have a sample of $x_{worst}$. Repeat that process a bunch of times to compute the 95% CI. (This is not the same as the worst component in 20M, but it's not any harder to code up.)

Next, the posterior distribution for $x$. This is going to be driven by two pieces: uncertainty in the parameters $\mu$ and $\sigma$, and random noise from the distribution itself. If the dataset is large enough, then the uncertainty in $\mu$ and $\sigma$ will be small, so the distribution itself will be the dominant term. In that case, we can just find the best-fit (i.e. maximum a-posteriori) estimates of $\mu$ and $\sigma$, and then declare that $P(x|model,data)$ is approximately the standard distribution (Weibull, log-normal, etc) with those exact parameters. Presumably we can sample from any of those distributions with known parameters, so we go do the order statistic part and we're done.

If the uncertainty in $\mu$ and $\sigma$ is not small enough to ignore, then the problem gets more complicated - we'll need to sample from the posterior distribution $P(\mu ,\sigma |model,data)$. At that point we're in the land of Laplace approximation and MCMC and all that jazz; I'm not going to walk through it here, because this comment is already really long.

So one last thing to wrap it up. I wrote all that out because it's a great example problem of how to Bayes, but there's still a big problem at the model level: the lifetime of the millionth-worst component is probably driven by qualitatively different processes than the vast majority of other components. If some weird thing happens one time in 10000, and causes problems in the components, then a best-fit model of the whole dataset probably won't pick it up at all. Nice-looking distributions like Weibull or log-normal just aren't good at modelling two qualitatively different behaviors mixed into the same dataset. There's probably some standard formal way of dealing with this kind of thing - I hear that "Rare Event Modelling" is a thing, although I know nothing about it - but the fundamental problem is just getting any relevant data at all. If we only have a hundred thousand data points, and we think that millionth-worst is driven by qualitatively different processes, then we have zero data on the millionth-worst, full stop. On the other hand, if we have a billion data points, then we can just throw out all but the worst few thousand and analyse only those.

Can you program? In that case I highly recommend using PyMC or Stan for this kind of work. There's a pretty rich literature and culture of how to iteratively improve these types of models, and some tool support around these specific toolkits.