Estimating distribution characteristics from characteristics of multiple samples Definition
Suppose that $X \sim D(\mu, \sigma)$, where $D$ is a 1D distribution (generating from $\mathbb{R}$) with mean $\mu$ and stddev $\sigma$.
If I sample $X$ (the random variable) $M\times N$ times, receiving a matrix $A$ (with $M$ rows and $N$ columns), then generate

*

*a vector $V_\mu$ from $A$ by reducing each row of $A$ to a single value by calculating the mean of that row, and


*a vector $V_\sigma$ from $A$ by reducing each row of $A$ to a single value by calculating the biased standard deviation of that row from the mean --> single number of $V_\sigma = \sqrt{\frac{1}{N}\sum_{i=0}^{N}(x_i - \overline{x})^2}$,
how do I properly estimate the original $\mu$ and $\sigma$ from $V_\mu$ and $V_\sigma$? It's clear that the estimate of $\mu$ is just the mean of $V_\mu$, but how do I handle the $\sigma$?
Sigma
I've tried experimenting, deciding that, for starters, let's assume that $D = \mathcal{N}$ and $X \sim \mathcal{N}(\mu=3, \sigma=5)$. I generated $A$, computed unbiased $\sigma(A)$ and then reduced $A$ to $V_\sigma$, from which I calculated the mean. This is the result:

Obviously, the result is biased, the estimate undershoots the true $\sigma$ more often than not. If I instead multiply each member of $V_\sigma$ by $\frac{N}{N-1}$ before calculating the mean of $V_\sigma$, I get:

which overshoots the true $\sigma$ more often than not, so I'm really at a loss here. I've also tried multiplying $V_\sigma$ by $\frac{MN}{MN-1}$ instead, but this yields

which is slightly better than the original, but still heavily biased.
EDIT: Thanks to Ryan for pointing out my mistake, of course, I forgot to square root the correction factor. Still, I had no idea that the $c_4$ factor should be taken into account as well. By multiplying the mean of $V_\sigma$ by $\sqrt{\frac{N}{N-1}}$ and also by $1/c_4(M)$, I've obtained:

which is unbiased, but distributed more uniformly (with larger standard deviation).
Mu
The mean estimate is good, as expected:

Question
After everything that I've shown, what I wonder is this:

*

*Given $D = \mathcal{N}$, what are these distributions I'm observing? They look normal, but aren't they t? EDIT I know now that the distribution for mean estimate, if we subtract the true $\mu$, is a $t$ distribution with $N-1$ degrees of freedom.


*SOLVED COMPLETELY Given $D = \mathcal{N}$, how do I correct the $\sigma$ estimate from $V_\sigma$? EDIT: answered by Ryan, see Sigma section.


*Given $D = \mathcal{N}$, after I correct the $\sigma$ estimate, is it okay to claim that $D$ is probably $\mathcal{N}(\mu, \sigma)$? Certainly, the larger $M$ and $N$ get, the more confident I can be while claiming such a fact, right? What is the proper statistical procedure I should execute after I get my estimate of $\mu$ and $\sigma$? For example, according to my experiment, I can see that the $\mu$ estimate falls into $[2, 4]$ about 95% of the time. EDIT: I know now that if I generate the interval as $\overline{x} \pm 1.96 \cdot \sigma \cdot \sqrt{N}$, the interval will contain $\mu$ 95% of the time. But what about $\sigma$? And when I finally settle on some interval estimates of $\mu$ and $\sigma$, can something about $P(X > c), X \sim \mathcal{N}(...)$ be said?


*Can this problem be solved for general $D$, that is, obtain $V_\mu$ and $V_\sigma$ from a bunch of samples ($M\times N$ to be specific) and conclude something about the true $\mu$ and $\sigma$?
 A: Welcome to CV!
(1) Sampling Distribution of Sample Statistics
Given $D = \mathcal{N}$, what are these distributions I'm observing? They look normal, but aren't they t?
a) Sample mean $\bar{X}$
$\frac{\bar{x}-\mu}{S/\sqrt{n}}\sim t_{n-1}$ - t distribution with $n-1$ degrees of freedom.
b) Sample variance $S^2$
$\frac{(n-1)}{\sigma^2}S^2 \sim\chi^2_{n-1}$ chi-squared distribution with $n-1$ degrees of freedom (see Sampling Distribution of Sample Variance)
c) Sample standard deviation $S$
$\sqrt{\frac{(n-1)}{\sigma^2}}S \sim\chi_{n-1}$ - chi distribution with $n-1$ degrees of freedom. This follows from the fact that if $X\sim \chi(n)$ then $X^2\sim \chi^2(n)$ (see Wikipedia: Chi Distribution)
(2) Unbiased Estimate of Population $\sigma$
Given $D = \mathcal{N}$, how do I correct the $\sigma$ estimate from $V_\sigma$?
The sample variance with Bessel's correction ($\tfrac{n}{n-1}$) provides an unbiased estimate for the population variance. Two reasons that statement doesn't help you.

*

*You are applying Bessel's correction $\frac{n}{n-1}$ to the sample standard deviation. In fact, you would want to multiply the sample standard deviation by $\sqrt{\frac{n}{n-1}}$ to apply the correction.

*Even then, you will not get an unbiased estimate of the sample standard deviation. The corrected variance is unbiased, but the square root of that value is not an unbiased estimate of population standard deviation. See Wikipedia and related question. In the case where $D = \mathcal{N}$, there is a correction factor ($c_4(n)$) you can apply. It is discussed in the wikipedia article linked above. For the case where $n=10$, the correction looks like $c_4(10)= \left(\frac{128}{105}\sqrt{\frac{2}{\pi}}\right)\approx 0.9726592741$.

