Monte Carlo with error on individual samples I'm performing a Monte Carlo integration where the individual samples have an error, and I'm wondering how to estimate the final error.
Some more detail: The integral E I'm after is estimated in the Monte Carlo as
$\langle E\ \rangle\approx\frac{\sum_ke_k\ F_k}{\sum_kF_k}\ ,\ \ \ \ \ \ $ (1) 
where the denominator comes from normalization. The individual $F_k$ have errors $\sigma_k$ of their own (actually they are computed with a secondary Monte Carlo), in the form of a standard deviation around the average. This error $\sigma_k$ depends on how many samples $m$ there are in the secondary MC. 
But it turns out that $\langle E\ \rangle$ "inherit" these errors in a nonsymmetric way. Specifically, taking a single one of the $F_k$ and exchanging $F_k\rightarrow F_k+\sigma_k$ or $F_k\rightarrow F_k-\sigma_k$, leading to a $\langle E\ \rangle_+$ and $\langle E\ \rangle_-$ respectively, I always find 
$\mathcal{D}=\frac{\big|\ \langle E\ \rangle_+-\ \langle E\ \rangle\ \big|}{\big|\ \langle E\ \rangle_--\ \langle E\ \rangle\ \big|}<1\ \ \ \ \ \ \ \ $ and $\ \ \ \ \ \ \ \ \langle E\ \rangle_-<\langle E\ \rangle_+\ .$
As one would expect $\mathcal{D}$ becomes smaller with $m$ (since the errors $\sigma_k$ get smaller) and with the number $n$ of samples in the main Monte Carlo (since both the numerator and denominator in (1) become bigger).
Now my question is how to estimate the errors on (1). Normally one would just use the standard deviation, which in this case could be found using blocking, jackknife or probably most stably bootstrap. But it seems like this won't work in this case, since there is a sort of systematic error involved (namely, that $\langle E\ \rangle$ will decrease with $m$).
To try to test this I tried running the main Monte Carlo several times for different kinds of parameters. And indeed I see that in general the outcome decreases with $m$, and more importantly, doing the same run (same parameters) a number of times I get a higher variation in the answers than one would expect from the standard deviation of the mean $\frac{\sigma}{\sqrt{n-1}}$ calculated from (1) using different methods including bootstrap.
I could of course use this variation of several calculated means as an error estimate but firstly this might not be robust either because of systematic error-like behaviour, and secondly this would mean I have to perform a lot of extra calculations just to get the error.
Any ideas on how to find the actual error on $\langle E\ \rangle$?
 A: You may want to check your sampling with regards to the asymmetry.  That suggests a bias is somehow being introduced into whatever Monte Carlo subprocesses are showing that asymmetric uncertainty, or that those subprocesses have yet to reach equilibrium.  
