I have a function $f:\mathbb{R}^+ \to \mathbb{R}^+$ that I want to prove is eventually concave - i.e. that there exists $\gamma _0 > 0$ such that for every $\gamma>\gamma_0$, $f(\gamma)$ is concave. I also know that *$f$ is real-analytic, increasing and upper bounded*, and it is therefore enough to show that $f''$ has a finite number of zeros.

A concrete definition of $f(\gamma)$ is as follows: $$f(\gamma)= \frac{1}{\sqrt{2\pi}}\sum_{x\in\mathcal{X},x'\in\mathcal{X}'}p_{x,x'}\int_{-\infty}^\infty{e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}\log\left(\frac{\sum_{x\in\mathcal{X}}p_{x|x'}e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}}{\sum_{x\in\mathcal{X}}p_{x}e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}}\right)dy}$$ where $\mathcal{X},\mathcal{X}'$ are two finite sets, $p:\mathcal{X}\times\mathcal{X}'\to[0,1]$ is a probability mass function and $p_{x|x'}=p_{x,x'}/\sum_{\xi\in\mathcal{X}}p_{\xi,x'}$ , $p_x=\sum_{\xi'\in\mathcal{X}'}p_{x,\xi'}$ are the conditional and marginal distributions on $\mathcal{X}$, respectively. In Information-Theoretic terms, $f(\gamma)$ is the mutual information between $X'$ and $Y_{\gamma}=\sqrt{\gamma}X+N$, where $N \sim \mathcal{N}(0,1)$ is independent of $X'$ and $X$, which are dependent and have finite alphabets. The variable $\gamma$ is the signal-to-noise ratio.

The function is known to be concave in the special case $X=X'$ ($p_{x|x'}=\delta_{x,x'}$). When $X\neq X'$, examples can be found in which $f$ is not always concave, but becomes concave eventually as $\gamma$ grows.

Intuitively, $f''$ must have finitely many zeros (=$f$ is eventually concave), as $f$ is a finite combination of "well-behaved" functions. I would like to formalize this by showing that $f$ belongs to class of functions that have a finite number of zeros, that is closed under differentiation. This will show that $f''$ also has finitely many zeros, which is just the thing we need.

There exist several such classes of functions including Hardy's class L of "orders of infinity" and the Pfaffian closure. See also this related question. However, the mathematical machinery here is a bit complicated for me, and I could not show $f$ belongs to either class due to the integration with resepect to $y$ in its definition. Perhaps this can be shown by recasting $f$ as a solution to an ODE? Or perhaps there exists another, more suitable class of functions that I am not aware of?

There is also an expression for $f''(\gamma)$ which is quite lengthy. Here it is for completeness: $$f''(\gamma) = -\frac{1}{2\sqrt{2\pi}}\sum_{x\in\mathcal{X},x'\in\mathcal{X}'}p_{x,x'}\int_{-\infty}^\infty{e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}[\phi_X^2(y;\gamma)-\phi_{X|X'}^2(y,x';\gamma)]dy}$$ where $$\phi_X(y;\gamma)=\frac{\sum_{x\in\mathcal{X}}x^2p_{x}e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}}{\sum_{x\in\mathcal{X}}p_{x}e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}}-\left(\frac{\sum_{x\in\mathcal{X}}xp_{x}e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}}{\sum_{x\in\mathcal{X}}p_{x}e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}}\right)^2$$ and similarly $$\phi_{X|X'}(y,x';\gamma)=\frac{\sum_{x\in\mathcal{X}}x^2p_{x|x'}e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}}{\sum_{x\in\mathcal{X}}p_{x|x'}e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}}-\left(\frac{\sum_{x\in\mathcal{X}}xp_{x|x'}e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}}{\sum_{x\in\mathcal{X}}p_{x|x'}e^{-\frac{1}{2}(y-\sqrt{\gamma}x)^2}}\right)^2$$ This has some estimation-theoretic meaning, but not anything that seems useful in establishing my desired result.