I'm working with a set of independent and identically distributed random variables $\{ x_i \}_{i=1}^N$, where each $x_i$ follows a Gaussian distribution $P_X(x) = \mathcal{N}(x; \mu, \sigma^2)$. This notation $\mathcal{N}(x; \mu, \sigma^2)$ represents the Gaussian probability density function (p.d.f.) with mean $\mu$ and variance $\sigma^2$. It's straightforward to show that, by defining $y_i \equiv \frac{x_i - \mu}{\sigma}$, we obtain $y_i \sim P_Y(y) = \mathcal{N}(y; 0, 1)$.

I'm interested in exploring a similar transformation but using empirical estimates of the mean and standard deviation instead:
$$z_i \equiv \frac{x_i - \hat{\mu}}{\hat{\sigma}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(1)$$
 where $\hat{\mu}$ and $\hat{\sigma}$ are the empirical mean and standard deviation, respectively, given by:
$$\hat{\mu} \equiv \frac{1}{N} \sum_{i=1}^N x_i\; \;\;\;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\;\;\;\;\;\;(2)$$
$$\hat{\sigma^2} \equiv \frac{1}{N} \sum_{i=1}^N (x_i - \hat{\mu})^2 \;\;\;\;\;\;\;\;\;\;\;\;\;\; (3)$$

**My question is: How can we derive the distribution of $z_i$ for a given fixed $N$?**.

It's important to note the differences from the initial example:

- Unlike $\mu$ and $\sigma$, which are constants, $\hat{\mu}$ and $\hat{\sigma}$ are random variables dependent on the set ${ x_i }_{i=1}^N$.
- Both $\hat{\mu}$ and $\hat{\sigma}$ are functions of the $x_i$'s.
- The distribution of $z_i$ is generally not normal. For instance, with $N=2$, any two distinct $x_1, x_2$ will be mapped to $z_i \in \{ -1, +1 \}$, resulting in a distribution for $z_i$ of $P_Z(z) = \frac{1}{2} \delta(z-1) + \frac{1}{2} \delta(z+1)$.