X = []; for i=1:100; X=[X;q4_8(1000,2.9)];   end
mean(X)
V=cov(X);
Z=ones(3,1);
V1=inv(V);
b=V1*Z/(Z'*V1*Z);
o=mean(X*b);
v=1/(Z'*V1*Z);  
X=[X X*b];
disp(['estimator means ' num2str(mean(X))])
disp('estimator variances ' )
var(X)