# Estimating an underlying pdf from binomial trials

I'm afraid I'm not an expert in statistics, but I have a particular problem I'm interested in solving. I'm pretty sure this area already has a lot of literature, but I'm having difficulty finding something directly applicable to what I'm doing, so it would be great if somebody could nudge me in the right direction.

I'll give an example to illustrate the problem: Say I have a machine which produces biased coins. It has some underlying continuous probability density function which it is using to pick a number 0 < p < 1, then creates a coin which will come up heads with probability p.

My task is to estimate the function which the machine uses to generate these coins. I'm allowed to flip the coins, and I have access to a very large number of coins, but each one after a certain, random, amount of flips is taken away from me. The number of flips is not necessarily large.

How would I go about doing this? My initial thought was to just sum each resulting binomial distribution and divide by the total number of coins tested. But I'm pretty sure this doesn't give a good result.

I'm faintly aware of kernal density estimation, but I don't have enough expertise to know if/how it can be used for this kind of task, or what I should know in terms of tailoring it to this task.

• Hmm, interesting Q... What's the distribution of the number of flips? Is it independent of $p$ ? – onestop Apr 30 '12 at 13:51
• Yeah, for simplicity's sake, say the distriubtion of the number of flips is known to be independent of p. Other than that it's unknown, and we're not really interested in finding it out (unless there's some reason we need to in order to achieve the main task) – Ben Aaronson Apr 30 '12 at 23:19

This is just a simple idea and is not something I have seen in the literature. I will take away the randomness of the flips by conditioning on the observed number for each coin. Take the usual estimate for pi (i.e. number of heads divided by the number of flips) for the ith randomly selected coin. This set of estimates forms a histogram and one can then use a kernel density method to approximate the continuous curve. The difficulty with this approach is that it ignores the uncertainty in the estimate of p which depends on the number of flips and the true p. If ni is the number of flips for the ith coin and is large for each i, ignoring this uncertainty will not matter. I think it complicates things a little that each coin has a different variance associated with its estimate of p. Maybe this can be taken into account by using a variable width kernel.

• Please consider editing to use $\LaTeX$. It's easier on the eyes (and prettier). :) – cardinal May 2 '12 at 18:21
• Thank you for the response. Unfortunately this skips most of the key detail that I need, about how to implement a kernal density method for this specific problem, or how exactly to use a variable width kernal. I also forgot to mention (and have now editted in) that we can not rely on the number of flips being large for a given coin (though the number of coins is large) – Ben Aaronson May 3 '12 at 2:19
• I don't think that it is my job to do a tutorial on kernal density estimation as part of providing an answer to your question. Bernie Silverman would do a better job in his Density Estimation text. You could probably read about it in wikipedia or plenty of other books including one by David Scott. – Michael R. Chernick May 4 '12 at 1:59
• @Stereotomy, the density function in R gives a kernel density estimate if you pass it the list of data points - this answer suggests using the list of $\hat{p}$'s as the data points. – Macro May 4 '12 at 13:34
• Okay, thank you for the recommendation of Silverman and Scott. I guess what I was hoping for as an answer was any considerations for a general approach like kernal density estimation which would be specific to this problem. Binomial distributions for low n are only parially approximated by normal distributions, so would there be some particular kernal function which would be more suitible? Does the problem suggest a particular way of choosing a bandwidth? And so on. I don't want to be coached in the basics of the kernal density method, but I want to know how to adapt it to this specifically. – Ben Aaronson May 5 '12 at 1:03

Your problem can be generalized as follows:

1. The parameter $$p$$ is generated according to an unknown probability density.
2. You can measure $$p$$ only with a measurement error $$\sigma$$, which yields an estimator $$\hat{p}$$ for $$p$$ that has a known distribution (binomial distribution).

A simple approach is suggested by Glen_b in this answer to a similar question: estimate a global bandwidth $$h$$ first as if the measurements $$\hat{p}_1,\ldots,\hat{p}_N$$ were exact, and then increase the bandwidth to $$\sqrt{h^2+\sigma^2}$$. In your case $$\sigma_i^2=\hat{p}_i(1-\hat{p}_i)/n_i$$. For more sophisticated methods, see

Achilleos, Delaigle: Local bandwidth selectors for deconvolution kernel density estimation. Statistics and Computing 22,2 pp. 563–577, 2012