Let $X$ be a Gaussian random variable with mean=variance=$g>0$, namely its probability density is $p(x)=\frac{1}{\sqrt{2\pi g}}\,e^{-\frac{(x-g)^2}{2g}}$. I am looking for proofs of the following fact: $$ I(g):=\mathbb E[(1-3\tanh^2 X)\,(1-\tanh^2 X)^2] \,\geq\,0 .$$
I have one proof but I am looking also for other ones. My proof relies on the fact that $$ i:=\int_0^1 (1-3y^2) (1-y^2)^\frac{1}{2} dy \,>\, 0$$ since it can be computed explicitly. Now changing variable $\tanh(x)=y$ and expanding the square in the gaussian density you can rewrite $$ I(g) = \frac{e^{-\frac{g}{2}}}{\sqrt{2\pi g}}\, \int_{-1}^1 (1-3y^2)\,(1-y^2)^2 \,e^{\operatorname{atanh}(y)}\,e^{-\frac{1}{2g}\operatorname{atanh}^2(y)} \frac{d y}{1-y^2}$$ then, using $e^{\operatorname{atanh}(y)}+e^{\operatorname{atanh}(-y)}=\frac{2}{\sqrt{1-y^2}}$, you get $$ I(g) = \frac{2\,e^{-\frac{g}{2}}}{\sqrt{2\pi g}}\, \int_{0}^1 (1-3y^2)\,(1-y^2)^\frac{1}{2} \,e^{-\frac{1}{2g}\operatorname{atanh}^2(y)} \,d y $$ Finally, using the fact that $e^{-\frac{1}{2g}\operatorname{atanh}(y)^2}$ is monotonically decreaing in $y\in[0,1]$, you can bound with the value at the point $y=\frac{1}{\sqrt{3}}$ where the integrand changes sign: $$ I(g) \geq \frac{2\,e^{-\frac{g}{2}}}{\sqrt{2\pi g}}\, e^{-\frac{1}{2g}\operatorname{atanh}^2(\frac{1}{\sqrt{3}})} \ \int_{0}^1 (1-3y^2)\,(1-y^2)^\frac{1}{2} \,d y$$ and you conclude since you know $i>0$.
