def gfit(X, sigma, p=<NUM_LIT:5>, nbin=<NUM_LIT:200>, unif_fraction=<NUM_LIT:0.1>):
min_x = min(min(X) - <NUM_LIT:2> * np.std(X, ddof=<NUM_LIT:1>), <NUM_LIT:0>)<EOL>max_x = max(max(X) + <NUM_LIT:2> * np.std(X, ddof=<NUM_LIT:1>),<EOL>np.std(X, ddof=<NUM_LIT:1>))<EOL>xvals = np.linspace(min_x, max_x, nbin)<EOL>binw = (max_x - min_x) / (nbin - <NUM_LIT:1>)<EOL>zero_idx = max(np.where(xvals <= <NUM_LIT:0>)[<NUM_LIT:0>])<EOL>noise_kernel = norm().pdf(xvals / sigma) * binw / sigma<EOL>if zero_idx > <NUM_LIT:0>:<EOL><INDENT>noise_rotate = noise_kernel[list(np.arange(zero_idx, len(xvals))) +<EOL>list(np.arange(<NUM_LIT:0>, zero_idx))]<EOL><DEDENT>else:<EOL><INDENT>noise_rotate = noise_kernel<EOL><DEDENT>XX = np.zeros((p, len(xvals)), dtype=np.float)<EOL>for ind, exp in enumerate(range(<NUM_LIT:1>, p+<NUM_LIT:1>)):<EOL><INDENT>mask = np.ones_like(xvals)<EOL>mask[np.where(xvals <= <NUM_LIT:0>)[<NUM_LIT:0>]] = <NUM_LIT:0><EOL>XX[ind, :] = pow(xvals, exp) * mask<EOL><DEDENT>XX = XX.T<EOL>def neg_loglik(eta):<EOL><INDENT>mask = np.ones_like(xvals)<EOL>mask[np.where(xvals <= <NUM_LIT:0>)[<NUM_LIT:0>]] = <NUM_LIT:0><EOL>g_eta_raw = np.exp(np.dot(XX, eta)) * mask<EOL>if ((np.sum(g_eta_raw) == np.inf) |<EOL>(np.sum(g_eta_raw) <=<EOL><NUM_LIT:100> * np.finfo(np.double).tiny)):<EOL><INDENT>return (<NUM_LIT:1000> * (len(X) + sum(eta ** <NUM_LIT:2>)))<EOL><DEDENT>g_eta_main = g_eta_raw / sum(g_eta_raw)<EOL>g_eta = ((<NUM_LIT:1> - unif_fraction) * g_eta_main +<EOL>unif_fraction * mask / sum(mask))<EOL>f_eta = fftconvolve(g_eta, noise_rotate, mode='<STR_LIT>')<EOL>return np.sum(np.interp(X, xvals,<EOL>-np.log(np.maximum(f_eta, <NUM_LIT>))))<EOL><DEDENT>eta_hat = minimize(neg_loglik,<EOL>list(itertools.repeat(-<NUM_LIT:1>, p))).x<EOL>g_eta_raw = np.exp(np.dot(XX, eta_hat)) * mask<EOL>g_eta_main = g_eta_raw / sum(g_eta_raw)<EOL>g_eta = ((<NUM_LIT:1> - unif_fraction) * g_eta_main +<EOL>unif_fraction * mask) / sum(mask)<EOL>return xvals, g_eta<EOL>
Fit empirical Bayes prior in the hierarchical model [Efron2014]_.

.. math::

    mu ~ G, X ~ N(mu, sigma^2)

Parameters
----------
X: ndarray
    A 1D array of observations.
sigma: float
    Noise estimate on X.
p: int
    Number of parameters used to fit G. Default: 5
nbin: int
    Number of bins used for discrete approximation.
    Default: 200
unif_fraction: float
    Fraction of G modeled as "slab". Default: 0.1

Returns
-------
An array of the posterior density estimate g.