I've observed something similar when accidentally non-uniform sampling point picking on a sphere.  If I naively pick an $\mathsf { x<r }$ and then pick y via the remainder $\mathsf { y<\sqrt{r^2-x^2} }$ angles , I'll most commonly get angles at 90 degrees an arbitrary polar axis (a point on the surface and the center of the sphere).  Hence if I'm Boltzmann sampling based on $E\left(\theta+\delta\right)=E\left(\theta-\delta\right)$ my energy will be symmetric, but it's moot as if my central angle $\mathsf {\theta_0=114°} $ thus $\mathsf {\theta=110°} $ and $\mathsf {\theta=118°} $ will produce equal acceptance rates (per the Boltzmann sampling scheme) but the bigger angle will be picked less often (which can be mentally visualized as it traces a smaller circular cross section of the sphere versus the more central angle).  Hence I've introduced a subtle bias that's pushing my equilibrium angle to something other than $\mathsf { \theta_{avg}<\theta_0=114° }$.
I'm fairly confident mathematically your result means your asymmetric processes either are sampling at nonequilibrium or in biased fashion.  The question is whether that's intentional.
If it's not, you can fix it.  Your probability should now follow standard propogation of errors:
The denominator should be given by the propogation of errors formula for an average [source]:
$$\mathsf \langle \sum_kF_k \rangle=\sqrt{\sum_k{\left(\sigma_k-\overline{\sigma}\right)}^2}\,\,\,\,\,\,\,\,\,\,\,\,\,\,(1)$$
I'm going to guess that $\mathsf e_k$ is some output with numerical uncertainty of its own (due to the stochastic nature of the sampling) which goes to zero as you sample more states.  Let's call the relative uncertainty in this quantity as: 
$$\mathsf \langle \delta e_k \rangle$$
So the error in the numerator per the product rule applied via the chain rule to the average form should be soemthing like:
$$ \mathsf \langle \sum_ke_kF_k \rangle =\sqrt{\sum_k{\left(\langle e_kF_k \rangle-\overline{\langle e_kF_k \rangle}\right)}^2}\,\,\,\,\,\,\,\,\,\,\,\,\,\,(2)$$
where
$$\mathsf \langle e_kF_k \rangle=\sqrt{{\left(\frac{\langle \delta e_k \rangle}{e_k}\right)}^2+{\left(\frac{\sigma_k}{F_k}\right)}^2}\,\,\,\,\,\,\,\,\,\,\,\,\,\,(3)$$
Finally putting Eqn (1) and (2) all this together, you have via the division rule for the numerator and denominator:
$$\mathsf \langle E \rangle =\sqrt{{\left(\frac{\langle \sum_ke_kF_k \rangle}{\sum_ke_kF_k}\right)}^2+{\left(\frac{\langle \sum_kF_k \rangle}{\sum_kF_k}\right)}^2}$$
If indeed your bias is due to error or non-equilbrium and your intent is to sample without bias and equilbrium, adjust your Monte Carlo simulations to properly equilibrate and to remove and biases, then use the above derivation. (If this the case and you use the formula in a published work, please cite/acknowledge my contribution and email me a note on that at jason.is.coding@gfilter.com (replace "filter" with "mail"), or simply leave me a note here in the comments.
If your intent is indeed to sample with an unremoved bias (which violates detailed balance if you're using the Metropolis-Hastings sampling scheme) I suggest referring to Barlow's 2003 paper (arXiv:physics/0306138v1) on a means of propogating/quantifying uncertainty in an asymmetric sampling scheme.  The basic gist is you add together the asymmetric errors in both directions and either divide by $\mathsf 2 $ or $\mathsf { \sqrt{2\pi} }$, i.e.
$$\mathsf {\sigma=\frac{\sigma^{+}+\sigma^{-}}{2}}$$
The real key is the introduction of an asymmetry factor:
$$\mathcal{A}\mathsf{=\frac{\sigma^{+}-\sigma^{-}}{\sigma^{+}+\sigma^{-}}}$$
There's also a difference factor ($\mathsf\alpha$). 
To combine the asymmetric and nonasymmetric, I assume you would just follow their combining rules with $\mathcal{A}\mathsf{=1}$ and $\mathsf {\alpha=0}$.  It looks like that's fairly convoluted/nontrivial, though, so you probably want to be sure that's really really what youre intending to do.
But assuming it is necessary, there should be enough there to provide at least one analytical solution for you wrt to the individual numerically solved error distributions.
There's probably more proposals of how to do that if you search "asymmetric errors".  I also found this MIT paper on the topic who suggests a slightly different scheme.  I'm not an expert in particular on that field of math, but based on a quick search on scholarly search portals it appears to be a somewhat new field with no clear consensus.  My best guess is that's due to the fact that in most cases this kind of distribution, like my experience with the angle picking, represents an error and not true intent to introduce a bias (and not a bias from some nonuniform source you're measuring).
If you go that route, seems like the best option is to pick one scheme and stick with it.  The Barlow one seems pretty well explained although, again, the combining rules sound pretty challenging.
