在 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)

编辑

一种直接（非迭代）解决方案，与上面给出的解决方案兼容。这个想法是从价值观开始$Z$永远不会隐藏的（低于$min(Y)$)，并根据比例估计它们的概率。之后，这两个值和$min(Y)$可以从问题中删除。因此，问题变得越来越小。

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)
}