I understand that the point is how much the variability of the $\tau$ parameter will affect the shape of the $M(t)=M_0\cdot E^{-t/\tau}$ function.
$\tau$ obay a certain probability distribution.
We nead to obtain $\tau$ values for the probabilities $p$ that set the boundaries of ever wider confidence intervals. This is done using InverseCDF[]
. I have selected the fallowing probabilities: {0.125,0.25,0.375,0.5,0.625,0.75,0.875}
.
For respectve $\tau$ values, we define the function $M(t)$.
Now, instead of a single $M(t)$ function, we plot a family of functions that differ in the $\tau$ parameter, so that the 'central'* function (the one with $p=0.5$) is marked in black and its less and less likely variants are gray.
Such a plots ware crated for the $\tau$ obaying the LogNormalDistribution[μ,σ]
** - where different values of $\mu$ and $\sigma$ were tried.
pp = Table[i, {i, 0.125, 1 - 0.125, 0.125}]
\[Mu]\[Mu] = Table[i, {i, 1, 5, 1}];
\[Sigma]\[Sigma] = {0.5, 1, 2};
dist1 := LogNormalDistribution[\[Mu], \[Sigma]]
\[Tau][p_] := InverseCDF[dist1, p]
M[t_, p_] := E^(-t/\[Tau][p])
MM := Table[M[t, p], {p, pp}]
rowh = Table[StringJoin["\[Mu]=", ToString[i]], {i, \[Mu]\[Mu]}];
colh = Table[
StringJoin["\[Sigma]=", ToString[i]], {i, \[Sigma]\[Sigma]}];
TableForm[
Table[Plot[MM, {t, 0, 25},
PlotStyle -> {Gray, Gray, Gray, Black, Gray, Gray, Gray},
Filling -> {1 -> {7}, 2 -> {6}, 3 -> {5}},
FillingStyle -> {Gray, Opacity[0.5]},
PlotRange -> {0,
1}], {\[Mu], \[Mu]\[Mu]}, {\[Sigma], \[Sigma]\[Sigma]}],
TableHeadings -> {rowh, colh}]

Similarly for NormalDistribution[μ,σ]
, but here some of $\tau$'s fell below $0$ causing exponential explosion of $M$.

*) Please note that this one corresponds to the median of the $\tau$ distribution and not to mode (the most probable value).
**) Please pay attention to the meaning of parameters in log-normal distribution - they are non-trivial!