The joint sampling distribution of $\bar{x}$ and $S$ is derived in integral form for probability density functions of doubly infinite range. This derivation is effected through the use of a transformation which transforms the sample probability element $f(x_1)f(x_2) \cdots f(x_n) dx_1 dx_2 \cdots dx_n$ into the element $f(x_1)f(x_2) \cdots f(x_{n - 2})f((n\bar{x} - \sum^{n - 2}_1 x_i \pm \Omega_1)/2)f((n\bar{x} - \sum^{n - 2}_1 x_i \mp \Omega_1)/2) |J| \\ \cdot dx_1 dx_2 \cdots dx_{n - 2}d\bar{x} dS,$ where $\bar{x} = (1/n) \sum^n_1 x_i, S^2 = (1/n) \sum^n_1(x_i - \bar{x})^2$, and $J$ is the Jacobian of the transformation. Bounds on $x_{n - r}, r = 2, 3, \cdots, n - 1$, are established in terms of $\bar{x}, S$, and $x_{n - r -j}, j = 1, 2, \cdots, n - r - 1$. The probability element $f(x)_1f(x_2) \cdots f(x_{n - 2})f((n\bar{x} - \sum^{n - 2}_1 x_i \pm \Omega_1)/2)f((n\bar{x}\\ - \sum^{n - 2}_1 x_i \mp \Omega_1)/2)|J|dx_1 dx_2 \cdots dx_{n - 2}d\bar{x} dS$ must then be integrated with respect to $x_{n - r}, r = 2, 3, \cdots, n - 1$, between these limits to obtain $F(\bar{x},S)d\bar{x}dS$, the joint probability element of $\bar{x}$ and $S$. These limits of integration of $x_{n - r}, r = 2, 3, \cdots, n - 1$ enable one to express $F(\bar{x}, S)$ in terms of quadratures when $f(x)$ is any probability density function of doubly infinite range. To illustrate the method, $F(\bar{x}, S)$ is obtained when $f(x)$ is the normal probability density function.