diff --git a/+estimator/penalized_mle_inverse_covariance.m b/+estimator/penalized_mle_inverse_covariance.m index f14a522..1730d96 100644 --- a/+estimator/penalized_mle_inverse_covariance.m +++ b/+estimator/penalized_mle_inverse_covariance.m @@ -64,26 +64,18 @@ warning('Only lasso_mle supported') end - if(~isempty(options.W)) - warning('options.refit is turned on. Overwriting options.W'); - end - if(ndims(ThetaHat0)==2) - options.W = options.lambda_max*100.*(ThetaHat0==0) + ... - 1.0.*(ThetaHat0~=0); - else - warning(['Need to perform model selection first. ...' ... - ' Cannot refit to a path of estimates']); - end - end + [SigHat ThetaHat options] = constrain_mle(S,ThetaHat0,options); - switch options.estimator + else - case 'lasso_mle' - [SigHat ThetaHat] = lasso_mle(S,options); - otherwise - warning('Only lasso_mle supported') + switch options.estimator + + case 'lasso_mle' + [SigHat ThetaHat] = lasso_mle(S,options); + otherwise + warning('Only lasso_mle supported') + end end - end @@ -108,9 +100,28 @@ end -function [SigHat ThetaHat] = constrain_mle(S,options) +function [SigHat ThetaHat options] = constrain_mle(S,ThetaHat0,options) + + if(~isempty(options.W)) + warning('options.refit is turned on. Overwriting options.W'); + end + if(ndims(ThetaHat0)==2) + options.W = options.lambda_max*100.*(ThetaHat0==0) + ... + 1.0.*(ThetaHat0~=0); + else + warning(['Need to perform model selection first. ...' ... + ' Cannot refit to a path of estimates']); + end + + switch options.estimator + + case 'lasso_mle' + [SigHat ThetaHat] = lasso_mle(S,options); + otherwise + warning('Only lasso_mle supported') + end end diff --git a/@GGM/adaptiveWeights.m b/@GGM/adaptiveWeights.m index 093b864..601efa5 100644 --- a/@GGM/adaptiveWeights.m +++ b/@GGM/adaptiveWeights.m @@ -3,40 +3,44 @@ narginchk(1,2); p = length(Theta); - if(nargin==3) + if(nargin>=2) Lambda = varargin{1}; else Lambda = .01*ones(p,p); end nonzeroW = 1.*(abs(Theta)>1e-5); - A = diag(diag(Theta)) - Theta; + laplacian = @(A)(eye(length(A))- ... + diag(1./sqrt(diag(A)))*A*diag(1./sqrt(diag(A)))); + A = laplacian(Theta); t_start = cputime; - if(exist('mexeig')) - [V,D] = mexeig(A); - else - [V,D] = eig(A); - D = diag(D); - end - expA = V*diag(exp(D))*inv(V); + % if(exist('mexeig')) + % [V,D] = mexeig(A); + % else + % [V,D] = eig(A); + % D = diag(D); + % end + % expA = V*diag(exp(D))*inv(V); + expA = expm(A); timing = cputime-t_start; disp(sprintf('expm took %.4f sec',timing)); - % % Convert expA into a distance metric - % diagA = diag(expA); - % % A_i,i - A_i,j - % Gpos = bsxfun(@minus,diagA,expA); - % % A_j,j - A_j,i - % Gneg = bsxfun(@minus,diagA',expA'); - % distA = (Gpos + Gneg)/2; - % sqrtD = 1./(diagA); - % distA = scale_cols(distA,sqrtD); - % distA = scale_rows(distA,sqrtD); + % % Convert expA into a distance metric + diagA = diag(expA); + % A_i,i - A_i,j + Gpos = bsxfun(@minus,diagA,expA); + % A_j,j - A_j,i + Gneg = bsxfun(@minus,diagA',expA'); + distA = (Gpos + Gneg)/2; + sqrtD = 1./(diagA); + distA = scale_cols(distA,sqrtD); + distA = scale_rows(distA,sqrtD); + W = distA; + + %distA = 1./abs(expA); + %W = ((nonzeroW==1) + (nonzeroW==0).*(distA)).*Lambda; - % distA = 1./abs(expA); - % W = ((nonzeroW==1) + (nonzeroW==0).*(distA)).*Lambda; - - W = (nonzeroW==0).*Lambda; + % W = (nonzeroW==0).*Lambda; W = setdiagLS(W,0); end \ No newline at end of file diff --git a/@GGM/stableEstimate.m b/@GGM/stableEstimate.m index 6e0be72..fedd908 100644 --- a/@GGM/stableEstimate.m +++ b/@GGM/stableEstimate.m @@ -62,7 +62,7 @@ end tmpGGM = GGM(self.Data(:,:,1),1,0); [tmpGGM mleresults] = tmpGGM.estimate(); - mleresults.sparsemle.sparsity' + mleresults.sparsemle.sparsity'; self.ThetaPath = tmpGGM.ThetaPath; self.Theta = tmpGGM.ThetaPath(:,:,best_lambda); warning('Change to take in alternative grid of lambdas') diff --git a/test/test_mle_inverse_covariance.m b/test/test_penalized_mle_inverse_covariance.m similarity index 98% rename from test/test_mle_inverse_covariance.m rename to test/test_penalized_mle_inverse_covariance.m index 2ba92d1..d05f8c7 100644 --- a/test/test_mle_inverse_covariance.m +++ b/test/test_penalized_mle_inverse_covariance.m @@ -27,6 +27,7 @@ options.nlambdas = 1; options.Lambda = .1; + options.refit = true; results = estimator.fit(Data,options);