function sd=sder(beta,y,x,c);
%ll=sum(-(1-c).*(x*beta)-y.*exp(-x*beta));
[n,k]=size(x);
yexb=y.*exp(-x*beta);
sd=zeros(k,k);
for j1=1:k,
    for j2=1:k,
        sd(j1,j2)=sum(x(:,j1).*x(:,j2).*yexb);
        end,
    end