@@ -8,28 +8,26 @@ import LinearAlgebra: diag, det
88import Distributions: median, cov, mean, quantile, sample, Chisq
99
1010import .. Basis:
11- RegressionSetting, @extractRegressionSetting , designMatrix, responseVector, applyColumns
11+ RegressionSetting,
12+ @extractRegressionSetting , designMatrix, responseVector, applyColumns, applyColumns!
13+
1214import .. Diagnostics: mahalanobisSquaredMatrix
1315
1416
1517function enlargesubset (initialsubset, data:: AbstractMatrix , h:: Int )
16- n, p = size (data)
18+ p = size (data, 2 )
1719
1820 basicsubset = copy (initialsubset)
1921
2022 meanvector = Array {Float64} (undef, p)
21- covmatrix = Matrix {Float64} (undef, p, p)
22- md2mat = Matrix {Float64} (undef, n, n)
23- md2 = Array {Float64} (undef, n)
24- md2sortedindex = Array {Int} (undef, n)
2523
2624 while length (basicsubset) < h
27- meanvector . = applyColumns ( mean, data[basicsubset, :])
28- covmatrix . = cov (data[basicsubset, :])
29- md2mat . =
25+ applyColumns! (meanvector, mean, data[basicsubset, :])
26+ covmatrix = cov (data[basicsubset, :])
27+ md2mat =
3028 mahalanobisSquaredMatrix (data, meanvector = meanvector, covmatrix = covmatrix)
31- md2 . = diag (md2mat)
32- md2sortedindex . = sortperm (md2)
29+ md2 = diag (md2mat)
30+ md2sortedindex = sortperm (md2)
3331 basicsubset = md2sortedindex[1 : (length (basicsubset)+ 1 )]
3432 end
3533 return basicsubset
@@ -55,21 +53,22 @@ function robcov(data::Matrix; alpha = 0.01, estimator = :mve)
5553 hsubset = Array {Int} (undef, h)
5654 besthsubset = Array {Int} (undef, h)
5755
58- covmatrix = Matrix {Float64} (undef, p, p)
56+
5957 meanvector = Array {Float64} (undef, p)
60- md2mat = Matrix {Float64} (undef, n, n )
58+ fill! (meanvector, 0.0 )
6159
62- md2 = Array {Float64} (undef, n)
6360
6461 for iter = 1 : maxiter
6562 goal = Inf
63+
64+
6665 try
67- initialsubset . = sample (indices, k, replace = false )
68- hsubset . = enlargesubset (initialsubset, data, h)
69- covmatrix . = cov (data[hsubset, :])
66+ initialsubset = sample (indices, k, replace = false )
67+ hsubset = enlargesubset (initialsubset, data, h)
68+ covmatrix = cov (data[hsubset, :])
7069 if estimator == :mve
71- meanvector . = applyColumns ( mean, data[hsubset, :])
72- md2mat . = mahalanobisSquaredMatrix (
70+ applyColumns! (meanvector, mean, data[hsubset, :])
71+ md2mat = mahalanobisSquaredMatrix (
7372 data,
7473 meanvector = meanvector,
7574 covmatrix = covmatrix,
@@ -82,17 +81,20 @@ function robcov(data::Matrix; alpha = 0.01, estimator = :mve)
8281 catch e
8382 # Possibly singularity
8483 end
84+
85+
8586 if goal < mingoal
8687 mingoal = goal
87- bestinitialsubset . = initialsubset
88- besthsubset . = hsubset
88+ bestinitialsubset = initialsubset
89+ besthsubset = hsubset
8990 end
9091 end
9192
9293
93- meanvector .= applyColumns (mean, data[besthsubset, :])
94- covmatrix .= cov (data[besthsubset, :])
95- md2 .= diag (
94+
95+ applyColumns! (meanvector, mean, data[besthsubset, :])
96+ covmatrix = cov (data[besthsubset, :])
97+ md2 = diag (
9698 mahalanobisSquaredMatrix (
9799 data,
98100 meanvector = meanvector,
0 commit comments