I want to average over all positive functions on the unit interval whose Lebesgue integral is one. Formally, I want to compute the mean of the following probability distribution defined over function space $L^1[0,1]$ \begin{align} P(f) = \begin{cases} \frac{1}{c} &\text{ when } \int_0^1 dx \; f(x) = 1\; ; f(x) \geq 0\\ 0 &\text{ otherwise} \end{cases} \end{align} where $c$ is a normalization factor.
Noticing the symmetry of the problem regarding variable $x$, we conclude that the result should be a constant. Then using the normalization, we find the mean \begin{equation} \bar{f}(x) = 1 \end{equation} I perform a Monte Carlo simulation on uniform grid in the interval $[0,1]$ and the result confirms the above conclusion...so far so good!
Now I do a variable transformation $ z = x^2 \Rightarrow \frac{dx}{dz} = \frac{1}{2\sqrt{z}} $. The probability distribution then becomes \begin{align} P(f) = \begin{cases} \frac{1}{c} &\text{ when } \int_0^1 dz \; \frac{f(\sqrt(z))}{2\sqrt{z}} = 1\; ; f(z) \geq 0\\ 0 &\text{ otherwise} \end{cases} \end{align} Denoting $g(z)=\frac{f(\sqrt(z))}{2\sqrt{z}}$, we have the same initial probability distribution but with different names \begin{align} P(g) = \begin{cases} \frac{1}{c} &\text{ when } \int_0^1 dz \; g(z) = 1\; ; f(z) \geq 0\\ 0 &\text{ otherwise} \end{cases} \end{align} This distrbution is again defined over $L^1[0,1]$ (because whenever $g(z)$ is integrable, so is $f(x(z))$ and vice versa).
This means that I should get the same result as above but for the function $g(z)$ i.e. \begin{align} \bar{g}(z) = 1 &\Rightarrow \frac{\bar{f}(\sqrt(z))}{2\sqrt{z}} = 1 \\ &\Rightarrow \bar{f}(x) = 2 x \end{align} which is different from the previous result $\bar{f}(x)=1$!
That was just one example. In general, when I do a variable transformation from $x$ to $z$, the result is \begin{align} \bar{f}(x) = \frac{dz}{dx} \text{ (up to a normalization factor)} \end{align} So the results depend on the variable that the problem is expressed in. When I performed Monte Carlo simulations on different nonuniform girds (which correspond to different variable transformations) I did get different results and the results confirms the above relation.
It looks like that "a uniform probability distribution over all positive functions on the unit interval whose Lebesgue integral is one" is not mathematically well-defined because its mean depends on the way I compute it.
I understand that the source of these different results is the hidden assumption about the symmetry of all the points in the interval $[0,1]$. Since the points get compressed or stretched upon variable transformation, I get different results.
Is there a mathematically rigorous way of expressing this assumption and including it in the definition of the probability distribution?