How do you calculate the uncertainty on a reconstructed probability distribution? This question may be very basic, but I have looked through some resources on Monte Carlo statistical analysis as well as other methods to reconstruct a distribution, but have not found an answer suitable so far. I describe my ansatz solution at the bottom, but am not satisfied with it for the reasons I state.
I have a continuous probability distribution $\mathcal{P}(x)$, for $x\in [0, 1]$, where $\mathcal{P}(x)$ is the probability of picking a $+1$ vs a $0$. 
I then randomly pick $N$ locations $x_i\in[0, 1]$, and get out $+1$ or $0$ at each location, at a rate proportional to the probability at each $x_i$. As in, if I picked $x_i=y$ for all $i$, and $N\rightarrow \infty$, then the number of $+1$ vs $0$ would approach $\mathcal{P}(y)$.
Given this, for a finite $N$, how do I best estimate $\mathcal{P}(x)$ for ALL $x$, and what is the uncertainty on that estimate at each point $x$?
My first idea was to place a Gaussian at each $x_i$ that received a $+1$, then divide by Gaussians placed at $+1$ and $0$. This seems intuitive, but the standard deviation is a parameter I don't think is relevant. Furthermore, I am uncertain of how to construct the error, $\delta\mathcal{P}(x)$ in my reconstruction. Lastly, I have not grounded this technique in any statistical formalism, which is my ultimate goal.
Thank you in advance.
Update 1
Another idea I had was to randomly drop groups of the measurements (maybe 50% at a time) at each reconstruction, and use the standard deviation about the mean reconstruction as the error.
Update 2
Thinking more about this, it seems to me that maybe continuous is not enough for this problem to be well defined. Therefore, let's say that the probability distribution is known to be Lipschitz continuous with Lipschitz constant $L$, such that:
$$\left \lvert\frac{\mathcal{P}(x_2)-\mathcal{P}(x_1)}{x_2-x_1}\right \rvert \leq L\ \forall\ x_1, x_2 \in [0, 1]$$
For the Gaussians, perhaps the standard deviation could be a function of $N$ (giving the density of samples), and $L$. 
 A: I think you figured it out yourself, but I can point you to a more general framework.
Basically, you can treat this as a nonparametric regression problem where you need to regress your outcomes +1 or 0 on $x$, so any of the several nonparametric regression techniques will do. Here are a few straightforward solutions in increasing order of sophistication:
Binning: Just group your $x_i$'s into discrete uniformly spaced bins and count the total fraction of +1 outcomes in each bin.
k-nearest neighbors: Choose a set of target points that are evenly spaced in $x$. For each point in that set, find the k nearest neighbors from among the set of all $x_i$'s and compute the fraction of +1 outcomes corresponding to those neighbors.
kernel regression: This is what you had in mind. Pick a kernel of your choice (Gaussian will do) and place it at each $x_i$ that produced a +1 and add them up. Divide this by the sum of kernels placed at each of your $x_i$'s regardless of outcome.
Local linear regression: This is a generealization of the above technique where you assume that $P(x)$ is locally linear so there exist some $a, b$ for each $x_i$ such that $P(x_i) = a + bx_i$. If your $x_i$'s are very non-uniformly spaced (if there are empty pockets with no samples), this will do better than naive kernel smoothing.
You'll find more such techniques here: http://www.stat.cmu.edu/~larry/=sml/nonpar.pdf
For each of these techniques, you can compute uncertainty of the estimator by bootstrapping. Rather than dropping 50% of your observations, you resample your data with replacement and recompute the estimator several times.
I cannot comment on the choice of k or bandwidth of the Gaussian as that would depend on the ensemble of your observations. One technique you can use if your $N$ is large enough, is evaluate many different choices of bandwidth using a cross-validation set. And pick the bandwidth with the least cross-validation error.
