In [1]:
versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, haswell)
Environment:
  JULIA_CMDSTAN_HOME = C:\CmdStan
  JULIA_NUM_THREADS = 4
  JULIA_PKGDIR = C:\JuliaPkg
In [2]:
using Distributions
using Random: randn!
using RandomNumbers.Xorshifts
const XOROSHIRO128PLUS = Xoroshiro128Plus()

using Plots
ENV["PLOTS_TEST"] = "true"
gr()
Out[2]:
Plots.GRBackend()
In [3]:
function summary_MF(M, F)
    MM = M[.!(isnan.(M))]
    FF = F[.!(isnan.(F))]
    @show median(MM), median(FF)
    @show mean(MM), mean(FF)
    @show std(MM), std(FF)
    @show count(isnan.(M)), count(isnan.(F))
    @show count(M .> F)
    @show count(M .< F)
end

function plot_MF(M, F;
        bins=range(max(2.0, mean(dist)-4std(FF)), mean(dist)+4std(FF), length=50),
        figsize=(500, 300),
        kwargs...)
    MM = M[.!(isnan.(M))]
    FF = F[.!(isnan.(F))]
    histogram([MM, FF], label=["M", "F"], alpha=0.7; bins=bins, size=figsize, kwargs...)
end
Out[3]:
plot_MF (generic function with 1 method)
In [4]:
function mfmean(; n=10^6, N=1000, p=0.3, t=2.0)
    Mmean = zeros(n)
    Fmean = zeros(n)
    M = zeros(N)
    F = zeros(N)
    for i in 1:n
        randn!(M)
        randn!(F)
        f = F[F .> t]
        Mmean[i] = mean(M[M .> t])
        Fmean[i] = mean(@view f[1:round(Int, length(f)*p)])
    end
    Mmean, Fmean
end

@time M, F = mfmean(; n=10^6, N=1000, p=0.3, t=2.0)
summary_MF(M, F)
plot_MF(M, F, bins=range(2.0, 3.0, length=50))
 29.046963 seconds (13.98 M allocations: 9.127 GiB, 4.87% gc time)
(median(MM), median(FF)) = (2.3692721938100814, 2.3600456745956366)
(mean(MM), mean(FF)) = (2.3731864361142145, 2.3732336567502634)
(std(MM), std(FF)) = (0.07253761787211951, 0.13238713490179752)
(count(isnan.(M)), count(isnan.(F))) = (0, 0)
count(M .> F) = 524741
count(M .< F) = 475259
Out[4]:
2.00 2.25 2.50 2.75 3.00 0 2.5×10 4 5.0×10 4 7.5×10 4 1.0×10 5 M F
In [5]:
function mfmean_Xoroshiro128Plus(; n=10^6, N=1000, p=0.3, t=2.0)
    Mmean = zeros(n)
    Fmean = zeros(n)
    M = zeros(N)
    F = zeros(N)
    for i in 1:n
        randn!(XOROSHIRO128PLUS, M)
        randn!(XOROSHIRO128PLUS, F)
        f = F[F .> t]
        Mmean[i] = mean(M[M .> t])
        Fmean[i] = mean(@view f[1:round(Int, length(f)*p)])
    end
    Mmean, Fmean
end

@time M, F = mfmean_Xoroshiro128Plus(; n=10^6, N=1000, p=0.3, t=2.0)
summary_MF(M, F)
plot_MF(M, F, bins=range(2.0, 3.0, length=50))
 15.061378 seconds (13.26 M allocations: 9.092 GiB, 7.26% gc time)
(median(MM), median(FF)) = (2.369391029041538, 2.3600583244468893)
(mean(MM), mean(FF)) = (2.3732821067912475, 2.3732173869949995)
(std(MM), std(FF)) = (0.07254981589802849, 0.13262921798920727)
(count(isnan.(M)), count(isnan.(F))) = (0, 0)
count(M .> F) = 524674
count(M .< F) = 475326
Out[5]:
2.00 2.25 2.50 2.75 3.00 0 2.5×10 4 5.0×10 4 7.5×10 4 1.0×10 5 M F
In [6]:
function mfmean_Xoroshiro128Plus_Ver2(; n=10^6, N=1000, p=0.3, t=2.0)
    Mmean = zeros(n)
    Fmean = zeros(n)
    M = zeros(N)
    F = zeros(N)
    for i in 1:n
        randn!(XOROSHIRO128PLUS, M)
        randn!(XOROSHIRO128PLUS, F)
        Mmean[i] = mean(M[k] for k in 1:N if M[k]  t)
        m = 0
        for k in 1:N
            if F[k]  t
                m += 1
                F[m] = F[k]
            end
        end
        Fmean[i] = mean(F[k] for k in 1:round(Int, p*m))
    end
    Mmean, Fmean
end

@time M, F = mfmean_Xoroshiro128Plus_Ver2(; n=10^6, N=1000, p=0.3, t=2.0)
summary_MF(M, F)
plot_MF(M, F, bins=range(2.0, 3.0, length=50))
 10.144046 seconds (6.18 M allocations: 176.866 MiB, 0.26% gc time)
