Learning about a hidden distribution behind Bernoulli trials Consider a process where a 'probability' $p$ is drawn from a fixed but unknown distribution $D$ supported on $[0,1]$ and then $y$ is drawn from Bernoulli distribution with parameter $p$. If $y_1,y_2..y_n$ are drawn by this process independently, then a reasonable estimator for $E_D[p]$ is $\sum_{}y_i/n$. Is there a distribution-free approach to estimate $E_D[f(p)]$ for any generic function $f$ ?

As $E[y^m]=E[y]$ for $m>0$, I doubt that $y$'s can tell us anything about $D$ apart from its expectation. What if we also know that $p_1\le p_2 \le . . \le p_n$?
 A: $\newcommand{\E}{\operatorname{E}}$This is nowhere close to a solution but it connects higher moments of $p$ to your data. Maybe it'll help if $f$ is a polynomial?
I'm assuming that $y_1,\dots,y_n\mid p \stackrel{\text{iid}}\sim \text{Bern}(p)$ with $n$ non-random. Let $S = \sum_{i=1}^n y_i$ so $S\mid p\sim\text{Bin}(n,p)$.
Let $G(z\mid p) = \E(z^S\mid p)$ be the conditional probability generating function of $S$. We have
$$
G(z\mid p) = \sum_{s=0}^n {n\choose s} (pz)^s(1-p)^{n-s} = (1-p + pz)^n
$$
by the binomial theorem in reverse. 
It is a standard result that
$$
G^{(k)}(1\mid p) = \E\left[\frac{S!}{(S-k)!}\mid p\right]
$$
(i.e. derivatives of $G$ give factorial moments of $S$) and I can see that
$$
G^{(k)}(1\mid p) = \frac{n!}{(n-k)!}p^k
$$
for $k \leq n$. This means that the marginal PGF of $S$ is
$$
G(z) = \E(z^S) = \E_D(G(z\mid p)) = \E_D((1-p+pz)^n).
$$
Differentiating here, and assuming I can exchange differentiation and integration, I have
$$
G^{(k)}(1) = \frac{n!}{(n-k)!}\E_D(p^n)
$$
so e.g.
$$
\E(S) = G'(1) = n \E_D(p)
$$
and
$$
\E(S(S-1)) = G''(1) = n(n-1)\E_D(p^2).
$$
Unfortunately it seems we only have one observation of $S$ so this might not be actually helpful, but it at least connects $S$ to higher moments of $p$.
If $n$ is sufficiently large we could partition the sample into $m$ groups of $k$ (so $n = mk$) and then take $S_1,\dots,S_m \mid p \stackrel{\text{iid}}\sim \text{Bin}(k, p)$ which would allow for some law of large numbers type estimators. 
A: This is similar to a Bernoulli factory problem. However, your problem has a solution only if f(p), in the interval [0, 1]—

*

*Is continuous, and

*is either constant on the interval [0, 1] and polynomially bounded away from 0 and 1

(Keane and O'Brien 1994).
Functions that meet these conditions are called factory functions, and only factory functions admit unbiased estimators when given an unlimited supply of Bernoulli(p) random variables (Łatuszyński et al. 2009/2011).
Merely being differentiable and concave is not enough for f to have an unbiased estimator. (Although being concave may simplify matters in designing an appropriate algorithm, depending on the function f in question, to Jensen's inequality.)
REFERENCES:

*

*Keane, M. S., and O'Brien, G. L., "A Bernoulli factory", ACM Transactions on Modeling and Computer Simulation 4(2), 1994.

*Łatuszyński, K., Kosmidis, I., Papaspiliopoulos, O., Roberts, G.O., "Simulating events of unknown probabilities via reverse time martingales", arXiv:0907.4018v2 [stat.CO], 2009/2011.

A: Take $f = -x^2 + x$. It is concave and differentiable, as specified. For any $\alpha \in [0,1]$, define $D_{\alpha} := \alpha\delta(p-\frac14) + (1-\alpha) \delta(p-\frac34)$, where $\delta$ is the Dirac delta. If we choose $p \sim D_{\alpha}$ we always get $f(p) = \frac{1}{4} + \frac{1}{16} = \frac{3}{16}$, for any $\alpha \in [0,1]$. Thus $E_{D_{\alpha}}[f(p)] = \frac{3}{16}, \forall \alpha \in [0,1]$.
Let's assume we can observe $n$ of your trials. We seek an estimator $E: \{0,1\}^n \to \mathbb{R}$. Whatever $E$ is, it will have to give similar answers for sequences distributed $B(n, \frac14)$ and $B(n, \frac34)$ (corresponding to $\alpha=1$ and $\alpha=0$, respectively). 
It seems to me that for that to hold, you will have to tailor your estimator to this given problem instance. I think you cannot find a general estimator (not without making further assumptions).
