Consider a random function $\phi(\underline{x})$, where $\underline{x}$ is a spatial position. For each independent realization $\phi^{(k)}(\underline{x})$ of the function there is a certain fraction $\xi^{(k)}$ of the domain that is mapped to the interval $\phi_0 < \phi < \phi_0 + \Delta \phi$. Since the function is random, this fraction is a random variable. My hypothesis is that the expected fraction is given by $$ \langle \xi \rangle = \Delta \phi f(\phi_0) + \mathcal{O}(\Delta \phi^2), \tag{1} $$ where $f(\psi)$ is the probability density function (PDF) of $\phi$; the probability density of the event $\phi(\underline{x}) = \psi$. Is the hypothesis correct? I would be glad if you could explain how this can be proven (or disproved).
The functions $\phi(\underline{x})$ that I am interested in are smooth and posses no singularities, so feel free to assume that random functions that are well-behaved in these senses. Furthermore, the functions are statistically homogeneous, which means that their statistics are independent of $\underline{x}$.
I have proven Equation (1) for a one-dimensional periodic random function, by expressing the function as a Fourier series. I provide the proof below. Please remark on it if you find it incorrect. The proof can be straightforwardly extended to periodic problems in higher dimensions. However, it seems non-intuitive to me that periodicity should be a necessary criterion for Equation (1) to hold. What about random functions $\phi(\underline{x})$ on arbitrary two- or three-dimensional domains, for example?
Proof for the periodic one-dimensional case
I denote the spatial position $\underline{x}$ by $x$ for the one-dimensional case considered in the following. Assuming that the random function is periodic over the period $L$, we have $\phi(x) = \phi(x + L)$. For a given realization $\phi^{(k)}(x)$ of the function, one may sample a random variable $y$ in the inter in the interval $0 \leq y < L$. The corresponding value of $\phi^{(k)}(y)$ is another random variable. I denote the cumulative distribution function of $\phi^{(k)}(y)$ by $F^{(k)}(\psi)$; it is the probability of sampling $\phi^{(k)}(y) < \psi$. I denote the corresponding PDF by $f^{(k)}(\psi) = \partial/\partial \psi F^{(k)}(\psi)$. It is clear that the fraction $\xi^{(k)}$ of the domain that is mapped to the interval $\phi_0 < \phi < \phi_0 + \Delta \phi$ is given by $$ \xi^{(k)} = F^{(k)}(\phi_0 + \Delta\phi) - F^{(k)}(\phi_0) = \Delta \phi f^{(k)}(\phi_0) + \mathcal{O}(\Delta \phi^2). \tag{2} $$
A given realization $\phi^{(k)}(x)$ of the random function can be expressed as a Fourier series,
$$ \phi^{(k)}(x) = \frac{A_0^{(k)}}{2} + \sum_{n = 1}^{\infty} A_n^{(k)}\cos \left ( \frac{2\pi x}{L} - B_n^{(k)} \right )\tag{3} $$
with random amplitudes $A_n^{(k)}$ and phases $B_n^{(k)}$. I denote the infinite sets of amplitudes and phases for the realization $\phi^{(k)}(x)$ by $\underline{A}^{(k)} = \{ A_0^{(k)}, A_1^{(k)}, A_2^{(k)}, \dots \}$ and $\underline{B}^{(k)} = \{ B_1^{(k)}, B_2^{(k)}, \dots \}$. These sets determine the function $\phi^{(k)}(x)$, and one may therefore write
$$ f^{(k)}(\psi) = f_{\phi|\underline{A}\underline{B}}(\psi|\underline{A}^{(k)},\underline{B}^{(k)}).\tag{4} $$
Here, $f_{\phi|\underline{A}\underline{B}}(\psi|\underline{\alpha},\underline{\beta})$ denotes the PDF of $\phi$, conditioned upon the that the sets $\underline{A}^{(k)}$ and $\underline{B}^{(k)}$ are given by $\underline{\alpha}$ and $\underline{\beta}$. This PDF is given by
$$ f_{\phi|\underline{A}\underline{B}}(\psi|\underline{\alpha},\underline{\beta}) = \frac{f_{\phi\underline{A}\underline{B}}(\psi,\underline{\alpha},\underline{\beta})}{f_{\underline{A}\underline{B}}(\underline{\alpha},\underline{\beta})},\tag{5} $$ where $f_{\phi\underline{A}\underline{B}}(\psi,\underline{\alpha},\underline{\beta})$ is the probability density of the event that $\phi = \psi$, and that the amplitudes and phases of the random function are given by $\underline{\alpha}$ and $\underline{\beta}$. Similarly, $f_{\underline{A}\underline{B}}(\underline{\alpha},\underline{\beta})$ is the probability density of the event that the random function has the amplitudes and phases $\underline{\alpha}$ and $\underline{\beta}$.
Using Equations (2) and (5), one finds $$ \xi^{(k)}f_{\underline{A}\underline{B}}(\underline{A}^{(k)},\underline{B}^{(k)}) = \Delta \phi f_{\phi\underline{A}\underline{B}}(\phi_0,\underline{A}^{(k)},\underline{B}^{(k)}) + \mathcal{O}(\Delta \phi^2). \tag{6} $$ One obtains Equation (1) by integrating this equation over all amplitudes $A_n^{(k)}$ and phases $B_n^{(k)}$: $$ \langle \xi \rangle = \int_{-\infty}^\infty \int_{-\infty}^\infty \int_{-\infty}^\infty \int_{-\infty}^\infty \int_{-\infty}^\infty \quad\dots \quad \xi^{(k)}f_{\underline{A}\underline{B}}(\underline{A}^{(k)},\underline{B}^{(k)})\quad \mathrm d A_0^{k} \mathrm d A_1^{k} \mathrm d B_1^{k} \mathrm d A_2^{k} \mathrm d B_2^{k}\quad \dots\\ = \Delta \phi \int_{-\infty}^\infty \int_{-\infty}^\infty \int_{-\infty}^\infty \int_{-\infty}^\infty \int_{-\infty}^\infty \quad\dots \quad f_{\phi\underline{A}\underline{B}}(\phi_0,\underline{A}^{(k)},\underline{B}^{(k)})\quad \mathrm d A_0^{k} \mathrm d A_1^{k} \mathrm d B_1^{k} \mathrm d A_2^{k} \mathrm d B_2^{k}\quad \dots \quad + \mathcal{O}(\Delta \phi^2) \\ = \Delta \phi f(\phi_0) + \mathcal{O}(\Delta \phi^2). \quad \square \tag{7} $$
After thinking for a while, I was able to prove the hypothesis [Equation (1)]. The breakthrough came when I started to think of the sampling of a value from the random function $\phi(\underline{x})$ as a two-stage process. First, one samples one realization of the random function. After that, one samples a value from the realization.
I now drop the superscript $^{(k)}$, and simply denote the fraction of the domain that a random function maps to $\phi_0 < \phi < \phi_0 + \Delta \phi$ by $\xi$. If the random function $\phi(\underline{x})$ maps a fraction $\xi$ of the domain to this interval, it holds that $$ \xi = F_{\phi|\xi}(\phi_0 + \Delta \phi|\xi) - F_{\phi|\xi}(\phi_0|\xi) = \Delta \phi f_{\phi|\xi}(\phi_0|\xi) + \mathcal{O}(\Delta \phi^2). \tag{P1} $$ Here, $F_{\phi|\xi}(\psi|\zeta)$ is the CDF of $\phi$ conditioned upon $\xi = \zeta$; the probability that $\phi < \psi$ if $\xi = \zeta$. Furthermore, $f_{\phi|\xi}(\psi|\zeta) = \partial/\partial \psi F_{\phi|\xi}(\psi|\zeta)$ is the the corresponding PDF. Equations (P1) is similar to Equation (2), but does not describe how $\xi$ is computed for a single realization of the random function. Instead, it describes how $\xi$ relates to the conditional PDF $f_{\phi|\xi}(\psi|\zeta)$.
I denote the PDF of $\xi$ by $f_\xi(\zeta)$; it gives the probability density of the event $\xi = \zeta$. The average fraction $\langle \xi \rangle$ is obtained by multiplying Equation (P1) by $f_\xi(\xi)$, and integrating over all fractions $\xi$:
$$ \langle \xi \rangle = \int \xi f_\xi(\xi) \mathrm d \xi\\ = \Delta \phi \int f_{\phi|\xi}(\phi_0|\xi)f_\xi(\xi) \mathrm d \xi + \mathcal{O}(\Delta \phi^2)\\ = \Delta \phi f(\phi_0) + \mathcal{O}(\Delta \phi^2). \quad \square \tag{P2} $$
This proof is much more elegant and general than my old proof for a one-dimensional periodic function. I believe that it is correct, but please correct me if I have made a mistake.
I find it interesting that I had been thinking about how to prove Equation (1) for some days, but it was only when I wrote about the question here that I had a breakthrough. I feel that writing something that an audience may read (you), catalyzed my thinking. So thank you! :)