the way i usually hear of this done is via a beta-binomial mixture distribution. it's a bit advanced, but i'll outline the concept anyway. you're going to want to hunt for a beta-binomial package in R or something to actually evaluate probabilities.
given probability p, X follows a binomial distribution. the beta distribution is continuous and assigns probability to a variable on the unit interval (0,1), so it's a popular choice for modeling probabilities (and can be used to model the probability P in this case).
X|P=p ~ bin(n, p)
P ~ beta(a, b),
the density of X, after removing the condition on knowing p, is a tad complex:
f(x) = nCx * beta(a+x, b+n-x) / beta(a, b)
- nCx is a binomial coefficient
- beta(a, b) is a beta function, using the form beta(a, b) = gamma(a)gamma(b)/gamma(a + b)
- and gamma(a) is a gamma function
you can use maximum likelihood estimates of a and b (when estimating p), and then your density is simply a function of x. a and b carry all the information about the beta-modeled binomial probability p.
if you have limited exposure to math/stats, i'm aware this might look a bit daunting. but with some experience it's really not too bad, and like i said the applied work could be as easy as finding the correct package in R.