Here is a very partial answer.
Theorem : if ${\mathbf{A}}$ is positive semidefinite and ${\mathbf{J}} = {\mathbf{A}}{{\mathbf{A}}^ + }{\mathbf{J}}$, then $\mathbb{E}{\mathbf{x}} = {{\mathbf{A}}^ + }{\mathbf{J}}$ , provided we allow ourself to cancel out terms like $\frac{a}{a}$ even if $a$ is infinite.
Proof : recall one of the usual proofs of the identity
$\int\limits_{{\mathbb{R}^n}} {{e^{ - \frac{1}{2}{{\mathbf{x}}^{\mathbf{T}}}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}}}}{\text{d}}{\mathbf{x}}} = \frac{{{{\left( {2\pi } \right)}^{\frac{n}{2}}}}}{{\sqrt {\left| {\mathbf{A}} \right|} }}{e^{\frac{1}{2}{{\mathbf{J}}^T}{{\mathbf{A}}^{ - 1}}{\mathbf{J}}}}$
for a positive definite matrix ${\mathbf{A}}$ .
Substitute ${\mathbf{x}}$ by ${\mathbf{y}} = {\mathbf{x}} - {{\mathbf{A}}^{ - 1}}{\mathbf{J}}$
${{\text{d}}^n}{\mathbf{x}} = {{\text{d}}^n}{\mathbf{y}}$
$ - \frac{1}{2}{{\mathbf{x}}^T}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}} = - \frac{1}{2}{\left( {{\mathbf{y}} + {{\mathbf{A}}^{ - 1}}{\mathbf{J}}} \right)^T}{\mathbf{A}}\left( {{\mathbf{y}} + {{\mathbf{A}}^{ - 1}}{\mathbf{J}}} \right) + {{\mathbf{J}}^T}\left( {{\mathbf{y}} + {{\mathbf{A}}^{ - 1}}{\mathbf{J}}} \right) = \\
- \frac{1}{2}\left( {{{\mathbf{J}}^T}{{\mathbf{A}}^{ - 1}}{\mathbf{J}} + {{\mathbf{J}}^T}{\mathbf{y}} + {{\mathbf{y}}^T}{\mathbf{J}} + {{\mathbf{y}}^T}{\mathbf{Ay}}} \right) + {{\mathbf{J}}^T}{{\mathbf{A}}^{ - 1}}{\mathbf{J}} + {{\mathbf{J}}^T}{\mathbf{y}} = \\
\frac{1}{2}{{\mathbf{J}}^T}{{\mathbf{A}}^{ - 1}}{\mathbf{J}} - \frac{1}{2}{{\mathbf{y}}^T}{\mathbf{Ay}} \\ $
Therefore
$\int\limits_{{\mathbb{R}^n}} {{e^{ - \frac{1}{2}{{\mathbf{x}}^T}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}}}}{\text{d}}{\mathbf{x}}} = {e^{\frac{1}{2}{{\mathbf{J}}^T}{{\mathbf{A}}^{ - 1}}{\mathbf{J}}}}\int\limits_{{\mathbb{R}^n}} {{e^{ - \frac{1}{2}{{\mathbf{y}}^T}{\mathbf{Ay}}}}{\text{d}}{\mathbf{y}}} = \frac{{{{\left( {2\pi } \right)}^{\frac{n}{2}}}}}{{\sqrt {\left| {\mathbf{A}} \right|} }}{e^{\frac{1}{2}{{\mathbf{J}}^T}{{\mathbf{A}}^{ - 1}}{\mathbf{J}}}}$
Then, by the Leibniz rule/Feynman trick
$ \frac{{\partial \int\limits_{{\mathbb{R}^n}} {{e^{ - \frac{1}{2}{{\mathbf{x}}^T}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}}}}{\text{d}}{\mathbf{x}}} }}{{\partial {{\mathbf{J}}_i}}} = \int\limits_{{\mathbb{R}^n}} {\frac{{\partial {e^{ - \frac{1}{2}{{\mathbf{x}}^T}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}}}}}}{{\partial {{\mathbf{J}}_i}}}{\text{d}}{\mathbf{x}}} = \int\limits_{{\mathbb{R}^n}} {{x_i}{e^{ - \frac{1}{2}{{\mathbf{x}}^T}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}}}}{\text{d}}{\mathbf{x}}} = \\
\frac{{{{\left( {2\pi } \right)}^{\frac{n}{2}}}}}{{\sqrt {\left| {\mathbf{A}} \right|} }}\frac{{\partial {e^{\frac{1}{2}{{\mathbf{J}}^T}{{\mathbf{A}}^{ - 1}}{\mathbf{J}}}}}}{{\partial {{\mathbf{J}}_i}}} = \frac{1}{2}\int\limits_{{\mathbb{R}^n}} {{e^{ - \frac{1}{2}{{\mathbf{x}}^T}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}}}}{\text{d}}{\mathbf{x}}} \frac{\partial }{{\partial {{\mathbf{J}}_i}}}{{\mathbf{J}}^T}{{\mathbf{A}}^{ - 1}}{\mathbf{J}} \\ $
Hence
$\mathbb{E}{x_i} = \frac{{\int\limits_{\,{\mathbf{x}}} {{x_i}{e^{ - \frac{1}{2}{{\mathbf{x}}^T}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}}}}{\text{d}}{\mathbf{x}}} }}{{\int\limits_{\,{\mathbf{x}}} {{e^{ - \frac{1}{2}{{\mathbf{x}}^T}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}}}}{\text{d}}{\mathbf{x}}} }} = \frac{1}{2}\frac{\partial }{{\partial {{\mathbf{J}}_i}}}{{\mathbf{J}}^T}{{\mathbf{A}}^{ - 1}}{\mathbf{J}} = {\left( {{{\mathbf{A}}^{ - 1}}{\mathbf{J}}} \right)_i}$
Now, for a positive semidefinite matrix ${\mathbf{A}}$, substitute ${\mathbf{x}}$ by ${\mathbf{y}} = {\mathbf{x}} - {{\mathbf{A}}^ + }{\mathbf{J}}$
$ - \frac{1}{2}{{\mathbf{x}}^T}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}} = - \frac{1}{2}{\left( {{\mathbf{y}} + {{\mathbf{A}}^ + }{\mathbf{J}}} \right)^T}{\mathbf{A}}\left( {{\mathbf{y}} + {{\mathbf{A}}^ + }{\mathbf{J}}} \right) + {{\mathbf{J}}^T}\left( {{\mathbf{y}} + {{\mathbf{A}}^ + }{\mathbf{J}}} \right) = \\
- \frac{1}{2}\left( {{{\mathbf{J}}^T}\underbrace {{{\mathbf{A}}^ + }{\mathbf{A}}{{\mathbf{A}}^ + }}_{{{\mathbf{A}}^ + }}{\mathbf{J}} + {{\mathbf{J}}^T}{{\mathbf{A}}^ + }{\mathbf{Ay}} + {{\mathbf{y}}^T}{\mathbf{A}}{{\mathbf{A}}^ + }{\mathbf{J}} + {{\mathbf{y}}^T}{\mathbf{Ay}}} \right) + {{\mathbf{J}}^T}{{\mathbf{A}}^ + }{\mathbf{J}} + {{\mathbf{J}}^T}{\mathbf{y}} = \\
\frac{1}{2}{{\mathbf{J}}^T}{{\mathbf{A}}^ + }{\mathbf{J}} - \frac{1}{2}{{\mathbf{y}}^T}{\mathbf{Ay}} + {{\mathbf{J}}^T}\left( {{\mathbf{I}} - {{\mathbf{A}}^ + }{\mathbf{A}}} \right){\mathbf{y}} \\ $
The integral
$\int\limits_{\,{\mathbf{x}}} {{e^{ - \frac{1}{2}{{\mathbf{y}}^T}{\mathbf{Ay}}}}{\text{d}}{\mathbf{y}}} $
now is infinite. But it is not a big deal because it cancels out in the Leibniz rule/Feynman trick above (please tell me).
Therefore, the term ${{\mathbf{J}}^T}\left( {{\mathbf{I}} - {{\mathbf{A}}^ + }{\mathbf{A}}} \right){\mathbf{y}}$ , where ${\mathbf{I}} - {{\mathbf{A}}^ + }{\mathbf{A}}$ is the orthogonal projector on $\ker {\mathbf{A}}$, is the main obstruction against the generalized formula.
So, if ${{\mathbf{J}}^T}\left( {{\mathbf{I}} - {{\mathbf{A}}^ + }{\mathbf{A}}} \right) = 0 \Leftrightarrow {\mathbf{J}} = {\mathbf{A}}{{\mathbf{A}}^ + }{\mathbf{J}}$ then
$\int\limits_{{\mathbb{R}^n}} {{e^{ - \frac{1}{2}{{\mathbf{x}}^T}{\mathbf{Ax}} + {{\mathbf{J}}^T}{\mathbf{x}}}}{\text{d}}{\mathbf{x}}} \propto {e^{\frac{1}{2}{{\mathbf{J}}^T}{{\mathbf{A}}^ + }{\mathbf{J}}}}$
and the generalized formula
$\mathbb{E}{\mathbf{x}} = {{\mathbf{A}}^ + }{\mathbf{J}}$
follows by the Leibniz rule/Feynman trick.
Perhaps this condition is fulfilled with my own ${\mathbf{A}}$ and ${\mathbf{J}}$, I need to check.