"Since $x$ is near-gaussian, its PDF can be written as..." Short question: Why is this true??
Long question:
Very simply, I am trying to figure out what justifies this first equation. The author of the book I am reading, (context here if you want it, but not necessary), claims the following:

Due to the assumption of near-gaussianity, we can write:
$$ p_0(\xi) = A  \; \phi(\xi) \; exp( a_{n+1}\xi + (a_{n+2} +
 \frac{1}{2})\xi^2 + \sum_{i=1}^{n} a_i G_i(\xi)) $$

Where $p_0(\xi)$ is the PDF of your observed data that has maximum entropy, given that you had only observed a series of expectations, (simple numbers) $c_i, i = 1 ... n$, where $c_i = \mathbb{E}\{G_i(\xi)\}$, and $\phi(\xi)$ is the PDF of a standardized gaussian variable, that is, 0 mean, and unit variance. 
Where all this is going is that he uses the above equation as a starting point for making the PDF, $p_0(\xi)$ simpler, and I get how he does it, but I do not get how he justifies the above equation, ie, the starting point. 
I have tried to keep it brief to as not to obfuscate anyone, but if you want additional details please let me know in the comments. Thanks!
 A: (Note: I've changed your $\xi$ to $x$.)
For a random variable $X$ with density $p$, if you have constraints
$$
  \int G_i(x)\,p(x)\,dx=c_i \, ,
$$
for $i=1,\dots,n$, the maximum entropy density is
$$
  p_0(x)=A\exp\left(\sum_{i=1}^n a_iG_i(x)\right) \, ,
$$
where the $a_i$'s are determined from the $c_i$'s, and $A$ is a normalization constant. 
In this context, the Gaussian approximation ("near-gaussianity") means two things:
1) You accept to introduce two new constraints: the mean of $X$ is $0$ and the variance is $1$ (say);
2) The corresponding $a_{n+2}$ (see bellow) is much bigger than the other $a_i$'s.
These additional constraints are represented as
$$
  G_{n+1}(x)=x \, , \qquad c_{n+1}=0 \, ,
$$
$$
  G_{n+2}(x)=x^2 \, , \qquad c_{n+2}=1 \, ,
$$
yielding
$$
  p_0(x)=A\exp\left(a_{n+2}x^2 + a_{n+1}x + \sum_{i=1}^n a_iG_i(x)\right) \, ,
$$
which can be rewritten as (just "add zero" to the exponent)
$$
  p_0(x)=A\exp\left(\frac{x^2}{2} - \frac{x^2}{2} + a_{n+2}x^2 + a_{n+1}x + \sum_{i=1}^n a_iG_i(x)\right) \, ,
$$
leading to what you want:
$$
  p_0(x)=A'\,\phi(x)\exp\left(a_{n+1}x + \left(a_{n+2}+\frac{1}{2}\right)x^2 + \sum_{i=1}^n a_iG_i(x)\right) \, ;
$$
ready to be Taylor expanded (using the second condition of the Gaussian approximation).
Doing the approximation like a Physicist (which means that we don't care about the order of the error term), using $\exp(t)\approx 1+t$, we have the approximate density
$$
  p_0(x) \approx A'\,\phi(x)\left(1+a_{n+1}x + \left(a_{n+2}+\frac{1}{2}\right)x^2 + \sum_{i=1}^n a_iG_i(x)\right) \, .
$$
To finish, we have to determine $A'$ and the values of the $a_i$'s. This is done imposing the conditions
$$
  \int p_0(x)\,dx=1 \, , \qquad \int x \,p_0(x)\,dx=0 \, , \qquad \int x^2 \,p_0(x)\,dx=1
$$
$$
\int G_i(x)\, p_0(x)\,dx=c_i \, , \quad i=1,\dots,n \, ,
$$
to obtain a system of equations, whose solution gives $A'$ and the $a_i$'s.
Without imposing additional conditions on the $G_i$'s, I don't believe that there is a simple solution in closed form.
P.S. Mohammad clarified during a chat that with additional orthogonality conditions for the $G_i$'s we can solve the system.
