Should Bayesian probability of model given the data coincide with intuitive estimation of this probability? I have $N$ real valued targets ordered by one single real valued feature. I would like to build a simple probabilistic model that states that all the targets corresponding to the features that are smaller than $f_{split}$ are coming from one normal distribution given by mean $\mu_1$ and dispersion $\sigma_1$ and all the other targets coming from another normal distributions with another mean $\mu_2$ and another dispersion $\sigma_2$.
So, my model is given by 5 parameters (split value, two means and two dispersion). When I search the best model by maximization of likelihood, the best model splits the data such that in one basket of targets I have just 1 point and the corresponding distribution has mean equal to the value of that single target ($\mu_1 = t_1$) and dispersion is equal to zero ($\sigma_1 = 0$).
I do understand why it is happening. The single target gets this way a huge (practically infinite) probability density and therefore the total likelihood is very large.
My question is: What is the best way to avoid this problem. From an intuitive point of view, if I see by eyes that there is a split (targets on the left hand side and on the right hand side clearly have different averages and dispersion), I would not believe that the real model is the one which contains just one target on the left. So, why my intuitive estimation of model given data do not coincide with the Bayesian one (I assume that all the $a\ priori$ models have the same probability).
I guess my intuition says that it is very unlikely that just by chance we had a model whose mean is equal to the target and whose $\sigma$ is so close to zero. But I do not know how to formalize it.

The intuitive notion that you are trying to formalize is exactly what a proper Bayesian analysis will provide. To see this, consider for a example the following simple model, which has similar features to the case you are considering:
$$ X \sim 0.9\mathcal N(\mu,\sigma^2) + 0.1 \mathcal N(0,1)$$
here $X$ is assumed to be generated from a mixture of a normal distribution with unknown parameters $\mu,\sigma$, and a small "contamination" having a known standard normal distribution.
Given a set of observations $x_i$, the likelihood is (up to irrelevant constants) :
$$\log \mathcal L (\mu,\sigma) = \sum_i \log \left( 0.1e^{-\frac{1}{2}x_i^2} + \frac{0.9}{\sigma}e^{-\frac{1}{2}(x_i-\mu)^2/\sigma^2} \right)$$
and it diverges to infinity at the points $(\mu=x_i,\sigma=0)$. However, despite the fact the likelihood is not defined at those points, the posterior probability of $\mu$ and $\sigma$
$$ P(\mu,\theta | x ) \propto \mathcal L (\mu,\sigma)\pi(\mu,\theta) $$
with any proper prior $\pi(\mu,\theta)$, is a perfectly valid probability distribution.
For example, this is a plot of a slice of  $\log \mathcal L (\mu,\sigma)$ along  the line $\mu=x_1$ when the data is $x=[1.0309,0.4307,-0.7551,0.3038,0.4514]$ (generated with a "true" $\sigma=1$ )

With any prior that is more or less uniform in the range of the plot, the posterior expectation of $\theta$ would be around 1, since the divergent spike at $\sigma=0$ contains only a negligible fraction of the posterior probability. This corresponds to your intuition that "it is very unlikely that just by chance we had a model whose mean is equal to the target and whose σ is so close to zero".
This is exactly the meaning of the infinity at $\sigma=0$ : if someone had a-priori a model that predicts exactly the value of $x_1$ we observed, namely if the prior distribution had a non-zero point mass at $(\mu=x_i,\sigma=0)$, then the posterior would be a $\delta$-function at this point: the probability of observing exactly this predicted value by chance is zero, so wo would have to conclude the this model is correct with absolute certainty. But if we consider a-priori all values of $\sigma$ and $\mu$ to be "equally" possible, then we are assigning zero prior probability to the exact point $\sigma=0$.
