quant=[]; n=100; nsim=1000;
thet=2.3
xp=gaminv(.999,10,1);
x=gamrnd(10,1,n,nsim);
z=thet*x;  
z=thet*x;  w=(thet^10)*exp(-z*(1-1/thet));
d=z<xp;
v=[.001*.999 mean(var(d.*w))];