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
using Distributions
using Random: randn!
using RandomNumbers.Xorshifts
const XOROSHIRO128PLUS = Xoroshiro128Plus()
using Plots
ENV["PLOTS_TEST"] = "true"
gr()
Plots.GRBackend()
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
plot_MF (generic function with 1 method)
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
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
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
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
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
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)
3.1421304
@time simpi(10^8)
0.356153 seconds (1.93 k allocations: 115.887 KiB)
3.14132636
@time simpi(10^9)
3.792615 seconds (6 allocations: 192 bytes)
3.141545184
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)
3.140606
@time simpi_Xoroshiro128Plus(10^8)
0.622216 seconds (6 allocations: 192 bytes)
3.14168644
@time simpi_Xoroshiro128Plus(10^9)
6.105230 seconds (6 allocations: 192 bytes)
3.141479288
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)
3.1410048
@time simpi_count(10^8)
1.045771 seconds (16 allocations: 1.502 GiB, 18.86% gc time)
3.14158204
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)
3.1412792
@time simpi_Xoroshiro128Plus_count(10^8)
1.168187 seconds (16 allocations: 1.502 GiB, 19.40% gc time)
3.141661