From 58dc9557d981978a6e6f974c1182a181a0a03321 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 30 May 2011 14:36:48 +0200 Subject: [PATCH] Fixes around point-estimate of sensitivity and identification strength --- matlab/identification_analysis.m | 38 +++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 11 deletions(-) diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m index e675149ba..61838f348 100644 --- a/matlab/identification_analysis.m +++ b/matlab/identification_analysis.m @@ -172,9 +172,14 @@ if info(1)==0, indok = find(max(ide_hess.indno,[],1)==0); cparam(indok,indok) = inv(AHess(indok,indok)); normaliz(indok) = sqrt(diag(cparam(indok,indok)))'; - cmm = siJ*((AHess)\siJ'); + cmm = NaN(size(siJ,1),size(siJ,1)); + ind1=find(ide_hess.ind0); + cmm = siJ(:,ind1)*((AHess(ind1,ind1))\siJ(:,ind1)'); + chh = siH(:,ind1)*((AHess(ind1,ind1))\siH(:,ind1)'); + ind1=ind1(ind1>offset); + clre = siLRE(:,ind1)*((AHess(ind1,ind1))\siLRE(:,ind1)'); rhoM=sqrt(1./diag(inv(tildaM(indok,indok)))); - deltaM = deltaM.*params'; +% deltaM = deltaM.*abs(params'); flag_score=1; catch, replic = max([replic, length(indJJ)*3]); @@ -193,33 +198,44 @@ if info(1)==0, indok = find(max(ide_hess.indno,[],1)==0); % rhoM=sqrt(1-1./diag(inv(tildaM))); % rhoM=(1-1./diag(inv(tildaM))); + ind1=find(ide_hess.ind0); + chh = siH(:,ind1)*((MIM(ind1,ind1))\siH(:,ind1)'); + ind1=ind1(ind1>offset); + clre = siLRE(:,ind1)*((MIM(ind1,ind1))\siLRE(:,ind1)'); rhoM(indok)=sqrt(1./diag(inv(tildaM(indok,indok)))); normaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))'; - deltaM = deltaM.*params'; +% deltaM = deltaM.*abs(params'); flag_score=0; end ide_strength_J(indok) = (1./(normaliz(indok)'./abs(params(indok)'))); ide_strength_J_prior(indok) = (1./(normaliz(indok)'./normaliz1(indok)')); ide_strength_J(params==0)=ide_strength_J_prior(params==0); + deltaM_prior = deltaM.*abs(normaliz1'); + deltaM = deltaM.*abs(params'); + deltaM(params==0)=deltaM_prior(params==0); quant = siJ./repmat(sqrt(diag(cmm)),1,nparam); siJnorm = vnorm(quant).*normaliz1; % siJnorm = vnorm(siJ(inok,:)).*normaliz; quant=[]; - inok = find((abs(TAU)<1.e-8)); - isok = find((abs(TAU)>=1.e-8)); - quant(isok,:) = siH(isok,:)./repmat(TAU(isok,1),1,nparam); - quant(inok,:) = siH(inok,:)./repmat(mean(abs(TAU)),length(inok),nparam); +% inok = find((abs(TAU)<1.e-8)); +% isok = find((abs(TAU)>=1.e-8)); +% quant(isok,:) = siH(isok,:)./repmat(TAU(isok,1),1,nparam); +% quant(inok,:) = siH(inok,:)./repmat(mean(abs(TAU)),length(inok),nparam); + quant = siH./repmat(sqrt(diag(chh)),1,nparam); siHnorm = vnorm(quant).*normaliz1; % siHnorm = vnorm(siH./repmat(TAU,1,nparam)).*normaliz; quant=[]; - inok = find((abs(LRE)<1.e-8)); - isok = find((abs(LRE)>=1.e-8)); - quant(isok,:) = siLRE(isok,:)./repmat(LRE(isok,1),1,np); - quant(inok,:) = siLRE(inok,:)./repmat(mean(abs(LRE)),length(inok),np); +% inok = find((abs(LRE)<1.e-8)); +% isok = find((abs(LRE)>=1.e-8)); +% quant(isok,:) = siLRE(isok,:)./repmat(LRE(isok,1),1,np); +% quant(inok,:) = siLRE(inok,:)./repmat(mean(abs(LRE)),length(inok),np); + quant = siLRE./repmat(sqrt(diag(clre)),1,np); siLREnorm = vnorm(quant).*normaliz1(offset+1:end); % siLREnorm = vnorm(siLRE./repmat(LRE,1,nparam-offset)).*normaliz(offset+1:end); ide_hess.ide_strength_J=ide_strength_J; ide_hess.ide_strength_J_prior=ide_strength_J_prior; + ide_hess.deltaM=deltaM; + ide_hess.deltaM_prior=deltaM_prior; ide_moments.siJnorm=siJnorm; ide_model.siHnorm=siHnorm; ide_lre.siLREnorm=siLREnorm;