function james_stein_risk(nsim,p)
theta=.1:.1:4; a=.1:.2:2*(p-2); ntheta=length(theta); na=length(a);
x=randn(nsim,p);
risk=zeros(ntheta,na);
for i=1:ntheta
    for j=1:na
        z=theta(i)+x;
        m=max(0,(1-a(j)./sum(z'.^2)));
        m=repmat(m',1,p);
        d=z.*m;
        risk(i,j)=sum(mean((d-theta(i)).^2));
    end
end
mesh(repmat(theta',1,na),repmat(a,ntheta,1),risk)
xlabel('theta')
ylabel('a')
zlabel('risk')