Estimating distribution from censored data

Question

$X$ is a positive variable with known support (assume discrete support, if that simplifies solution).

$Y$ is another variable with the same support.

$X$ and $Y$ are independent.

$Z$ is equal to $X$ if $X < Y$, and equal to $0$ otherwise.

$Y$ and $Z$ are observed, $X$ is not. How to estimate distribution of $X$?

I realize that there may be no single objectively optimal answer. Assume prior estimates, if needed for Bayesian inference.

Answer 1

In R:

estimate=function(y,z,u=1e-9){
  ys=sort(unique(y))
  # Inf signifies x's never observed (as they are higher than max y)
  zs=c(sort(unique(z))[-1],Inf)
  counts=xtabs(~z+y)
  observed=rbind(counts[-1,],rep(0,length(ys)))
  marginalHidden=counts[1,]
  m=sapply(seq(ys),function(i)zs>ys[i])
  d=rep(1/length(zs),length(zs))
  while(T){
    # allocate hidden data according to current parameters
    p=apply(m*d,2,function(v)v/sum(v))
    # can result in fractional counts
    hidden=sweep(p,2,marginalHidden,'*')
    total=observed+hidden
    d2=apply(total,1,sum)/sum(total)
    msd=mean((d2-d)^2)
    if(msd<u^2)
      break;
    d=d2
  }
  d
}

xSupport=c(3,5,7)
xDistribution=c(1/4,1/2,1/4)
x=sample(xSupport,1000,replace=T,prob=xDistribution)
ySupport=c(4,6)
yDistribution=c(1/2,1/2)
y=sample(ySupport,length(x),replace=T,prob=yDistribution)
z=ifelse(x<y,x,0)

estimate(y,z)
table(x)

---

Edit

A direct (non-iterative) solution, compatible with the one given above. The idea is to start with the values of $Z$ that are never hidden (lower than $min(Y)$), and estimate their probability from proportions. After that, both these values and $min(Y)$ can be removed from the problem. Thus, the problem becomes smaller and smaller.

estimate=function(y,z,u=1e-9){
  ys=sort(unique(y))
  # Inf signifies x's never observed (as they are higher than max y)
  z[z==0]=Inf
  zs=sort(unique(z))
  counts=xtabs(~z+y)
  s=c()
  r=1
  while(ncol(counts)>0){
    # zs < min(ys) are all observed, so can be estimated from counts
    mzi=which(zs<min(ys))
    ds=r*apply(counts[mzi,,drop=F],1,sum)/sum(counts)
    s=c(s,ds)

    # reduce probability remaining for the hidden cases
    r=r-sum(ds)

    # reduce the problem by removing the solved levels of zs, and the min(ys)
    counts=counts[-mzi,-1,drop=F]
    zs=zs[-mzi]
    ys=ys[-1]
  }
  c(s,r)
}