The question is about applying a transformation to a given sample to obtain the corresponding sample with continuous standard uniform distribution.

I'll start by describing a natural and convenient approach to perform this type of transformation, and provide an example of the method with data considered by the OP using R. Then, I'll briefly discuss the approach considered by the OP using a quantile transformation based on subsamples.

Non-parametric transformation

A natural way to apply the probability integral transform without specifying a parametric distribution is to use the empirical distribution function (ECDF).

To avoid boundary problems, it is more convenient to use the modified empirical distribution function $$ F_n(x) = \frac{1}{n+1}\sum_{i = 1}^n I_{\{x_i\leq x\}} , $$ where $I(\cdot)$ denotes the indicator function. Applying the probability integral transform to the dataset $X_1, \ldots, X_n$ using $F_n$ corresponds to compute the normalized ranks $$ U_i = \frac{R_i}{n+1} , $$ where $R_i$ denotes the rank of $X_i$ in increasing order (i.e., $R_k = 1$ if $X_k$ is the smallest observation).

Example

The following figure shows a histogram of the simulated dataset of $10^4$ values from the standard normal distribution (left), a histogram of the transformed data (middle), and the transformed values $U_i$ versus the original values $X_i$ (right). By construction the distribution of the transformed variables is uniform between 0 and 1.

The R code used to run this simulation is shown below.

## Parameters
n <- 10000
## Data set generation
set.seed(345)
x <- rnorm(n)
## Transformation
u <- rank(x) / (n + 1)
## Result visualization
par(mfrow = c(1, 3))
hist(x, breaks = 101, main = "")
hist(u, breaks = 101, main = "")
plot(x, u)
par(mfrow = c(1, 1))

Comments

Alternative estimates of the underlying cumulative distribution function could also be used (see here for a discussion).

The quantile-based transformation with interpolation considered by the OP is an option which corresponds to a "smoothed" version of the ECDF. However, this approach discards information contained in the dataset since it uses $m<n$ observations, resulting in a less efficient estimator of the underlying distribution. Increasing $m$ improves the quality of the estimator, so chosing $m=n$ seems most appropriate.

The question is about applying a transformation to a given sample to obtain the corresponding sample with continuous standard uniform distribution.

I'll start by describing a natural and convenient approach to perform this type of transformation, and provide an example of the method with data considered by the OP using R. Then, I'll briefly discuss the approach considered by the OP using a quantile transformation based on subsamples.

Non-parametric transformation

A natural way to apply the probability integral transform without specifying a parametric distribution is to use the empirical distribution function (ECDF).

To avoid boundary problems, it is more convenient to use the modified empirical distribution function $$ F_n(x) = \frac{1}{n+1}\sum_{i = 1}^n I_{\{x_i\leq x\}} , $$ where $I(\cdot)$ denotes the indicator function. Applying the probability integral transform to the dataset $X_1, \ldots, X_n$ using $F_n$ corresponds to compute the normalized ranks $$ U_i = \frac{R_i}{n+1} , $$ where $R_i$ denotes the rank of $X_i$ in increasing order (i.e., $R_k = 1$ if $X_k$ is the smallest observation).

Example

The following figure shows a histogram of the simulated dataset of $10^4$ values from the standard normal distribution (left), a histogram of the transformed data (middle), and the transformed values $U_i$ versus the original values $X_i$ (right). By construction the distribution of the transformed variables is uniform between 0 and 1.

The R code used to run this simulation is shown below.

## Parameters
n <- 10000
## Data set generation
set.seed(345)
x <- rnorm(n)
## Transformation
u <- rank(x) / (n + 1)
## Result visualization
par(mfrow = c(1, 3))
hist(x, breaks = 101, main = "")
hist(u, breaks = 101, main = "")
plot(x, u)
par(mfrow = c(1, 1))

Comments

Alternative estimates of the underlying cumulative distribution function could also be used (see here for a discussion).

The quantile-based transformation with interpolation considered by the OP is an option which corresponds to a "smoothed" version of the ECDF. However, this approach discards information contained in the dataset since it uses $m<n$ observations, resulting in a less efficient estimator of the underlying distribution. Increasing $m$ improves the quality of the estimator, so chosing $m=n$ seems most appropriate.