How can I estimate the following expectation \begin{equation} \Big\langle\log\big(\Gamma(x)\big)\Big\rangle_{q(x)} \end{equation} let's assume that the distributions for $\alpha$ is a gamma distribution $q(x)=\mathrm{Gamma}(a,b)$?
If we first assume $b=1$ then we have \begin{equation} \int \ln\Gamma(x)\frac{1}{\Gamma(a)}x^a e^{-x}\frac{dx}{x} \end{equation} and maybe we can substitute $x$ with $e^z$ so we get \begin{equation} \int \ln\Gamma(e^z)\frac{1}{\Gamma(a)} e^{az-e^z}dz \end{equation} Is there any analytical way to compute this integral?
I'll assume the region of integration is $[0,\infty)$. The previous form you had for the integral was actually better since we can identify it as $$\frac{1}{\Gamma(s)}\int_0^\infty x^{s-1}e^{-x}\ln\Gamma(x) \mathrm{d}x=\frac{1}{\Gamma(s)}\mathcal{M}\{e^{-x}\ln\Gamma(x)\}(s)$$ However this is unfortunately still quite bothersome. So perhaps we can use the series expansion of the log Gamma: $$\ln\Gamma(x)=-\gamma x-\ln x+\sum_{k=1}^\infty \left[\frac{x}{k}-\ln\left(1+\frac{x}{k}\right)\right]$$ With of course $\gamma$ the Euler-Mascheroni constant. Now we can express our integral as $$\mathcal{M}\{e^{-x}\ln\Gamma(x)\}(s)$$ $$=\int_0^\infty-\gamma x^{s-1}xe^{-x}\mathrm{d}x-\int_0^\infty x^{s-1}e^{-x}\ln x~\mathrm{d}x+\int_0^\infty x^{s-1} e^{-x}\sum_{k=1}^\infty \left[\frac{x}{k}-\ln\left(1+\frac{x}{k}\right)\right]\mathrm{d}x$$ The first integral is very simple (just the Gamma function). The second integral can also be dealt with fairly simply. Notice that $$\int_0^\infty x^{s-1}e^{-x}\ln(x)\mathrm{d}x=\int_0^\infty \partial_s[x^{s-1}e^{-x}]\mathrm{d}x=\Gamma'(s)=\Gamma(s)\psi^{[0]}(s)$$ By the Leibniz rule, and using the usual definition of the polygamma. Thus we now have (wow, a lot of "gamma" here... the Gamma, the Euler gamma, the polygamma, oh my!) $$\mathcal{M}\{e^{-x}\ln\Gamma(x)\}(s)=-\gamma\cdot \Gamma(s+1)-\Gamma(s)\psi^{[0]}(s)+\int_0^\infty\sum_{k=1}^\infty \left[\frac{x^s e^{-x}}{k}-x^{s-1}e^{-x}\ln\left(1+\frac{x}{k}\right)\right]\mathrm{d}x$$ Using the recurrence relation of the Gamma and hoping and praying we can interchange summation and integration, $$\mathcal{M}\{e^{-x}\ln\Gamma(x)\}(s)=-\Gamma(s)(\gamma s+\psi^{[0]}(s))+\sum_{k=1}^\infty \int_0^\infty \frac{x^s e^{-x}}{k}-x^{s-1}e^{-x}\ln\left(1+\frac{x}{k}\right) \mathrm{d}x$$ $$\mathcal{M}\{e^{-x}\ln\Gamma(x)\}(s)+\Gamma(s)(\gamma s+\psi^{[0]}(s))$$ $$=\sum_{k=1}^\infty \left[\frac{1}{k}\Gamma(s+1)-\int_0^\infty x^{s-1}e^{-x}\ln\left(1+\frac{x}{k}\right)\mathrm{d}x\right]$$ The last integral is still pretty nasty though. I put it into Mathematica and here's what it came up with: $$I(s,k)=\int_0^\infty x^{s-1}e^{-x}\ln\left(1+\frac{x}{k}\right)\mathrm{d}x$$ $$=\frac{1}{(-k^2)^s}\bigg[k^s \csc(\pi s)\gamma_{\text{low}}(s,-k)+(-k)^s \bigg(k\Gamma(s-1){}_2 F_2([1,1],[2,2-s];k)+\Gamma(s)\big(\psi^{[0]}(s)-\ln(k)\big)\bigg)\bigg]$$ Which uses the lower incomplete Gamma and the generalized hypergeometric function.
There are probably some good ways to reach fast approximations of the above. This brings us to $$\frac{1}{\Gamma(s)}\mathcal{M}\{e^{-x}\ln\Gamma(x)\}(s)=-\gamma s-\psi^{[0]}(s)+\sum_{k=0}^\infty \left[\frac{s}{k}-\frac{1}{\Gamma(s)}I(s,k)\right]$$ Unfortunately this is pretty much the best we can do (I think). Perhaps there are some more elaborate Mellin transform identities you can use from the start, but I'm absolutely no expert. Hopefully this is a good start.