You can go a bit faster with integer arithmetic:

```
function mean_loop_int(p::AbstractArray, q::AbstractArray)
z = 0
for x in q
z += count(>(x), p)
end
return z / length(p)
end
```

```
julia> x = rand(10^4);
julia> y = rand(10^4);
julia> @btime mean_loop($x, $y)
19.034 ms (0 allocations: 0 bytes)
4968.276699999993
julia> @btime mean_loop_int($x, $y)
15.797 ms (0 allocations: 0 bytes)
4968.2767
```

But if you change the algorithm, you can get significant speedups:

```
function mean_sorted(p, q)
ps = sort(p)
Lq = length(q)
z = 0
for x in q
ind = searchsortedfirst(ps, x)
z += Lq - ind + 1
end
return z / length(p)
end
```

```
julia> @btime mean_sorted($x, $y)
1.512 ms (2 allocations: 78.20 KiB)
4968.2767
```

Possibly, it can be even faster if we sort both arrays and search for the next larger at each iteration, but I couldnâ€™t immediately find a `searchsortednext`

function.

This also scales well for longer vectors. For 10^6 length I get:

```
julia> @btime mean_sorted(x, y) setup=(x=rand(10^6); y=rand(10^6))
395.748 ms (2 allocations: 7.63 MiB)
500344.792228
```

That is, N \log N complexity, rather than N^2.