This generic function fits a semiparametric model, which consists of parametric and nonparametric, for density estimation (Lenk, 2003):
$$f(x | \beta, Z) = \frac{\exp[h(x)^\top\beta + Z(x)]}{\int_\mathcal{X} \exp[h(y)^\top\beta + Z(y)]dG(y)}$$
where \(Z\) is a zero mean, second-order Gaussian process with bounded, continuous covariance function. i.e.,
$$E[Z(x), Z(y)] = \sigma(x,y), \quad \int_\mathcal{X}ZdG = 0 ~~(a.s.)$$
Using the Karhunen-Loeve Expansion, \(Z\) is represented as infinite series with random coefficients
$$Z(x) = \sum_{j=1}^\infty \theta_j\varphi_j(x), $$
where \(\{\varphi_j\}\) is the cosine basis, \(\varphi_j(x)=\sqrt{2}\cos[j\pi G(x)]\).
For the random Fourier coefficients of the expansion, two smoother priors are assumed (optional),
$$\theta_j | \tau, \gamma \sim N(0, \tau^2\exp[-j\gamma]), ~ j \ge 1 ~ (geometric ~smoother)$$
$$\theta_j | \tau, \gamma \sim N(0, \tau^2\exp[-ln(j+1)\gamma]), ~ j \ge 1 ~ (algebraic ~smoother)$$
The coefficient \(\beta\) have the popular normal prior,
$$\beta | m_{0,\beta}, V_{0,\beta} \sim N(m_{0,\beta}, V_{0,\beta})$$
To complete the model specification, independent hyper priors are assumed,
$$\tau^2 | r_0, s_0 \sim IGa(r_0/2, s_0/2)$$
$$\gamma | w_0 \sim Exp(w_0)$$
Note that the posterior algorithm is based on computing a discrete version of the likelihood over a fine mesh on \(\mathcal{X}\).