function Cvec = DualDiscreteImperfect(w,ghat,sigma2w,sigma2e,p,lambda,R,Eta)
Eta = Eta(:);
p = p(:);
R = R(:);
Lp = length(p);
gbar = sigma2e/sigma2w;
Etavec = repmat(Eta(2:end)',1,Lp-1);
Rvec =  repmat(R',1,Lp-1);
Pvec = reshape(repmat(p(2:end),1,Lp-1)',1,(Lp-1)*(Lp-1));
Pvecall = reshape(repmat(p(2:end),1,Lp)',1,(Lp-1)*(Lp));
m = Etavec./Pvec;
% CDF = 1-marcumq(sqrt(2*ghat/gbar)*ones(1,length(m)),sqrt(2*m/gbar));
CDF = zeros(size(m));
for i = 1 : length(m)
  a = sqrt(2*ghat/gbar)
  b = 
  CDF(i) = 1-marcumq(sqrt(2*ghat/gbar),sqrt(2*m/gbar));
end;
CDFvec = reshape([reshape(CDF,Lp-1,Lp-1);ones(1,Lp-1)],1,Lp*(Lp-1));
CDFvec2 = [0,CDFvec(1:end-1)];
CDFvec2(1:Lp:end)=0;
Cmat = reshape(w*Rvec.*(CDFvec-CDFvec2)-lambda*Pvecall,Lp,Lp-1);
Cvec = [0,sum(Cmat)];

