function [m,s,p]=LRexample(n)
m=[]; s=[]; p=.6:.03:1.2;
for psi=p; 
X=gamrnd(n,psi,1,30000);
m=[m mean(X.*exp(X*(1-1/psi)))];
s=[s sqrt(var(X.*exp(X*(1-1/psi))))];
end
m=m./(n*p.^n); s=s./(n*p.^n);
plot(p,s,'k')
xlabel('psi')
ylabel('standard deviation(T)')
figure
plot(p,m)
xlabel('psi')
ylabel('mean(T)')