@@ -13,15 +13,23 @@ import ..Diagnostics: mahalanobisSquaredMatrix
1313
1414
1515function enlargesubset (initialsubset, data:: DataFrame , dataMatrix:: AbstractMatrix , h:: Int )
16- n, _ = size (dataMatrix)
16+ n, p = size (dataMatrix)
17+
1718 basicsubset = copy (initialsubset)
19+
20+ 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)
25+
1826 while length (basicsubset) < h
19- meanvector = applyColumns (mean, data[basicsubset, :])
20- covmatrix = cov (dataMatrix[basicsubset, :])
21- md2mat =
27+ meanvector . = applyColumns (mean, data[basicsubset, :])
28+ covmatrix . = cov (dataMatrix[basicsubset, :])
29+ md2mat . =
2230 mahalanobisSquaredMatrix (data, meanvector = meanvector, covmatrix = covmatrix)
23- md2 = diag (md2mat)
24- md2sortedindex = sortperm (md2)
31+ md2 . = diag (md2mat)
32+ md2sortedindex . = sortperm (md2)
2533 basicsubset = md2sortedindex[1 : (length (basicsubset)+ 1 )]
2634 end
2735 return basicsubset
@@ -38,20 +46,30 @@ function robcov(data::DataFrame; alpha = 0.01, estimator = :mve)
3846 indices = collect (1 : n)
3947 k = p + 1
4048 mingoal = Inf
41- bestinitialsubset = []
42- besthsubset = []
49+
4350 maxiter = minimum ([p * 500 , 3000 ])
44- initialsubset = []
45- hsubset = []
51+
52+ initialsubset = Array {Int} (undef, k)
53+ bestinitialsubset = Array {Int} (undef, k)
54+
55+ hsubset = Array {Int} (undef, h)
56+ besthsubset = Array {Int} (undef, h)
57+
58+ covmatrix = Matrix {Float64} (undef, p, p)
59+ meanvector = Array {Float64} (undef, p)
60+ md2mat = Matrix {Float64} (undef, n, n)
61+
62+ md2 = Array {Float64} (undef, n)
63+
4664 for iter = 1 : maxiter
4765 goal = Inf
4866 try
49- initialsubset = sample (indices, k, replace = false )
50- hsubset = enlargesubset (initialsubset, data, dataMatrix, h)
51- covmatrix = cov (dataMatrix[hsubset, :])
67+ initialsubset . = sample (indices, k, replace = false )
68+ hsubset . = enlargesubset (initialsubset, data, dataMatrix, h)
69+ covmatrix . = cov (dataMatrix[hsubset, :])
5270 if estimator == :mve
53- meanvector = applyColumns (mean, data[hsubset, :])
54- md2mat = mahalanobisSquaredMatrix (
71+ meanvector . = applyColumns (mean, data[hsubset, :])
72+ md2mat . = mahalanobisSquaredMatrix (
5573 data,
5674 meanvector = meanvector,
5775 covmatrix = covmatrix,
@@ -66,25 +84,25 @@ function robcov(data::DataFrame; alpha = 0.01, estimator = :mve)
6684 end
6785 if goal < mingoal
6886 mingoal = goal
69- bestinitialsubset = initialsubset
70- besthsubset = hsubset
87+ bestinitialsubset . = initialsubset
88+ besthsubset . = hsubset
7189 end
7290 end
73- meanvector = applyColumns (mean, data[besthsubset, :])
74- covariancematrix = cov (dataMatrix[besthsubset, :])
75- md2 = diag (
91+ meanvector . = applyColumns (mean, data[besthsubset, :])
92+ covmatrix . = cov (dataMatrix[besthsubset, :])
93+ md2 . = diag (
7694 mahalanobisSquaredMatrix (
7795 data,
7896 meanvector = meanvector,
79- covmatrix = covariancematrix ,
97+ covmatrix = covmatrix ,
8098 ),
8199 )
82100 outlierset = filter (x -> md2[x] > chisqcrit, 1 : n)
83101 result = Dict {String, Any} ()
84102 result[" goal" ] = mingoal
85103 result[" best.subset" ] = sort (besthsubset)
86104 result[" robust.location" ] = meanvector
87- result[" robust.covariance" ] = covariancematrix
105+ result[" robust.covariance" ] = covmatrix
88106 result[" squared.mahalanobis" ] = md2
89107 result[" chisq.crit" ] = chisqcrit
90108 result[" alpha" ] = alpha
0 commit comments