Seeking a continuous, parametric, bimodal sampling distribution for proportions I am seeking a parametric probability model whose pdf has the following characteristics: (1) it is supported on a variate axis that is bounded between 0 and 1; (2) it is continuous; and (3) it is capable of having two (or more) modes that are not necessarily located at 0 or 1 (e.g., in contrast to the beta distribution whose a and b are both <1).  I am thinking of some sort expansion on the beta or Kumaraswamy distribution analogous to the bi-Weibull expansion on the Weibull distribution.  My intention is to use such a model as a sampling distribution for sample proportions that vary continuously.  My simulated data tend toward the low side of the variate range, typically ranging between 0 and 0.01, with a center of mass that in the neighborhood of 0.001.  See the illustration below for examples.  Has anyone run across such an obscure distribution that meets this challenging combination of criteria?

 A: Why not consider a finite mixture distribution?
In this case, probably one with two components, though some of those plots might just be better approximated by three-component mixtures. 
Each component would be something suitable for a continuous proportion.
A typical example would be to use two different beta distributions.
So if the observed proportion is 
$$Y=wX_1+(1-w)X_2\,,$$ 
where the components are say
$$X_i\sim\text{Beta}(\alpha_i,\beta_i),\quad i=1,2$$ 
then the density is
$$f_Y(y) = wf_{X_1}(y)+(1-w)f_{X_2}(y),\quad 0\leq y\leq 1\,.$$
Here's an example:

The parameters of this example:
$\quad X_1\sim \text{Beta}(1,2500)$
$\quad X_2\sim \text{Beta}(4,3500)$
$\quad w=0.35$
The distribution of $Y$ in a two component mixture of betas has 5 parameters, and is very flexible.
A: I like the finite mixture approach, a lot.  Such an approach didn't occur to me at all, because I am not very familiar with finite mixture models.  Instead, I was envisioning a multi-component expansion on a beta or Kumaraswamy distribution, analogous to the multi-component models I am familiar with from survival analysis in demography, like the bi-Weibull or Siler distributions.  I decided to see if I had it in me to work with the Kumaraswamy distribution, because it has an analytically tractable CDF (unlike the beta distribution):



The challenge here, though, is to fit each component's a and b so that the pdf looks right.  Just tinkering with it, I came up with

Based on the R code
d.multi.kumar=function(y,a,b){
    k=length(a)

    p.per.k=rep(0,k)
    for (i in 1:k){
        p.per.k[i]=(1-y^a[i])^b[i]
    }
    S.multi.kumar=prod(p.per.k)

    h.per.k=rep(0,k)
    for (i in 1:k){
        h.per.k[i]=(a[i]*b[i]*y^(a[i]-1))/(1-y^a[i])
    }
    h.multi.kumar=sum(h.per.k)
    return(S.multi.kumar*h.multi.kumar)
}

a=c(1.01,4.2)
b=c(150,1000000000)
y=seq(0,0.2,0.000001)

pdf.multi.kumar=sapply(X=y,FUN=d.multi.kumar,a,b)
plot(y,pdf.multi.kumar,"h",xlim=c(0,0.1),xlab="y",ylab="f(y)",col="light blue")
    abline(h=0)

I like the visual effect of the finite mixture model more.  Wrestling sharp and horizontally separated peaks out of a multi-Kumaraswamy distribution seems quite difficult.  In terms of first principles, I don't have any great guidance, so I am not sure that a mixture model has any more theoretical warrant than the other approach, but whatever works ...
