Hi:

I am mostly a R user, but am currently developing a package that has the primary functions written in Julia. This package is for estimating very large (inverse) co-variance matrices. I have considerable speed gains compared to R, already, but as a novice I was looking for some input to make the Julia code even more efficient.

```
function col_means(x)
mean(x, 1)
end
function center_data(x)
x .- col_means(x)
end
function post_prob(x)
grt_zero = mean(repeat([0], inner = size(x, 1)) .> x)
ls_zero = 1 - grt_zero
1 - min(ls_zero, grt_zero)
end
function BB_est(x, sims, prob)
X = center_data(x)
n = size(X, 1)
p = size(X, 2)
prob = convert(Float64, prob)
offdiag = (p * (p - 1)) / 2
loop_max = convert(Int64, offdiag)
temp = Exponential(1)
sims= convert(Int64, sims)
exp_drw = rand(temp, sims, n)
dir_weights = exp_drw ./sum(exp_drw, 2)
A = zeros(sims, offdiag)
B = zeros((1, offdiag))
C = zeros(sims, p)
mat_point = zeros(p,p)
mat_sig = zeros(p, p)
idx = tril(trues(size(mat_point)), -1);
for i = 1:sims
t = dir_weights[i, 1:n]
d = X .* sqrt.(t)
bb_cov = (d' * d)
inv_cov = inv(bb_cov)
d = diag(inv_cov)
inv_cov[diagind(inv_cov)] = 0
test = inv_cov[idx]
rmv_0 = filter!(x->x≠ 0,test)
A[i,:] = rmv_0
C[i,:] = d
end
for i in 1:loop_max
temp1 = post_prob(A[:,i])
B[1,i] = ifelse(temp1 > prob, 1, 0)
end
mat_sig[idx] = B
mat_point[idx] = mean(A, 1)
mat_point[diagind(mat_point)] = mean(C, 1)
force_sym(mat_sig), force_sym(mat_point), A, B
end
##########################
###### Run model ##########
##########################
using Distributions
p = 500
n = 20
my_mvn = MvNormal(eye(p))
X = rand(my_mvn, n)'
@time test = BB_est(X, 500, .95)
```