Should I standardize data using sample statistics? Suppose I sample 1,000 $\mathcal{N}(0,1)$ observations
>>> import numpy as np
>>> x = np.random.randn(1000)

The sample mean of x is
>>> x.mean()
-0.034215203795331606

The (biased) sample standard deviation of x is
>>> x.std(ddof = 1)
0.9923790554909595

Now suppose that I standardize these observations using these sample statistics
>>> y = (x - x.mean()) / x.std(ddof = 1)

The sample mean and (biased) sample standard deviation of y are
>>> y.mean()
1.5987211554602253e-17
>>> y.std(ddof = 1)
0.9999999999999999

However, I feel that this standardization procedure is a bit misleading, since y is no longer normally distributed. More formally, if
\begin{align}
X &\sim \mathcal{N}(0,1) \\
\hat{\mu} &= \frac{1}{1000} \sum_{i=1}^{1000} X_i \\
\hat{\sigma} &= \sqrt{\frac{1}{999} \sum_{i=1}^{1000} (X_i - \hat{\mu})^2}
\end{align}
then
$$
Y = \frac{X - \hat{\mu}}{\hat{\sigma}}
$$
My questions are:

*

*If $Y$ is not normally distributed, then when is standardizing data using sample statistics appropriate?

*How is $Y$ distributed?

*Does this standardization procedure using sample statistics have a name?

Update
I was not precise enough in my original question. Suppose that
$$
X_1,X_2,...,X_{1000} \sim \mathcal{N}(0,1)
$$
are independent and identically distributed, and suppose that
\begin{align}
\hat{\mu} &= \frac{1}{1000} \sum_{i=1}^{1000} X_i \\
\hat{\sigma} &= \sqrt{\frac{1}{999} \sum_{i=1}^{1000} (X_i - \hat{\mu})^2}
\end{align}
If
$$
\forall i \in \{1,2,...,1000\}, Y_i = \frac{X_i - \hat{\mu}}{\hat{\sigma}}
$$
then are
$$
Y_1,Y_2,...,Y_{1000} \sim \mathcal{N}(\mu_{Y},\sigma_{Y})
$$
for some mean $\mu_{Y}$ and standard deviation $\sigma_{Y}$? It seems that the answer is no, since each $Y_i$ is a non-linear function of $X_1,X_2,...,X_{1000}$ via $\hat{\sigma}$. If I am wrong, then what would $\mu_{Y}$ and $\sigma_{Y}$ be?

If $X$ is a standard normal, than $Y = (X - \mu)/\sigma$ is normally distributed with mean $\mu$ and standard deviation $\sigma$, as normal distribution is of location-scale family.

However, I feel that this standardization procedure is a bit
misleading, since y is no longer normally distributed.

What is y in here? You took 1000 samples from a normal distribution and calculated the sample mean and sample standard deviation of the samples. In statistics, you treat those 1000 observed values as 1000 i.i.d. normally distributed random variables (not as a random variable). You transformed the samples. The statistics you used are estimates of the true parameters, so you can think of the result of the operation as approximation of what you would observe if you knew the population parameters and transformed the random variables.
The calculated statistics are numbers, you used those numbers to do the transformation. You can use any numbers to scale a variable. Subtracting sample mean and dividing by sample standard deviation is called normalization, standardization, or converting to $z$-scores. You use the sample statistics here, not the population parameters, because you don't know the population parameters.
Moreover, you usually care about transforming the samples, not the underlying random variables. Standardization is done mostly for numerical reasons, to ease up the computations. Having samples with a non-zero mean and standard deviation other than one, even if the samples were taken from a standard distribution, would not help with the computations. Don't confuse samples and results of the computations with random variables.
