I think the easiest way to to what you want is to use confidence intervals (statistical inference).
In other words, assuming the population has a true variance $\sigma$, the sampling distribution of the variance $s^2$ of an $n$-sample verifies:
$$ \frac{s^2(n-1)}{\sigma^2}\sim \chi^2_{n-1}$$
You can exploit this result to build an $1-\alpha$ confidence interval for the population variance ($\alpha \in [0,1]$, typically $\alpha=5\%$).
Indeed, for a confidence level $1-\alpha$, the following equality holds:
$$ z_{\alpha/2} \leq \frac{s^2(n-1)}{\sigma^2} \leq z_{1-\alpha/2} $$
where $z_q$ figures the quantile $q$ of a chi-squared distribution with $n-1$ degrees of freedom i.e.
$$ X \sim \chi^2_{n-1},\quad \Bbb{P}(X \leq z_q) = q $$
Given a sample variance $\tilde{s}^2$, one can therefore turn the inequality on its head, to write, for a confiedence level $1-\alpha$:
$$ \frac{\tilde{s}^2(n-1)}{z_{1-\alpha/2}} \leq \sigma \leq \frac{ \tilde{s}^2(n-1)}{z_{\alpha/2}} $$
Hence the upper and lower bounds of your $1-\alpha$ confidence interval for the (unobserved) population variance:
\begin{align}
\sigma^+ = \frac{ \tilde{s}^2(n-1)}{z_{\alpha/2}},\quad \sigma^- = \frac{\tilde{s}^2(n-1)}{z_{1-\alpha/2}}
\end{align}
This could then help you construct $1-\alpha$ confidence bounds on the BS option price given the measure sample variance $\tilde{s}^2$:
$$ V^+ = \text{BSCall}(\sigma^+),\quad V^- = \text{BSCall}(\sigma^+)$$
[Edit]
Given your desire to obtain a full distribution, why not opt for a Bayesian approach?
Assume the true population variance $\sigma^2$ follows a certain prior distribution with hyperparameter $\alpha$, $p(\sigma;\alpha)$ over $\Bbb{R}^+$.
Suppose that, for a specific sample, you measure a sample variance $s^2$ and wish to compute the posterior of the population variance. Bayes' rule gives:
$$ p(\sigma^2 \vert s^2, \alpha) = \frac{p(s^2 \vert \sigma^2)}{\int_0^\infty p(s^2 \vert \sigma^2) p(\sigma^2;\alpha) d\sigma^2 } p(\sigma^2; \alpha) $$
Now you know:
- The prior distribution $p(\sigma^2; \alpha)$: you postulated it.
- The sampling distribution $p(s^2 \vert \sigma^2)$: $\quad s^2 \sim \sigma^2/(n-1) \chi^2_{n-1}$
Hence you have everything you need to compute the posterior distribution.
Obviously, if you stick with the Maxium A Posteriori (MAP) estimator, once again you'll have a pointwise estimate, so I suggest you to perform the full integration. Off the top of my head chi-squared distributions does not allow for conjugate priors so you might have to resort to numerical integration (e.g. adaptive quadrature and the likes).
Finally, the choice of hyper-parameter $\alpha$ will have an impact on the resulting posterior: you might want to set $\alpha$ so that the prior distribution is centered around the sample variance for instance?