# Estimating distribution from censored data

$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.

• I have an idea how to solve this using MCMC, but I would prefer something less computationally intensive. – Andris Birkmanis Dec 9 '13 at 23:10
• – Andris Birkmanis Dec 17 '13 at 16:53
• The basic idea is to initialize estimate of distribution from observed data, then iteratively: 1. sample censored data according to current estimate 2. pick distribution that maximizes likelihood given both observed and "guessed" censored data – Andris Birkmanis Dec 17 '13 at 23:56

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

• An explanation of what this code is doing would be welcome. – whuber Dec 18 '13 at 23:16
• I followed the idea mentioned in Parameter Estimation from Censored Samples using the Expectation-Maximization Algorithm, Section 3. I estimate parameters (in my case, the discrete distribution itself) from observed data, then iteratively "sample" (actually allocate proportionally) hidden data according to current parameters, then re-estimate parameters given both observed and "guessed" hidden data, until parameters stop changing significantly. Now, I would be glad to find a method without iterations - not sure if it's possible. – Andris Birkmanis Dec 18 '13 at 23:18
• One problem with current code is it allocates fractional counts of hidden data, which may produce unrealistic parameters. I leave this fix as an exercise :) – Andris Birkmanis Dec 18 '13 at 23:27