Density estimation using (different) order statistics I need to estimate a univariate distribution $F$ as flexibly as possible. However, I do not observe draws from $F$ directly. Each observation $x_i$ is the minimum of $a_i$ draws from $F$, where $a_i$ is a known number from 2 to 5 that varies across observations.
I know how to solve this problem with $a_i=a$ fixed across observations. Kernel density estimation could find $F_{(1:a)}$, from which $F$ follows by
$F_{(1:a)}(x)^{1/a}=F(x)$
So the question is how $F$ can be estimated specifically when $a_i$ varies.
 A: One option is to estimate the CDF by maximum likelihood. The output is a step function which can be kernel smoothed.
Estimation
For every $x$ choose $F(x)$ to maximize the log-likelihood
$\sum_a A_{ax} \log((1-F(x))^a) + B_{ax} \log(1-(1-F(x))^a)$
where for observations with $a_i=a$ $A_{ax}$ is the number of such observations above $x$ and $B_{ax}$ is the number below x.
The first order condition is
$0=\sum_a A_{ax} a (1-F(x))^{-1} - B_{ax} a \frac{(1-F(x))^{a-1}}{1-(1-F(x))^a}$
which can be solved numerically (I used Julia's fzero). The solutions $F(x)$ form a step function akin to an empirical CDF. Optionally, follow up with a kernel smoother to to recover a smooth estimate of $F$.
Results
Assume $F \sim$ Weibull(2,2) for testing. With 20,000 observations the maximum likelihood estimate is close to the true distribution.

Results are tolerable with 200 observations.

Julia code
using Distributions
using Roots
using DataFrames
using Plots
srand(1111)

n = 200
a = rand(2:5,n)

dist=Weibull(2,2)

all_draws = rand(dist,sum(a))

x = fill(NaN,n)
k = 1
for i=1:n
  draws = all_draws[k:(k+a[i]-1)]
  k += a[i]
  x[i] = minimum(draws)
end

x=reshape(x,length(x),1)
df = convert(DataFrame,x)
rename!(df, :x1, :x)
df[:a] = a
sort!(df, cols=:x)

above = collect(countmap(df[:a]))
sort!(above, by=x->x[1])
above = vcat(map(t -> [t[2]], values(above))...)
below = convert(Array{Int,1},zeros(size(above)))

function foc(one_min_F,above,below)
  ret=0
  for a=2:5
    ret += above[a-1]*a/one_min_F
    ret -= below[a-1]*a*one_min_F^(a-1)/(1-one_min_F^a)
   end
  return ret
end

one_min_F=fill(NaN,n-1)
for i=1:n-1
  tmp = df[i,:a]-1
  above[tmp] -= 1
  below[tmp] += 1
  one_min_F[i] = fzero(x->foc(x,above,below),1e-30,1-1e-30)
end

F=1.-one_min_F

scatter(df[1:end-1,:x],F,label="estimated CDF")
plot!(df[1:end-1,:x],t->cdf(dist,t),label="theoretical CDF")

