1-28-2016 also look at MOM estimator that is based on E [ 1/X ] fx = t / x^2 (* x>t *) Integrate[ fx, {x, t, Infinity}, Assumptions -> t > 0 ] Integrate[ x fx, {x, t, Infinity}, Assumptions -> t > 0 ] Integrate[ (1/x) fx, {x, t, Infinity}, Assumptions -> t > 0 ] (* E[1/X] = 1/(2 theta) ==> theta.mom = 1 / (2 value) *) Fx = Integrate[ fx, {x,t,s}, Assumptions -> {t>0,s>t} ] Solve[ Fx==u, s ] my.sim = function(n=100,t=1,seed=129) { set.seed(seed) u = runif(n); xx = seq(0,50,,201); yy=rep(0,length(xx)) yy[xx>t] = t / xx[xx>t]^2 x = t / (1-u); hist(x[x<50],seq(0,50,1),col=2,xlim=c(0,50),prob=T) lines(xx,yy,col=4,lwd=3) mom = 1 / (2 * mean(1/x)); print(mom); browser() } my.sim() # MOM is 3 Xb Q: how often is x_(n) > 3 Xb ? test.sim = function(n=100,t=1,seed=129,nrep=1000) { set.seed(seed) ans=matrix(0,nrep,2) for(i in 1:n) { if(i%%50==0) cat(".") u = runif(n); x = t / (1-u) ans[i,1] = 3*mean(x); ans[i,2] = max(x) } }