In general, an unbiased estimate for the population standard deviation where $D = \mathcal{N}$ is given by
$$\hat{\sigma}=\frac{1}{c_4(n)}\sqrt{\frac{\sum_{i=1}^N(x_i-\bar{x})^2}{N-1}}$$
Here is a quick plot to show the difference in the estimated standard deviation using the two corrections on samples from a normal distribution as well as Python code to reproduce the plot.

from math import gamma
import seaborn as sns
import pandas as pd

SIGMA = 5
MU = 3
m = 10000
# calculate correction
def c4(n):
    return np.sqrt(2/(n-1)) * gamma(n/2) / gamma((n-1)/2)

# calculate statistics for various N
results_dict = {x:[] for x in ['N','correction','s']}
for N in range(3, 25):
    A = np.random.normal(loc=MU, scale=SIGMA, size=[m,N])
    df_i = pd.DataFrame()
    results_dict['N'] += [N]*m*3
    results_dict['correction'] += ['None']*m+['Bessel']*m+['Bessel + c4']*m
    results_dict['s'] += list(np.std(A, axis=1))
    results_dict['s'] += list(np.std(A, axis=1)* ((N/(N-1))**0.5) )
    results_dict['s'] += list(np.std(A, axis=1)* ((N/(N-1))**0.5) / c4(N))

# create dataframe
results_df = pd.DataFrame(results_dict)

# plot results
plt.figure(figsize=(8,6))
sns.pointplot(
    data=results_df,
    x='N',
    y='s',
    hue='correction',
    ci=None
)
plt.title("Comparison of statistics for estimating $\sigma$")
plt.axhline(5, c='k', linestyle='--', label= "$\sigma$")
plt.show()

(3) Confidence Intervals - Normal
Given $D = \mathcal{N}$, after I correct the $\sigma$ estimate, is it okay to claim that $D$ is probably $\mathcal{N}(\mu, \sigma)$? Certainly, the larger $M$ and $N$ get, the more confident I can be while claiming such a fact, right? What is the proper statistical procedure I should execute after I get my estimate of $\mu$ and $\sigma$? For example, according to my experiment, I can see that the $\mu$ estimate falls into $[2, 4]$ about 95% of the time.
You of course can't say that $D$ is probably exactly $\mathcal{N}(\bar{X}, S)$ but you can construct confidence intervals for $\mu$ and $\sigma$.
As an aside, the maximum likelihood estimator for the variance is actually the uncorrected version $s^2=\frac{1}{N}\sum_{i=1}^N(X_i-\bar{X})^2$ (see MLE Biased). This is regardless of the fact that the uncorrected estimate tends to underestimate the true value. And if $S^2$ is the MLE estimate for $\sigma^2$ then $\sqrt{S^2}=S$ is the MLE estimate for $\sqrt{\sigma^2}=\sigma$ (see Maximum Likelihood Estimation) We can also see, using our simulation, that the average squared difference between our estimate $S^2$ and the population variance $\sigma^2$ is lowest for the uncorrected estimate.

# Variance
results_df['s2'] = results_df['s']**2
# Variance error
results_df['s2_mse'] = (results_df['s2']-SIGMA**2)**2

plt.figure(figsize=(8,6))
sns.pointplot(
    data=results_df,
    x='N',
    y='s2_mse',
    hue='correction',
    ci=None
)
plt.ylabel("$(S^2-\sigma^2)^2$")
plt.title("Squared Error of statistics for estimating $\sigma^2$")
plt.axhline(5, c='k', linestyle='--', label= "$\sigma$")
plt.show()

You can construct the following confidence intervals for your sample statistics.
a) Population mean $\bar{X}$
A $(1-\alpha)%$ confidence interval for the population mean is
$$\left(
\bar{X}-\frac{S}{\sqrt{n}}t_{n-1,\alpha/2}
\leq 
\mu
\leq 
\bar{X}+\frac{S}{\sqrt{n}}t_{n-1,\alpha/2}
\right)$$
see: Confidence Intervals with σ unknown
b) Population variance $\sigma^2$
A $(1-\alpha)%$ confidence interval for the population variance is
$$\left(\frac{(n-1)s^2}{\chi^2_{\alpha/2,n-1}} 
\leq 
\sigma^2 
\leq 
\frac{(n-1)s^2}{\chi^2_{1-\alpha/2,n-1}}
\right)$$
(see Confidence Intervals for Variances)
c) Population standard deviation $\sigma$
A $(1-\alpha)%$ confidence interval for the population standard deviation is
$$\left(\sqrt{\frac{(n-1)s^2}{\chi^2_{\alpha/2,n-1}}}
\leq 
\sigma^2 
\leq 
\sqrt{\frac{(n-1)s^2}{\chi^2_{1-\alpha/2,n-1}}}
\right)$$
See Confidence Intervals for Variances again or this related question
(4) General Case
Can this problem be solved for general $D$, that is, obtain $V_\mu$ and $V_\sigma$ from a bunch of samples ($M\times N$ to be specific) and conclude something about the true $\mu$ and $\sigma$?
@Ben's answer seems most relevant to this question. A similar procedure you may be able to apply in the general case and avoid any analytical computation is to utilize the bootstrap. You can certainly use the bootstrap procedure to estimate the distribution of the sample mean, but I am not sure how well it applies to estimating sample variance. I'm not finding a lot of specific discussion on this question, but this thesis seems to discuss the issue in depth. Evaluation of Using the Bootstrap Procedure to Estimate the Population Variance 
