How to force splits in decision tree to be distributed uniformelly in case of no dependency on feature? I have targets ordered by a feature. I want to find a single split that minimizes a squared deviation (RMSE). For example, I have 100 values (targets) and it might be the case that, if I split them as 30 + 70, I get minimal squared deviation (as compared with all the other splits like 10+90 or 83+17). Here I assume that 30 targets "on the left" get one prediction (average of those 30 targets) and 70 targets on the right get another prediction (average over 70 values).
To test my split function, I have generated 100 values using a single normal distribution and for those values I found an optimal split given by its location (for example location 15 means 15 values on the left and 85 values on the right). When I do it many times I find that split locations are distributed non uniformly. For example very small locations (3, 4, 5) or very large locations (97, 98) are much more frequent than locations close to the middle (48, 50, 52). In other words, good RMSEs are more likely to be found for the splits close to the two edges of the whole interval.
So, it looks to me as an inherent undesirable bias of the method. In the case of no pattern in the data (no dependency on the feature, as in my test), I would expect splits to be distributed uniformly between 1 and 100. Otherwise, when I do have a dependency on the feature, mixed with large noise, I still would get close to the edges of the intervals and I do not want it.
By the way, my explanation of the effect is that root means squares deviations corresponding to the splits close to the edges are less correlated with each other and therefore the chances that we find minimal root mean square there is larger.

One way to approach dataset splitting is as a model selection problem. Consider the following set of hypotheses :
$$ H_k : y_i \sim 
\begin{cases}
    \mathcal N(\mu_1,\sigma^2),& \text{if } i < k\\
    \mathcal N(\mu_2,\sigma^2),              & \text{if } i \ge k
\end{cases}$$
Namely, the set of data points $y_i$ are assumed to belong to two populations with different (unknown) means. Now calculate the posterior probabilities of each hypothesis :
$$ P(y|H_k) = \int d \mu_1 P(y^{(1,k)} | \mu_1,\sigma) \pi(\mu_1) \times \int d \mu_2 P(y^{(2,k)} | \mu_2,\sigma) \pi(\mu_2)$$
where $y^{(1,k)} = \{y_i | i<k\}$ and $y^{(2,k)} = \{y_i | k \le i \le N\}$.
Assuming uniform priors for the unknown means, we have
$$\int d \mu_1 P(y^{(1,k)} | \mu_1,\sigma) =  \frac{1}{(2\pi)^{k/2}\sigma^k}\int d \mu_1 \prod_{i=1}^k e^{-\frac{1}{2\sigma^2}(y_i-\mu_1)^2} \\ = \frac{1}{(2\pi)^{k/2}\sigma^k}\int d \mu_1 \exp\left({-\frac{1}{2\sigma^2} \sum_{i=1}^k(y_i-\mu_1)^2} \right) \\ = \frac{1}{(2\pi)^{k/2}\sigma^k}\int d \mu_1 \exp\left( -\frac{k}{2\sigma^2} (\mu_1 - \bar y^{(1,k)})^2 \right) \times \exp \left( -\frac{1}{2\sigma^2} \text{SSE}^{(1,k)}\right) \\ = \frac{\sqrt{2\pi}\sqrt{k}}{(2\pi)^{k/2}\sigma^{k-1} }\exp \left( -\frac{1}{2\sigma^2} \text{SSE}^{(1,k)}\right)$$
where
$$ \text{SSE}^{(1,k)} = \sum_{i=1}^k (y_i - \bar y^{(1,k)})^2 $$
is the sum of squared errors. Combined with a similar expression from the second integral we get
$$ P(y|H_k) \propto \sqrt{k(N-k)} \exp \left( -\frac{1}{2\sigma^2} \text{SSE}^{(k)}\right)$$
where $\text{SSE}^{(k)} = \text{SSE}^{(1,k)} + \text{SSE}^{(2,k)}$ is the total sum of squared errors. If we assign the same prior probability for each hypothesis then $P(H_k|y) \propto P(y|H_k)$, and if we take the $\log$ of the posterior probabilities we can see that
$$ -\sigma^2\log P(H_k | y ) = \frac{1}{2}\text{SSE}^{(k)} - \frac{\sigma^2}{2} \log k(N-k) + const.$$
So, the splitting which has the maximum a posteriori probability is the one that minimizes  $\text{SSE}^{(k)} - \sigma^2\log k(N-k)$ . Note that the second term effectively penalizes small leaves.
Note that the model so far assumed that the variance $\sigma^2$ is known and common to both leaves. If  $\sigma^2$ is not known we can use the sample estimate $\hat \sigma^2 = \frac{1}{N}SSE^{(k)}$, so the minimization objective becomes $\text{SSE}^{(k)}(1  - \frac{1}{N}\log k(N-k))$. (In principle if the variances are not known we should marginalize over them, similar to way we treated the unknown means, which might lead to slightly different results) .
To see the effect of choosing a splitting based on the posterior probabilities we can run a test similar to yours (here with $\sigma^2=1$ for simplicity) :

for n=1:10000
  y = randn(100,1) ;
  for k=1:99
    SSE(k) = k*var(y(1:k),1) + (100-k)*var(y(k+1:end),1) ;
    logP(k) = SSE(k) - log(k*(100-k)) ;
  end
[~,loc1(n)] = min(SSE) ; %location that minimizes SSE
[~,loc2(n)] = min(logP) ; %location that minimizes posterior probability
end
figure, histogram(loc1,50)

Splitting that minimizes $SSE$ :

figure, histogram(loc2,50)

Splitting that maximizes $P(H_k|y)$ :

This is indeed distributed uniformly as we would expect.
