You can do it by essentially iterating (2). First you get your sigma estimate using (2). This is not so great as you've said. Now that you have a (bad) sigma estimate, you can use this to get a first-approximation posterior distribution over the -1, +1 values. So instead of assigning them hard values according to their sign, you can use the gaussian pdf to probabilistically assign them -1 or 1 based on the bad sigma estimate. Then you can re-estimate the posterior expected sigma based on these fuzzy -1, +1 re-estimates. You can iterate between re-estimating sigma and re-estimating the soft -1, +1 distribution.
Technically this
This is probably an application of the expectation-maximization principle, so this algorithm will converge to a local maximum likelihood estimate of sigma. Or if it is Furthermore the likelihood landscape of sigma seems plain enough to not have misleading local optima, then you so this should probably use also converge to the real expectation-maximization algorithm instead global maximum likelihood (and therefore statistically consistent) estimate of the iterative procedure I'm describingsigma.