(median(MM), median(FF)) = (2.369333090186332, 2.3598376135652868)
(mean(MM), mean(FF)) = (2.373258705036385, 2.3731435435046486)
(std(MM), std(FF)) = (0.07251480489113193, 0.13242613895355249)
(count(isnan.(M)), count(isnan.(F))) = (0, 0)
count(M .> F) = 525047
count(M .< F) = 474953
Out[6]:
2.00 2.25 2.50 2.75 3.00 0 2.5×10 4 5.0×10 4 7.5×10 4 1.0×10 5 M F
In [7]:
function sim_tnd_MF(; t=2.0, N=1000, p=0.3, n=10^4)
    M = zeros(n)
    F = zeros(n)
    tnd = TruncatedNormal(0.0, 1.0, t, Inf)
    p_bin = ccdf(Normal(0.0, 1.0), t)
    bin = Binomial(N, p_bin)
    for i in 1:n
        M[i] = mean(rand(tnd, rand(bin)))
        F[i] = mean(rand(tnd, round(Int, p*rand(bin))))
    end
    tnd, M, F
end
sim_tnd_MF(n=10^3) # compile
GC.gc()

@time tnd, M, F = sim_tnd_MF(t=2.0, N=1000, p=0.3, n=10^6)
summary_MF(M, F)
plot_MF(M, F, bins=range(2.0, 3.0, length=50))
  2.736865 seconds (2.02 M allocations: 403.416 MiB, 1.34% gc time)
(median(MM), median(FF)) = (2.3693389346883755, 2.360191465060133)
(mean(MM), mean(FF)) = (2.373175129868477, 2.373357467179366)
(std(MM), std(FF)) = (0.07249434710518572, 0.13247871435315237)
(count(isnan.(M)), count(isnan.(F))) = (0, 0)
count(M .> F) = 524058
count(M .< F) = 475942
Out[7]:
2.00 2.25 2.50 2.75 3.00 0 2.5×10 4 5.0×10 4 7.5×10 4 1.0×10 5 M F
In [8]:
function sim_tnd_MF_Ver2(; t=2.0, N=1000, p=0.3, n=10^4)
    M = zeros(n)
    F = zeros(n)
    tnd = TruncatedNormal(0.0, 1.0, t, Inf)
    p_bin = ccdf(Normal(0.0, 1.0), t)
    bin = Binomial(N, p_bin)
    for i in 1:n
        M[i] = mean(k->rand(tnd), 1:rand(bin))
        F[i] = mean(k->rand(tnd), 1:round(Int, p*rand(bin)))
    end
    tnd, M, F
end
sim_tnd_MF_Ver2(n=10^3) # compile
GC.gc()

@time tnd, M, F = sim_tnd_MF_Ver2(t=2.0, N=1000, p=0.3, n=10^6)
summary_MF(M, F)
plot_MF(M, F, bins=range(2.0, 3.0, length=50))
  2.045361 seconds (9.46 k allocations: 15.780 MiB)
(median(MM), median(FF)) = (2.369385472835822, 2.3600769716008676)
(mean(MM), mean(FF)) = (2.3732029684624014, 2.373334869690733)
(std(MM), std(FF)) = (0.07249116004204485, 0.13263651358155584)
(count(isnan.(M)), count(isnan.(F))) = (0, 0)
count(M .> F) = 524815
count(M .< F) = 475185
Out[8]:
2.00 2.25 2.50 2.75 3.00 0 2.5×10 4 5.0×10 4 7.5×10 4 1.0×10 5 M F
In [9]:
function simpi(n)
    s = 0
    for i in 1:n
        s += ifelse(rand()^2 + rand()^2 < 1.0, 1, 0)
    end
    4s/n
end
simpi(10^5) # compile

@time simpi(10^7)
  0.038209 seconds (1.93 k allocations: 115.887 KiB)
Out[9]:
3.1421304
In [10]:
@time simpi(10^8)
  0.356153 seconds (1.93 k allocations: 115.887 KiB)
Out[10]:
3.14132636
In [11]:
@time simpi(10^9)
  3.792615 seconds (6 allocations: 192 bytes)
Out[11]:
3.141545184
In [12]:
function simpi_Xoroshiro128Plus(n)
    s = 0
    for i in 1:n
        s += ifelse(rand(XOROSHIRO128PLUS)^2 + rand(XOROSHIRO128PLUS)^2 < 1.0, 1, 0)
    end
    4s/n
end
simpi_Xoroshiro128Plus(10^5) # compile

@time simpi_Xoroshiro128Plus(10^7)
  0.061850 seconds (6 allocations: 192 bytes)
Out[12]:
3.140606
In [13]:
@time simpi_Xoroshiro128Plus(10^8)
  0.622216 seconds (6 allocations: 192 bytes)
Out[13]:
3.14168644
In [14]:
@time simpi_Xoroshiro128Plus(10^9)
  6.105230 seconds (6 allocations: 192 bytes)
Out[14]:
3.141479288
In [15]:
simpi_count(n) = 4count(rand(n).^2 .+ rand(n).^2 .< 1.0)/n
simpi_count(10^5) # compile

@time simpi_count(10^7)
  0.191046 seconds (16 allocations: 153.785 MiB, 55.34% gc time)
Out[15]:
3.1410048
In [16]:
@time simpi_count(10^8)
  1.045771 seconds (16 allocations: 1.502 GiB, 18.86% gc time)
Out[16]:
3.14158204
In [17]:
simpi_Xoroshiro128Plus_count(n) = 4count(rand(XOROSHIRO128PLUS,n).^2 .+ rand(XOROSHIRO128PLUS,n).^2 .< 1.0)/n
simpi_Xoroshiro128Plus_count(10^5) # compile

@time simpi_Xoroshiro128Plus_count(10^7)
  0.327133 seconds (16 allocations: 153.785 MiB, 70.27% gc time)
Out[17]:
3.1412792
In [18]:
@time simpi_Xoroshiro128Plus_count(10^8)
  1.168187 seconds (16 allocations: 1.502 GiB, 19.40% gc time)
Out[18]:
3.141661
In [ ]: