Given a discrete random variable $\xi\colon\Omega\rightarrow X$, $\left| X \right| < \infty$ with a probability distribution depending on a parameter $\theta\in\Theta=\Theta_1\times...\times\Theta_n$, $n\in\mathbb{N}_+$ and a function $f\colon X\times\Theta\rightarrow \mathbb{R} \cup \left\{ -\inf \right\}$ that is double differentiable by the second parameter, namely $$ P_{\xi}(x; \theta)=\frac{\exp{f\left(x, \theta\right)}}{Z\left(\theta\right)}, \\ Z\left(\theta\right) = \sum_{x\in X} \exp{f\left(x, \theta\right)}, $$ I would like to know whether its Shannon entropy is a convex function of $\theta$.
If my calculations are correct, then the second derivative of the entropy is $$ \frac{\partial^2 H_{\xi}\left(\theta\right)} {\partial\theta_j \partial\theta_i} = cov_{\xi}\left( \frac{\partial f\left( \xi, \theta \right)}{\partial\theta_j}, \frac{\partial f\left( \xi, \theta \right)}{\partial\theta_i} \cdot \left( E_{\xi}f\left( \xi, \theta \right) + f\left( \xi, \theta \right) + 1 \right) + f\left( \xi, \theta \right) \cdot E_{\xi} \frac{\partial f\left( \xi, \theta \right)}{\partial\theta_i} \right) \\ - cov_{\xi}\left( \frac{\partial^2 f\left( \xi, \theta \right)} {\partial\theta_j\partial\theta_i}, f\left( \xi, \theta \right) \right). $$ If the second derivative of $f$ is a constant, I still cannot tell anything about positive semidefiniteness unless I can tell in which cases $cov_{\xi}\left(\frac{\partial f\left( \xi, \theta \right)}{\partial\theta_i}, f\left( \xi, \theta \right) \right)$ and $cov_{\xi}\left(\frac{\partial f\left( \xi, \theta \right)}{\partial\theta_i}, \frac{\partial f\left( \xi, \theta \right)}{\partial\theta_i}\cdot f\left( \xi, \theta \right)\right)$ form positive semidefinite covariance matrices, because in this case I obtain a sum of positive semidefinite matrices, which is a positive semi-definite matrix.
Are there researches on this topic, am I losing something obvious, or is this a really complicated question depending on the specific case?