Unbiased estimator when only the magnitude can be measured given an random variable $x$, which is drawn from a normal distribution $f(x| \mu_x, \sigma_x) = \frac{1}{\sqrt{2 \pi \sigma_x^2}} \exp\left(-\frac{(x-\mu_x)^2}{2 \sigma_x^2}\right)$.
We are drawing $N$ samples $x_i$ from this distribution, but unfortunately, we can only measure their magnitude $|x_i|$.
Now, it is pretty obvious that we will not be able to estimate $\mu_x$ from our data, but I'd like to estimate $|\mu_x|$. The arithmetic mean of the $|x_i|$ is a pretty good estimator when $|\mu_x| \gg 0$, but when $\mu_x$ is close to zero (e.g. exactly zero), the arithmetic mean will be pretty biased (if $\mu_x$ is exactly zero, the arithmetic mean will be $\propto \sigma_x$).
Unfortunately, looking at the data I suspect that $\mu_x$ is actually pretty close to zero.
So, my question is: Where should I look for an unbiased estimator? I tried googling it, but I guess I didn't use the right search terms.
 A: The distribution to use for this problem is the folded normal distribution (https://en.wikipedia.org/wiki/Folded_normal_distribution). From this paper (https://www.tandfonline.com/doi/abs/10.1080/03610927408827103) we can get the likelihood function for this distribution, which is
$f(x; \mu_x, \sigma_x) = \frac{2}{\sqrt{2\pi}\sigma_x} \exp\left(-\frac{x^2+\mu_x^2}{2\sigma_x^2}\right) \cosh(x \mu_x / \sigma_x^2)$.
Using this likelihood function, I did a maximum likelihood estimation simply by brute-forcing the $\mu_x$ and $\sigma_x$. With sufficiently dense points for $\mu_x$ and $\sigma_x$, I could also calculate credible intervals.
Python code I used in the end (for the confidence intervals):
import numpy as np
import xarray as xr
import numba

@numba.jit(nopython=True)
def log_likelihood(x, mu, sigma):
    return np.log(2/(np.sqrt(2*np.pi)*sigma) * np.exp(-(x**2+mu**2)/(2*sigma**2)) * np.cosh(x*mu/sigma**2))

@numba.jit(nopython=True)
def _log_likelihood_sample(xi, mu, sigma):
    return log_likelihood(xi, mu, sigma).sum()

@numba.jit(nopython=True)
def _log_likelihood_map(xi, mus, sigmas, res):
    for i, mu in enumerate(mus):
        for j, sigma in enumerate(sigmas):
            res[i, j] = _log_likelihood_sample(xi, mu, sigma)

def log_likelihood_map(xi):
    mus = np.linspace(0, 10, 501)
    sigmas = np.linspace(0.1, 10, 500)
    res = np.empty((len(mus), len(sigmas)), dtype=np.float)
    _log_likelihood_map(xi, mus, sigmas, res)
    return xr.DataArray(res, coords=[mus, sigmas], dims=['mu', 'sigma'])

llm = log_likelihood_map(np.array([0.5, 0.7, 3., 2.1, 1.5))
llm = llm.where(np.isfinite(llm))
llm -= llm.max()
lm = np.exp(llm)
lm /= lm.sum()
cdf = lm.sum('sigma').cumsum('mu')
ci = np.interp(np.array([0.33, 0.5, 0.67]), cdf, cdf.mu)
print(ci)

