From e5c1056688c7ed4473985589de9e79cd0266ddb4 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 3 Jan 2024 11:49:01 +0100 Subject: [PATCH 1/8] ignore settings file --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 29126e4..9b79786 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,5 @@ docs/site/ # committed for packages, but should be committed for applications that require a static # environment. Manifest.toml + +.vscode/settings.json From edf4853ae8c753840fc13cb689c416f75c2311e9 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 3 Jan 2024 16:34:35 +0100 Subject: [PATCH 2/8] add fit_beta_mixture functionality --- Project.toml | 3 ++ src/SafetySignalDetection.jl | 9 +++- src/fit_beta_mixture.jl | 26 +++++++++++ src/fit_mle.jl | 87 +++++++++++++++++++++++++++++++++++ test/runtests.jl | 41 ++--------------- test/test_fit_beta_mixture.jl | 49 ++++++++++++++++++++ test/test_fit_mle.jl | 32 +++++++++++++ test/test_helpers.jl | 15 ++++++ test/test_meta_analytic.jl | 21 +++++++++ 9 files changed, 245 insertions(+), 38 deletions(-) create mode 100644 src/fit_beta_mixture.jl create mode 100644 src/fit_mle.jl create mode 100644 test/test_fit_beta_mixture.jl create mode 100644 test/test_fit_mle.jl create mode 100644 test/test_helpers.jl create mode 100644 test/test_meta_analytic.jl diff --git a/Project.toml b/Project.toml index 5c4b503..afd41c8 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,10 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ExpectationMaximization = "e1fe09cc-5134-44c2-a941-50f4cd97986a" FreqTables = "da1fdf0e-e0ff-5433-a45f-9bb5ff651cb1" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0" diff --git a/src/SafetySignalDetection.jl b/src/SafetySignalDetection.jl index 9fabb6a..3c4cba1 100644 --- a/src/SafetySignalDetection.jl +++ b/src/SafetySignalDetection.jl @@ -3,10 +3,17 @@ module SafetySignalDetection using Turing using StatsPlots using Distributions +using SpecialFunctions +using Statistics +using LinearAlgebra +using ExpectationMaximization export - meta_analytic + meta_analytic, + fit_beta_mixture include("meta_analytic.jl") +include("fit_mle.jl") +include("fit_beta_mixture.jl") end \ No newline at end of file diff --git a/src/fit_beta_mixture.jl b/src/fit_beta_mixture.jl new file mode 100644 index 0000000..7ee76c6 --- /dev/null +++ b/src/fit_beta_mixture.jl @@ -0,0 +1,26 @@ +function init_beta_dists(n_components::Int) + zero_one_range = range(start = 0, stop = 1, length = n_components + 2) + alpha_range = zero_one_range[2:(n_components + 1)] + [Beta(alpha, 1) for alpha in alpha_range] +end + +""" + Fit a beta mixture to a vector of prior samples + +This function returns a beta mixture of `n_components` components approximating +the distribution of the sample vector `x`. +""" +function fit_beta_mixture(x::AbstractArray{T}, n_components::Int) where T<:Real + 0 < n_components || throw(DomainError(n_components, "there must be at least one component")) + # Remove outliers to stabilize the fitting. + lower_quant = quantile!(x, 0.001) + upper_quant = quantile!(x, 0.999) + x = filter(y -> y > lower_quant && y < upper_quant, x) + + # We initialize here with Beta distributions that are not identical but have increasing alpha parameters. + beta_dists = init_beta_dists(n_components) + mix_guess = MixtureModel(beta_dists) + + # Fit the MLE with the classic EM algorithm. + fit_mle(mix_guess, x) +end diff --git a/src/fit_mle.jl b/src/fit_mle.jl new file mode 100644 index 0000000..ab1a671 --- /dev/null +++ b/src/fit_mle.jl @@ -0,0 +1,87 @@ +import Distributions: fit_mle, suffstats, varm + +# Weighted MLE for beta distribution +# This supplements Distributions.jl and is needed for the classic EM algorithm to work for the beta distribution. + +# sufficient statistics - these capture everything of the data we need +struct BetaStats <: SufficientStats + sum_log_x::Float64 # (weighted) sum of log(x) + sum_log_1mx::Float64 # (weighted) sum of log(1 - x) + tw::Float64 # total sample weight + x_bar::Float64 # (weighted) mean of x + v_bar::Float64 # (weighted) variance of x +end + +function suffstats(::Type{<:Beta}, x::AbstractArray{T}, w::AbstractArray{T}) where T<:Real + + tw = 0.0 + sum_log_x = 0.0 * zero(T) + sum_log_1mx = 0.0 * zero(T) + x_bar = 0.0 * zero(T) + + for i in eachindex(x, w) + @inbounds xi = x[i] + 0 < xi < 1 || throw(DomainError(xi, "samples must be larger than 0 and smaller than 1")) + @inbounds wi = w[i] + wi >= 0 || throw(DomainError(wi, "weights must be non-negative")) + isfinite(wi) || throw(DomainError(wi, "weights must be finite")) + @inbounds sum_log_x += wi * log(xi) + @inbounds sum_log_1mx += wi * log(one(T) - xi) + @inbounds x_bar += wi * xi + tw += wi + end + sum_log_x /= tw + sum_log_1mx /= tw + x_bar /= tw + v_bar = varm(x, x_bar) + + BetaStats(sum_log_x, sum_log_1mx, tw, x_bar, v_bar) +end + +# without weights we just put weight 1 for each observation +function suffstats(::Type{<:Beta}, x::AbstractArray{T}) where T<:Real + + w = ones(Float64, length(x)) + suffstats(Beta, x, w) + +end + +# generic fit function based on the sufficient statistics, on the log scale to be robust +function fit_mle(::Type{<:Beta}, ss::BetaStats; + maxiter::Int=1000, tol::Float64=1e-14) + + # Initialization of parameters at the moment estimators (I guess) + temp = ((ss.x_bar * (1 - ss.x_bar)) / ss.v_bar) - 1 + α₀ = ss.x_bar * temp + β₀ = (1 - ss.x_bar) * temp + + g₁ = ss.sum_log_x + g₂ = ss.sum_log_1mx + + θ= [log(α₀) ; log(β₀)] + + converged = false + t=0 + while !converged && t < maxiter + t += 1 + α = exp(θ[1]) + β = exp(θ[2]) + temp1 = digamma(α + β) + temp2 = trigamma(α + β) + temp3 = g₁ + temp1 - digamma(α) + grad = [temp3 * α + (temp1 + g₂ - digamma(β)) * β] + hess = [((temp2 - trigamma(α)) * α + temp3) * α temp2 * β * α + temp2 * α * β ((temp2 - trigamma(β)) * β + temp1 + g₂ - digamma(β)) * β ] + Δθ = hess\grad #newton step + θ .-= Δθ + converged = dot(Δθ,Δθ) < 2*tol #stopping criterion + end + + α = exp(θ[1]) + β = exp(θ[2]) + return Beta(α, β) +end + +# then define methods for the original data +fit_mle(::Type{<:Beta}, x::AbstractArray{T}, w::AbstractArray{T}; maxiter::Int=1000, tol::Float64=1e-14) where T<:Real = fit_mle(Beta, suffstats(Beta, x, w)) diff --git a/test/runtests.jl b/test/runtests.jl index f1c396d..d8bab6f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,40 +6,7 @@ using DataFrames using Turing using SafetySignalDetection -# Helper function for numerical tests. -# Taken from https://github.com/TuringLang/Turing.jl/blob/master/test/test_utils/numerical_tests.jl#L41 for now. -function check_numerical(chain, - symbols::Vector, - exact_vals::Vector; - atol=0.2, - rtol=0.0) - for (sym, val) in zip(symbols, exact_vals) - E = val isa Real ? - mean(chain[sym]) : - vec(mean(chain[sym], dims=1)) - @info (symbol=sym, exact=val, evaluated=E) - @test E ≈ val atol=atol rtol=rtol - end -end - -@testset "meta_analytic.jl" begin - rng = StableRNG(123) - - n_trials = 5 - n_patients = 50 - df = DataFrame( - y = rand(rng, Bernoulli(0.2), n_trials * n_patients), - time = rand(rng, Exponential(1), n_trials * n_patients), - trialindex = repeat(1:n_trials, n_patients) - ) - - chain = sample( - rng, - meta_analytic(df.y, df.time, df.trialindex, Beta(2, 8), Beta(9, 10)), - HMC(0.05, 10), - 1000 - ) - - check_numerical(chain, [:a], [0.223], rtol=0.001) - check_numerical(chain, [:b], [0.485], rtol=0.001) -end +include("test_helpers.jl") +include("test_meta_analytic.jl") +include("test_fit_beta_mixture.jl") +include("test_fit_mle.jl") diff --git a/test/test_fit_beta_mixture.jl b/test/test_fit_beta_mixture.jl new file mode 100644 index 0000000..137add2 --- /dev/null +++ b/test/test_fit_beta_mixture.jl @@ -0,0 +1,49 @@ +@testset "init_beta_dists(1)" begin + result = SafetySignalDetection.init_beta_dists(1) + @test typeof(result) <: Array + @test length(result) == 1 + @test typeof(result[1]) <: Beta + @test result[1].α ≈ 0.5 + @test result[1].β == 1 +end + +@testset "init_beta_dists(4)" begin + result = SafetySignalDetection.init_beta_dists(4) + @test typeof(result) <: Array + @test length(result) == 4 + @test result[1].α ≈ 0.2 + @test result[1].β == 1 + @test result[2].α ≈ 0.4 + @test result[2].β == 1 + @test result[3].α ≈ 0.6 + @test result[3].β == 1 + @test result[4].α ≈ 0.8 + @test result[4].β == 1 +end + +@testset "fit_beta_mixture" begin + + Random.seed!(123) + + N = 100_000 + α₁ = 10 + β₁ = 5 + α₂ = 5 + β₂ = 10 + π = 0.6 + + # Draw large sample from true beta mixture model. + mix_true = MixtureModel([Beta(α₁, β₁), Beta(α₂, β₂)], [π, 1 - π]) + x = rand(mix_true, N) + + # Fit the beta mixture. + mix_fit = fit_beta_mixture(x, 2) + + # Compare parameter estimates with true parameter values. + @test mix_fit.components[1].α ≈ α₂ rtol = 0.1 + @test mix_fit.components[1].β ≈ β₂ rtol = 0.1 + @test mix_fit.components[2].α ≈ α₁ rtol = 0.1 + @test mix_fit.components[2].β ≈ β₁ rtol = 0.1 + @test mix_fit.prior.p[1] ≈ 1 - π rtol = 0.01 + @test mix_fit.prior.p[2] ≈ π rtol = 0.01 +end diff --git a/test/test_fit_mle.jl b/test/test_fit_mle.jl new file mode 100644 index 0000000..e6484b7 --- /dev/null +++ b/test/test_fit_mle.jl @@ -0,0 +1,32 @@ +@testset "fit_mle.jl" begin + + Random.seed!(1234) + + # Try to fit a beta mixture. + N = 50_000 + α₁ = 10 + β₁ = 5 + α₂ = 5 + β₂ = 10 + π = 0.3 + + # Mixture Model of two betas. + mix_true = MixtureModel([Beta(α₁, β₁), Beta(α₂, β₂)], [π, 1 - π]) + + # Generate N samples from the mixture. + y = rand(mix_true, N) + + # Initial guess. + mix_guess = MixtureModel([Beta(0.33, 1), Beta(0.66, 1)], [0.5, 1 - 0.5]) + + # Fit the MLE with the classic EM algorithm: + mix_mle = fit_mle(mix_guess, y) + + # Compare parameter estimates with true parameter values. + @test mix_mle.components[1].α ≈ α₂ rtol = 0.1 + @test mix_mle.components[1].β ≈ β₂ rtol = 0.1 + @test mix_mle.components[2].α ≈ α₁ rtol = 0.1 + @test mix_mle.components[2].β ≈ β₁ rtol = 0.1 + @test mix_mle.prior.p[1] ≈ 1 - π rtol = 0.01 + @test mix_mle.prior.p[2] ≈ π rtol = 0.01 +end \ No newline at end of file diff --git a/test/test_helpers.jl b/test/test_helpers.jl new file mode 100644 index 0000000..fa0de91 --- /dev/null +++ b/test/test_helpers.jl @@ -0,0 +1,15 @@ +# Helper function for numerical tests. +# Taken from https://github.com/TuringLang/Turing.jl/blob/master/test/test_utils/numerical_tests.jl#L41 for now. +function check_numerical(chain, + symbols::Vector, + exact_vals::Vector; + atol=0.2, + rtol=0.0) + for (sym, val) in zip(symbols, exact_vals) + E = val isa Real ? + mean(chain[sym]) : + vec(mean(chain[sym], dims=1)) + @info (symbol=sym, exact=val, evaluated=E) + @test E ≈ val atol=atol rtol=rtol + end +end diff --git a/test/test_meta_analytic.jl b/test/test_meta_analytic.jl new file mode 100644 index 0000000..4674f86 --- /dev/null +++ b/test/test_meta_analytic.jl @@ -0,0 +1,21 @@ +@testset "meta_analytic.jl" begin + rng = StableRNG(123) + + n_trials = 5 + n_patients = 50 + df = DataFrame( + y = rand(rng, Bernoulli(0.2), n_trials * n_patients), + time = rand(rng, Exponential(1), n_trials * n_patients), + trialindex = repeat(1:n_trials, n_patients) + ) + + chain = sample( + rng, + meta_analytic(df.y, df.time, df.trialindex, Beta(2, 8), Beta(9, 10)), + HMC(0.05, 10), + 1000 + ) + + check_numerical(chain, [:a], [0.223], rtol=0.001) + check_numerical(chain, [:b], [0.485], rtol=0.001) +end From 67f01c3747fec9d9024a192eaaf7d2a6e7da0633 Mon Sep 17 00:00:00 2001 From: Kristian Brock Date: Wed, 10 Jan 2024 11:25:01 +0000 Subject: [PATCH 3/8] Added tests of fit_beta_mixture and meta_analytic to reproduce known examples from R --- .gitignore | 2 + Project.toml | 3 +- design/ReconcileWithR.jl | 6 +- test/Project.toml | 1 + test/large_historic_dataset.csv | 1001 +++++++++++++++++++++++++++++++ test/runtests.jl | 1 + test/small_historic_dataset.csv | 101 ++++ test/test_fit_beta_mixture.jl | 54 ++ test/test_meta_analytic.jl | 39 ++ 9 files changed, 1204 insertions(+), 4 deletions(-) create mode 100644 test/large_historic_dataset.csv create mode 100644 test/small_historic_dataset.csv diff --git a/.gitignore b/.gitignore index 9b79786..c4f7916 100644 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,5 @@ docs/site/ Manifest.toml .vscode/settings.json + +.DS_Store \ No newline at end of file diff --git a/Project.toml b/Project.toml index afd41c8..a9ca883 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" -StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0" diff --git a/design/ReconcileWithR.jl b/design/ReconcileWithR.jl index 2376f7b..e599846 100644 --- a/design/ReconcileWithR.jl +++ b/design/ReconcileWithR.jl @@ -335,7 +335,7 @@ Turing = "~0.29.3" PLUTO_MANIFEST_TOML_CONTENTS = """ # This file is machine-generated - editing it directly is not advised -julia_version = "1.9.3" +julia_version = "1.9.1" manifest_format = "2.0" project_hash = "4359e7ffd037ffcea031434e2d83a0afe89b7541" @@ -684,7 +684,7 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.5+0" +version = "1.0.2+0" [[deps.CompositionsBase]] git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad" @@ -1623,7 +1623,7 @@ version = "0.42.2+0" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.9.2" +version = "1.9.0" [[deps.PlotThemes]] deps = ["PlotUtils", "Statistics"] diff --git a/test/Project.toml b/test/Project.toml index 8e949d1..e1663e4 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [deps] +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/test/large_historic_dataset.csv b/test/large_historic_dataset.csv new file mode 100644 index 0000000..9f65459 --- /dev/null +++ b/test/large_historic_dataset.csv @@ -0,0 +1,1001 @@ +"","y","trial","time" +"1",0,1,0.508060037142824 +"2",0,7,0.5208363096442 +"3",0,2,0.556331578741567 +"4",0,10,0.563791380921675 +"5",0,8,0.571754800768225 +"6",0,7,0.590052150745623 +"7",0,5,0.591040918682293 +"8",0,10,0.595221235355885 +"9",0,7,0.597784269735739 +"10",0,1,0.607071721277385 +"11",0,5,0.611330809521675 +"12",1,4,0.613947846980951 +"13",0,2,0.618876835716143 +"14",0,9,0.640354126490284 +"15",0,1,0.645504716144507 +"16",0,7,0.648776558221881 +"17",0,2,0.653263780100834 +"18",0,9,0.664869826653364 +"19",0,10,0.671171271026457 +"20",0,2,0.674994900783847 +"21",0,2,0.679323147659513 +"22",0,2,0.681681905904048 +"23",0,7,0.693983536523541 +"24",1,3,0.713401719711471 +"25",0,8,0.715330274447549 +"26",0,7,0.715554267597332 +"27",0,4,0.724320719842258 +"28",0,2,0.743751407473829 +"29",0,9,0.744621157756557 +"30",0,7,0.749986377628643 +"31",1,7,0.756438894458238 +"32",0,8,0.769038580022455 +"33",0,10,0.781380031033633 +"34",0,7,0.782826549365607 +"35",0,3,0.78411009277818 +"36",1,5,0.791492904824573 +"37",1,6,0.801637802299991 +"38",0,8,0.811791839104605 +"39",1,3,0.837362781844232 +"40",1,6,0.841242652393076 +"41",0,1,0.862928944447668 +"42",0,1,0.869636503880502 +"43",0,5,0.886621511367525 +"44",0,2,0.889937450386028 +"45",0,3,0.89081983407345 +"46",0,1,0.894065672679249 +"47",0,3,0.895545855571809 +"48",0,3,0.909240301239177 +"49",0,6,0.919524707524235 +"50",0,7,0.924650129098574 +"51",0,5,0.945304040451132 +"52",0,4,0.946978089384413 +"53",0,3,0.951859002214535 +"54",0,2,0.957347076003118 +"55",0,2,0.960843806894175 +"56",1,6,0.965779375707409 +"57",0,3,0.966733577783565 +"58",0,4,0.96968774609861 +"59",0,5,0.972657807600261 +"60",0,9,0.979699684335007 +"61",0,9,0.979834431354125 +"62",0,7,0.990429732491942 +"63",0,4,0.993043940746905 +"64",0,5,0.9954099067624 +"65",0,7,0.997159593415083 +"66",0,1,1.00392896736284 +"67",0,7,1.01355110456085 +"68",0,8,1.04904199831044 +"69",0,9,1.06484850264858 +"70",0,8,1.06794773782012 +"71",0,2,1.0712293948771 +"72",0,3,1.08236479130322 +"73",0,3,1.0874704878659 +"74",0,5,1.08815656886374 +"75",1,8,1.09446561418411 +"76",0,7,1.11639969562923 +"77",0,4,1.13149712565995 +"78",1,5,1.15531456740156 +"79",1,10,1.17462593214085 +"80",0,8,1.18276276532017 +"81",0,10,1.18776627686728 +"82",0,2,1.20208448878763 +"83",0,7,1.21019494540132 +"84",0,10,1.21124169832635 +"85",0,10,1.2113759552136 +"86",0,10,1.21897941682307 +"87",0,1,1.22070786681776 +"88",0,10,1.24966304313974 +"89",0,3,1.2527988999829 +"90",0,5,1.25518082340978 +"91",0,4,1.25891242750253 +"92",0,5,1.26357294799558 +"93",0,3,1.26469715973161 +"94",0,2,1.27593270482697 +"95",0,4,1.28273087209758 +"96",0,3,1.29175769180613 +"97",0,4,1.31240415124692 +"98",0,6,1.31259288685395 +"99",0,9,1.33657055292152 +"100",0,1,1.34760702875274 +"101",0,1,1.35092539497948 +"102",0,4,1.35149998950137 +"103",1,4,1.36380995980719 +"104",0,1,1.3699089969072 +"105",0,5,1.39162108824477 +"106",0,9,1.41320690895193 +"107",0,4,1.42429481301177 +"108",0,3,1.43756269467673 +"109",0,6,1.45867396387275 +"110",0,9,1.4843040235794 +"111",0,7,1.50413489949682 +"112",0,6,1.50476468195466 +"113",0,2,1.50792893654395 +"114",1,2,1.51427706286971 +"115",0,2,1.51990337487504 +"116",0,1,1.52660963020016 +"117",1,9,1.54658357808942 +"118",0,2,1.55312719334958 +"119",0,9,1.56168943916997 +"120",0,7,1.56724044201979 +"121",0,2,1.56960279979827 +"122",0,4,1.57061575538562 +"123",0,8,1.58913724848069 +"124",0,9,1.63340894904666 +"125",1,7,1.63453285789541 +"126",0,6,1.63815252316143 +"127",0,8,1.63820153185591 +"128",0,8,1.63833771414927 +"129",0,7,1.64856892239365 +"130",0,1,1.64864410723622 +"131",0,2,1.65843753930201 +"132",0,10,1.67489939937304 +"133",0,6,1.67929099241295 +"134",0,6,1.6844599054337 +"135",0,10,1.68751908536558 +"136",0,2,1.6878555520203 +"137",0,2,1.69072741575521 +"138",0,1,1.70688129819827 +"139",0,6,1.70770595916445 +"140",1,7,1.73290502698305 +"141",0,6,1.74700937115527 +"142",0,6,1.76256142666764 +"143",0,3,1.77626084866161 +"144",0,4,1.78275493509761 +"145",0,9,1.79451785950573 +"146",0,1,1.79683551605908 +"147",0,4,1.79896815308863 +"148",0,10,1.82751915366714 +"149",0,3,1.82966782722829 +"150",0,1,1.83625485336615 +"151",0,1,1.83886171092779 +"152",0,6,1.83909025857822 +"153",0,8,1.8945210277692 +"154",0,10,1.89926617033151 +"155",0,9,1.90673089127317 +"156",1,1,1.91815551343426 +"157",1,3,1.92465540812868 +"158",0,3,1.93794650956532 +"159",0,1,1.95599599109422 +"160",0,10,1.96806005083641 +"161",0,1,1.99085121198116 +"162",0,10,2.02172043292174 +"163",0,8,2.04005428627372 +"164",0,8,2.0400546436315 +"165",0,1,2.04660221272852 +"166",0,7,2.05084673782118 +"167",1,7,2.05530288144633 +"168",0,10,2.05650863301627 +"169",1,8,2.05721580272336 +"170",0,9,2.06629428834408 +"171",0,10,2.07433926135219 +"172",0,1,2.07517725225008 +"173",0,7,2.07702938700236 +"174",1,9,2.08349003238744 +"175",0,5,2.08545468252293 +"176",0,6,2.08843742971915 +"177",0,5,2.09125745762232 +"178",1,2,2.09573347421896 +"179",0,3,2.09584005296818 +"180",1,2,2.11687577984 +"181",1,6,2.12263871705795 +"182",0,4,2.13031087978406 +"183",0,9,2.14804437796164 +"184",0,3,2.1556787161676 +"185",0,3,2.165854880378 +"186",0,8,2.17743488218256 +"187",0,6,2.17963428347458 +"188",1,4,2.1827845232248 +"189",1,2,2.18705152533963 +"190",1,8,2.20666564138536 +"191",1,2,2.22075810728068 +"192",1,4,2.23536079702288 +"193",0,8,2.23732688275952 +"194",1,8,2.24674480541422 +"195",0,9,2.27111616242821 +"196",0,3,2.29707248764131 +"197",1,4,2.30269893956389 +"198",0,6,2.31149708800925 +"199",0,5,2.31759316307434 +"200",0,4,2.31800139362399 +"201",0,5,2.32487089927715 +"202",0,4,2.32491497969662 +"203",0,6,2.3379459227658 +"204",0,9,2.34455429837007 +"205",0,9,2.34955306482866 +"206",1,6,2.35146895254396 +"207",0,7,2.3746241570288 +"208",0,10,2.39080330209961 +"209",0,9,2.39597046922116 +"210",0,8,2.40253279864943 +"211",1,3,2.45109375973492 +"212",0,9,2.45658902595546 +"213",0,9,2.46734397009312 +"214",1,9,2.48540336882281 +"215",0,6,2.4910158662084 +"216",0,5,2.50290516720628 +"217",0,4,2.51808304840159 +"218",0,6,2.53566151463734 +"219",0,8,2.54047373802708 +"220",0,8,2.55586121898709 +"221",0,10,2.5642575815931 +"222",1,1,2.57043802591018 +"223",1,1,2.57599809542275 +"224",0,10,2.5798014720853 +"225",0,4,2.60188409788506 +"226",0,4,2.60494448204963 +"227",0,9,2.61898731847099 +"228",0,9,2.63257886103118 +"229",1,1,2.63513932952421 +"230",0,8,2.64346689883851 +"231",0,9,2.65533703684569 +"232",0,6,2.65946364219529 +"233",0,10,2.68330744792809 +"234",0,10,2.6837927643965 +"235",1,6,2.69504299648166 +"236",0,4,2.69750349377808 +"237",1,5,2.71813140199961 +"238",0,10,2.72540949646783 +"239",0,6,2.7337994036673 +"240",1,1,2.76804693913277 +"241",1,3,2.77140597380422 +"242",0,9,2.77290830136001 +"243",0,5,2.7735271743052 +"244",0,10,2.78258184561654 +"245",1,2,2.80129256006394 +"246",0,9,2.81827811257893 +"247",0,9,2.81983426739367 +"248",1,10,2.82804965538767 +"249",0,9,2.82810524432317 +"250",0,3,2.84639902886248 +"251",0,1,2.84925574623759 +"252",0,9,2.87971451375761 +"253",1,10,2.87987368242929 +"254",0,9,2.88013422791494 +"255",0,8,2.89971238157276 +"256",0,1,2.92061296746504 +"257",1,8,2.92074277455822 +"258",0,6,2.9233737225914 +"259",0,7,2.94507460116186 +"260",0,2,2.9473715831224 +"261",0,5,2.95134313029946 +"262",0,8,2.95272876072452 +"263",0,4,2.97551104427898 +"264",0,7,2.99366881401777 +"265",0,1,3.01431784646069 +"266",0,6,3.01519259335422 +"267",0,9,3.02951112653925 +"268",0,9,3.03581230573687 +"269",1,6,3.05083710885716 +"270",0,6,3.05854022654005 +"271",1,4,3.05966864689951 +"272",0,4,3.06121246394746 +"273",1,4,3.06423588166989 +"274",1,6,3.06465768529147 +"275",0,6,3.13595963650939 +"276",0,7,3.14544568191796 +"277",0,1,3.16693777066273 +"278",0,2,3.17207865273663 +"279",0,6,3.17525053680026 +"280",0,9,3.18935173354711 +"281",1,10,3.1938271652874 +"282",0,5,3.21713316557471 +"283",1,10,3.22868540310733 +"284",1,6,3.25131989888221 +"285",0,8,3.25983650723736 +"286",0,8,3.28717629056562 +"287",0,7,3.29116387667877 +"288",1,7,3.31084625825945 +"289",1,4,3.31215629099596 +"290",1,7,3.33162725328111 +"291",0,3,3.34626008483491 +"292",0,5,3.35002290078349 +"293",0,6,3.35019346835467 +"294",0,4,3.37220488421868 +"295",0,6,3.3771032069369 +"296",1,1,3.40380940987852 +"297",1,6,3.45694791973398 +"298",0,2,3.47837917918132 +"299",0,7,3.48030374443093 +"300",1,5,3.48075126875126 +"301",0,3,3.48585896527794 +"302",0,9,3.49310717629626 +"303",1,9,3.50248671890673 +"304",1,9,3.50566223939412 +"305",0,1,3.52810345527389 +"306",0,8,3.5328647274957 +"307",0,6,3.53512367981159 +"308",0,4,3.54032524262228 +"309",1,2,3.54436489179869 +"310",1,9,3.54780218228837 +"311",0,10,3.59278868023193 +"312",0,5,3.5960429071484 +"313",1,1,3.60462725324165 +"314",1,8,3.60586257139295 +"315",0,5,3.62119302156328 +"316",1,9,3.63762410260716 +"317",0,3,3.6402407037744 +"318",0,10,3.64716014479731 +"319",0,7,3.65361189287135 +"320",1,9,3.65734039349926 +"321",1,7,3.65747485266977 +"322",0,2,3.66299238756744 +"323",0,5,3.66509032289476 +"324",0,3,3.67355190086115 +"325",0,1,3.67683449374449 +"326",0,5,3.69425526747333 +"327",0,4,3.69629858776314 +"328",0,2,3.69663932437861 +"329",1,8,3.71885898836615 +"330",0,3,3.71970924137353 +"331",0,10,3.75454565528763 +"332",0,1,3.76037113388517 +"333",1,7,3.76305392502324 +"334",0,10,3.76464380312228 +"335",0,5,3.7733639724407 +"336",1,10,3.77884087786485 +"337",0,6,3.78616970414435 +"338",1,4,3.79712185757332 +"339",0,3,3.79832820022576 +"340",1,8,3.80290510730449 +"341",0,6,3.80309731377684 +"342",0,9,3.83051218454064 +"343",0,3,3.83198435304956 +"344",0,6,3.84907601731083 +"345",1,8,3.86925456274402 +"346",0,10,3.88142107523106 +"347",1,9,3.89739223632005 +"348",1,7,3.9021709533377 +"349",0,5,3.91106849498692 +"350",0,10,3.91596833802871 +"351",0,10,3.92942391141387 +"352",1,5,3.93193865735904 +"353",0,10,3.93950661153379 +"354",0,3,3.95891881185346 +"355",0,8,3.9637291424139 +"356",0,8,3.99006790130967 +"357",1,7,4.00296279334058 +"358",0,7,4.0051040585641 +"359",0,5,4.01808977629241 +"360",1,10,4.01972032540948 +"361",1,10,4.02040480726869 +"362",1,6,4.02931781186219 +"363",0,8,4.04039035803837 +"364",0,2,4.04988518995531 +"365",0,1,4.06648601581557 +"366",0,2,4.06988519546324 +"367",0,5,4.08460971920416 +"368",1,5,4.09561695125508 +"369",0,3,4.09899210281293 +"370",0,5,4.11111030852934 +"371",0,9,4.11243553992873 +"372",0,5,4.11960386644587 +"373",0,5,4.13000077403578 +"374",1,4,4.13350078996852 +"375",1,10,4.14909292236468 +"376",1,4,4.16182628188738 +"377",0,10,4.16227189032703 +"378",1,6,4.16514378847896 +"379",1,5,4.16701930391479 +"380",0,2,4.18173941674039 +"381",1,4,4.19512187843497 +"382",0,2,4.20872042230574 +"383",1,3,4.23196637332461 +"384",0,9,4.23729535646757 +"385",1,1,4.24211593027673 +"386",0,6,4.24345785780325 +"387",1,8,4.26320547810187 +"388",0,7,4.27504713804278 +"389",0,4,4.28760679555832 +"390",0,6,4.29156240667922 +"391",0,7,4.29560829018017 +"392",1,9,4.31802014073937 +"393",0,2,4.31895055373761 +"394",1,5,4.33547468995517 +"395",1,2,4.34726702270829 +"396",0,4,4.34739736173568 +"397",0,10,4.35418548980993 +"398",1,8,4.35664665515559 +"399",0,4,4.36931000473989 +"400",0,2,4.3751571972257 +"401",0,8,4.40813038605642 +"402",1,9,4.42734543073284 +"403",1,5,4.42841078658416 +"404",0,6,4.43029730579178 +"405",1,1,4.44051998568018 +"406",1,9,4.45637944088484 +"407",0,7,4.45666368851904 +"408",0,1,4.46714097828096 +"409",0,8,4.46726219624704 +"410",1,6,4.47468479264947 +"411",0,3,4.49722221294711 +"412",0,2,4.50677848891267 +"413",1,6,4.51043304013 +"414",0,9,4.51439605141902 +"415",1,7,4.54728639679705 +"416",1,8,4.54816542314516 +"417",1,2,4.55663539113586 +"418",1,9,4.56794953162849 +"419",0,3,4.57729355364131 +"420",0,6,4.58597459245456 +"421",0,7,4.59869165172741 +"422",1,9,4.61065589394319 +"423",1,10,4.6138947467798 +"424",1,10,4.62103158135052 +"425",1,1,4.64469648165712 +"426",1,5,4.64953727676617 +"427",0,8,4.65244470117698 +"428",1,7,4.70014799998905 +"429",1,5,4.70793099753348 +"430",1,6,4.72796639645178 +"431",1,10,4.72927502950359 +"432",1,7,4.74986681574577 +"433",1,9,4.75290466232313 +"434",0,5,4.76022427043506 +"435",0,7,4.78595927756713 +"436",0,6,4.79543832419832 +"437",1,2,4.79931617277134 +"438",1,2,4.79968980932534 +"439",1,8,4.81564458242504 +"440",1,7,4.84126812831259 +"441",1,5,4.84777773040573 +"442",1,3,4.86859810804117 +"443",1,9,4.86898236614901 +"444",1,8,4.87879169278985 +"445",0,7,4.8888839742961 +"446",1,7,4.8915373496759 +"447",1,9,4.91536076265436 +"448",0,10,4.92551905689768 +"449",0,3,4.94430989850235 +"450",1,1,4.94650637648923 +"451",1,5,4.96040773521487 +"452",0,9,4.96125598259366 +"453",1,3,4.98455381197825 +"454",1,2,4.99415913254882 +"455",1,2,4.99814815998819 +"456",0,3,5.00753996048426 +"457",1,1,5.01482058648948 +"458",0,3,5.01965704577972 +"459",1,2,5.02195174519435 +"460",1,7,5.03237780481786 +"461",0,3,5.03963740799743 +"462",0,2,5.04090728749223 +"463",0,3,5.06982290182728 +"464",1,1,5.07097486058715 +"465",0,6,5.07260015944795 +"466",1,3,5.08396319784959 +"467",1,1,5.10150934536656 +"468",1,10,5.10700652853766 +"469",0,2,5.10924989571509 +"470",1,4,5.11774118590772 +"471",1,3,5.14585436086606 +"472",1,1,5.15656349451901 +"473",1,7,5.1674320587108 +"474",1,3,5.17099888593383 +"475",0,7,5.17557663877919 +"476",1,10,5.18732034485238 +"477",1,3,5.21305349205877 +"478",1,10,5.2195541103629 +"479",0,9,5.23472346899692 +"480",1,8,5.23490989191119 +"481",1,3,5.24936787634665 +"482",0,4,5.27201364443964 +"483",0,2,5.27745435968833 +"484",1,7,5.28242102179071 +"485",1,4,5.28423480070988 +"486",0,1,5.28630745664436 +"487",0,7,5.2896281768867 +"488",0,2,5.37070022581812 +"489",1,7,5.37855917775697 +"490",1,7,5.38528195868692 +"491",1,10,5.39124680714662 +"492",1,1,5.44078192898858 +"493",0,4,5.45730074589842 +"494",1,4,5.47042568596739 +"495",0,8,5.47105650072645 +"496",1,7,5.47169257553022 +"497",0,8,5.49879760004182 +"498",1,4,5.49986094396164 +"499",0,9,5.50072378966862 +"500",0,7,5.50946431542942 +"501",0,2,5.5116423277623 +"502",0,8,5.51637331811906 +"503",0,10,5.54162393920193 +"504",0,6,5.54234474881621 +"505",1,1,5.54446032376163 +"506",0,3,5.55870300148664 +"507",0,4,5.5608766827784 +"508",0,5,5.56195585462337 +"509",0,8,5.56901440484479 +"510",1,1,5.57432649399349 +"511",1,7,5.58622223174352 +"512",1,4,5.58955941164088 +"513",1,9,5.59268920765877 +"514",0,8,5.59305585976194 +"515",1,7,5.61088087414743 +"516",1,1,5.61115534101794 +"517",1,9,5.63304386369008 +"518",1,1,5.63592613955084 +"519",1,2,5.64217747753076 +"520",1,5,5.64989018552484 +"521",1,3,5.66145276907167 +"522",1,5,5.6778854876967 +"523",0,2,5.68335496549735 +"524",1,9,5.69229081797612 +"525",1,2,5.70685198561403 +"526",1,5,5.71278508550137 +"527",1,6,5.72101738121899 +"528",1,2,5.72681474441735 +"529",1,5,5.72798995837176 +"530",1,7,5.74457760406926 +"531",1,3,5.74537408167764 +"532",0,10,5.7484570182002 +"533",0,3,5.76756871911416 +"534",1,3,5.79382520787942 +"535",1,4,5.797074770006 +"536",0,3,5.79752263499482 +"537",1,7,5.81890333551846 +"538",1,10,5.84756833288583 +"539",1,8,5.85040079438823 +"540",0,3,5.8528831589602 +"541",1,3,5.86299406621639 +"542",0,5,5.877209509486 +"543",0,4,5.8984263347562 +"544",1,2,5.90132497815951 +"545",1,7,5.92697740253471 +"546",1,8,5.93228751148041 +"547",1,5,5.93482327178273 +"548",1,1,5.94197185075021 +"549",0,3,5.94856410386151 +"550",1,2,5.95394166126314 +"551",1,7,5.96131596799365 +"552",1,5,5.97345647779607 +"553",1,6,5.98307930524733 +"554",1,3,5.99331972946643 +"555",1,6,6.02685088832925 +"556",0,5,6.04842200553912 +"557",0,3,6.05805513312205 +"558",0,7,6.06738513144718 +"559",1,8,6.06979805426128 +"560",0,1,6.08544456162793 +"561",0,7,6.10219821728054 +"562",0,10,6.10526478683796 +"563",1,10,6.10883280962875 +"564",1,4,6.11260202788917 +"565",1,1,6.11743007708045 +"566",1,9,6.1224602198478 +"567",0,1,6.12348464375531 +"568",1,1,6.12524147150993 +"569",0,8,6.13082185571968 +"570",1,10,6.15106207994615 +"571",1,7,6.15463170201389 +"572",1,6,6.15608513980201 +"573",1,1,6.18178259548633 +"574",1,9,6.19106482802571 +"575",0,8,6.19281789308733 +"576",1,7,6.19764108751754 +"577",1,4,6.20636938798302 +"578",1,8,6.2164933958372 +"579",1,9,6.21886006319914 +"580",0,8,6.22709666074785 +"581",1,9,6.22812848190842 +"582",0,6,6.25214449443794 +"583",1,4,6.25273827375512 +"584",1,2,6.25440602530859 +"585",0,7,6.27871420416057 +"586",1,1,6.28044144613676 +"587",1,1,6.28859874294663 +"588",0,2,6.31953475516885 +"589",1,8,6.32430273963268 +"590",0,2,6.34942667264495 +"591",1,5,6.37757456395611 +"592",1,5,6.39548320967456 +"593",1,9,6.41162784948613 +"594",0,3,6.42685854766102 +"595",0,1,6.4415304446499 +"596",0,6,6.45652572551936 +"597",1,2,6.46267159130153 +"598",1,4,6.47045373955428 +"599",0,1,6.48390834702755 +"600",1,6,6.49101842861225 +"601",1,4,6.50856954717527 +"602",0,7,6.51473787777576 +"603",0,4,6.51977671257204 +"604",1,2,6.53348894396739 +"605",1,4,6.59267601037004 +"606",1,10,6.6001621705168 +"607",1,2,6.61614417335541 +"608",1,5,6.63855812859186 +"609",1,5,6.64666744029504 +"610",1,7,6.66737408371913 +"611",1,7,6.67362010066668 +"612",1,5,6.68247931663175 +"613",1,6,6.68726597153937 +"614",0,8,6.7106974057783 +"615",1,4,6.73104727015976 +"616",0,5,6.73197305687613 +"617",0,6,6.73311444291131 +"618",0,1,6.73906212041372 +"619",1,10,6.745648126148 +"620",0,4,6.74797439384492 +"621",0,3,6.78696551390931 +"622",0,9,6.78791648056831 +"623",1,3,6.79423769738067 +"624",1,2,6.80302944960504 +"625",0,3,6.80828803935985 +"626",0,5,6.80887244042433 +"627",1,3,6.8118559683017 +"628",1,7,6.81409027945674 +"629",0,2,6.83865800081119 +"630",1,4,6.85356852924665 +"631",0,3,6.8721033150902 +"632",0,4,6.87444968239534 +"633",1,3,6.87998226445764 +"634",0,5,6.88000141856789 +"635",0,6,6.88083361907924 +"636",0,5,6.88147788888837 +"637",1,7,6.88830928375541 +"638",1,5,6.89569035520332 +"639",0,10,6.94073649266037 +"640",0,3,6.94193095056108 +"641",1,9,6.9419495675224 +"642",1,10,6.9420818155163 +"643",0,9,6.96800694399359 +"644",1,8,6.97161459137066 +"645",1,4,6.99405301648469 +"646",1,2,7.00686403741363 +"647",0,4,7.02237613837385 +"648",1,6,7.02744361316936 +"649",0,5,7.03949411676374 +"650",1,4,7.03990150285062 +"651",0,2,7.06119778103178 +"652",1,4,7.06761330960767 +"653",0,10,7.09111519489759 +"654",1,5,7.10434780936037 +"655",1,6,7.11022266105802 +"656",0,5,7.11580935192348 +"657",1,1,7.12576509910207 +"658",1,3,7.13509479559488 +"659",1,5,7.13656199326734 +"660",1,3,7.14249064255967 +"661",0,1,7.16062369544436 +"662",0,7,7.18431675328193 +"663",0,1,7.19089744946239 +"664",1,1,7.19973092277024 +"665",1,7,7.20039856318805 +"666",0,9,7.20776582065603 +"667",1,2,7.22842523583692 +"668",0,6,7.23751919283039 +"669",0,8,7.2388265166592 +"670",0,7,7.24698782900213 +"671",0,2,7.24898790252189 +"672",1,9,7.26425906172664 +"673",1,4,7.27905431774071 +"674",1,4,7.28582217853608 +"675",1,9,7.28587325681164 +"676",1,6,7.29212851145468 +"677",1,8,7.2957114282746 +"678",0,10,7.29694643892072 +"679",0,7,7.29923831795998 +"680",0,3,7.30121113799841 +"681",1,9,7.30913644192707 +"682",0,9,7.31256624045748 +"683",1,3,7.32807709602704 +"684",0,1,7.32983511657519 +"685",1,1,7.33044432573239 +"686",1,1,7.33793513369313 +"687",1,3,7.35371518274228 +"688",1,3,7.35432196648716 +"689",1,1,7.35548573603716 +"690",1,10,7.35907829484765 +"691",1,9,7.36456607312292 +"692",1,8,7.36494202421841 +"693",1,8,7.37015086711247 +"694",0,5,7.38384602965608 +"695",1,4,7.38688083084586 +"696",1,2,7.38826103712354 +"697",1,8,7.39294512937016 +"698",0,8,7.39652476176325 +"699",1,6,7.40652461316279 +"700",1,1,7.41579990515286 +"701",1,10,7.42819672796634 +"702",1,2,7.4340501835501 +"703",0,4,7.44313182033563 +"704",1,2,7.45351504685279 +"705",1,7,7.49311339685601 +"706",1,5,7.49814938984719 +"707",1,9,7.53064154550173 +"708",0,1,7.55110260080202 +"709",0,5,7.59194077537772 +"710",1,7,7.59986609291444 +"711",1,2,7.60323896488333 +"712",1,2,7.6351457699201 +"713",0,8,7.64260643282083 +"714",1,6,7.65334569466368 +"715",0,5,7.67424731559061 +"716",1,2,7.6842571544544 +"717",0,10,7.69295039083666 +"718",0,7,7.70137112256073 +"719",1,6,7.70179872715849 +"720",0,6,7.70481276080015 +"721",1,8,7.71332770438826 +"722",0,9,7.71638329250233 +"723",1,4,7.72841194940927 +"724",0,7,7.73881301434865 +"725",1,8,7.75243180896267 +"726",0,2,7.76485023474391 +"727",0,7,7.77310304357698 +"728",1,8,7.78676063444451 +"729",0,4,7.78842926426184 +"730",1,4,7.79212278347787 +"731",1,3,7.802399852437 +"732",1,5,7.81769934282124 +"733",1,4,7.84816053688776 +"734",0,2,7.85866216698869 +"735",0,8,7.87339243988971 +"736",0,3,7.88558913270434 +"737",1,6,7.88878101440104 +"738",0,3,7.91623116951267 +"739",1,5,7.92460416875569 +"740",0,5,7.92945862720044 +"741",1,9,7.93295615901033 +"742",0,7,7.93436945867519 +"743",1,8,7.93553499740686 +"744",1,7,7.94160670027275 +"745",0,9,7.94520362720963 +"746",1,6,7.96527626563849 +"747",1,9,7.9720735023083 +"748",1,3,7.97309516912846 +"749",1,8,7.97816455768475 +"750",1,3,7.99002325641935 +"751",1,6,7.9917673605542 +"752",1,2,7.99180043902412 +"753",1,7,8.02830900709001 +"754",0,7,8.0488980129804 +"755",0,10,8.05455903667729 +"756",1,10,8.06000253519281 +"757",1,1,8.09272471277113 +"758",1,9,8.09508041644045 +"759",0,5,8.11033549380912 +"760",1,2,8.11110705778426 +"761",0,9,8.11483718105084 +"762",1,10,8.11484695630572 +"763",1,9,8.13703728385567 +"764",1,6,8.15976403233874 +"765",0,8,8.18974166837631 +"766",1,3,8.20129380571724 +"767",1,1,8.20149816728044 +"768",0,2,8.20239239251554 +"769",0,1,8.20689482922513 +"770",1,8,8.21071226585503 +"771",0,9,8.21571429477386 +"772",0,9,8.24300839479268 +"773",1,4,8.2567555190451 +"774",1,2,8.2613302847059 +"775",1,1,8.26202060470037 +"776",1,10,8.26730947267409 +"777",1,7,8.27590112022231 +"778",1,5,8.31760680715408 +"779",1,5,8.34371126257543 +"780",1,3,8.35247169620044 +"781",0,2,8.36016904043172 +"782",1,9,8.38953395659074 +"783",0,3,8.40166611176844 +"784",1,1,8.40927771135119 +"785",1,5,8.41988254886896 +"786",1,2,8.4255096829191 +"787",1,5,8.45237879392306 +"788",0,8,8.45420366263011 +"789",1,9,8.47479141877874 +"790",1,6,8.49209979772394 +"791",1,3,8.50493261661165 +"792",0,8,8.51114865656261 +"793",0,10,8.53287157503081 +"794",1,4,8.53468113443052 +"795",1,4,8.54588873100594 +"796",0,1,8.54591401747555 +"797",0,1,8.55759722240354 +"798",1,4,8.55815653445284 +"799",1,4,8.57573886204512 +"800",1,5,8.57614736296238 +"801",0,9,8.57820050904063 +"802",1,1,8.58127449082678 +"803",1,9,8.5835402505231 +"804",1,2,8.59058347728442 +"805",1,7,8.59062162067665 +"806",1,3,8.5919189298113 +"807",1,4,8.59484634127069 +"808",1,3,8.60417502964406 +"809",1,8,8.60522835392423 +"810",0,4,8.60865239795578 +"811",1,6,8.61154498080072 +"812",1,8,8.61288802904357 +"813",0,1,8.64194058326185 +"814",1,2,8.64494620737402 +"815",0,7,8.64641524506094 +"816",1,2,8.64790682893159 +"817",0,9,8.66148712254 +"818",1,7,8.66363743162119 +"819",0,5,8.67071762902957 +"820",0,9,8.67563223561567 +"821",1,8,8.6959577761239 +"822",1,4,8.70251732417127 +"823",1,10,8.72722072339236 +"824",1,9,8.73415120061106 +"825",1,10,8.73502901138602 +"826",1,9,8.76997049223175 +"827",0,2,8.7945349207276 +"828",1,6,8.8251758919395 +"829",1,1,8.82823498157965 +"830",0,8,8.83538109136884 +"831",1,1,8.83896594216255 +"832",1,6,8.84938170900624 +"833",1,7,8.85312843483462 +"834",1,10,8.8579191642912 +"835",0,1,8.85985264699641 +"836",1,6,8.86253942630501 +"837",1,7,8.86720015853891 +"838",0,8,8.88550578370854 +"839",0,2,8.91482346743341 +"840",1,5,8.93130797209511 +"841",1,10,8.96753994060481 +"842",0,4,8.96836158128855 +"843",1,10,8.96851365589583 +"844",0,3,8.97255537239415 +"845",0,4,8.97972660676114 +"846",1,7,8.9806762067648 +"847",1,2,8.98094284660618 +"848",1,5,8.98573019597649 +"849",1,10,8.99890482271017 +"850",1,7,9.00161035315061 +"851",1,9,9.0210382980372 +"852",1,2,9.0390046639437 +"853",1,9,9.04129916140569 +"854",1,1,9.06826070423802 +"855",0,9,9.07408435050497 +"856",1,7,9.07544092012727 +"857",1,4,9.07648413614527 +"858",1,4,9.08505513360839 +"859",0,5,9.09449844049363 +"860",1,3,9.10410477135362 +"861",0,5,9.10432989137279 +"862",0,7,9.11420001370822 +"863",1,9,9.12571458920784 +"864",1,10,9.1267388411187 +"865",1,2,9.12758909810866 +"866",1,5,9.13242605614565 +"867",1,8,9.13406300553396 +"868",1,3,9.14240695479559 +"869",1,2,9.15002997824182 +"870",1,8,9.16794751806131 +"871",1,6,9.17525471810353 +"872",1,6,9.17548135139615 +"873",0,7,9.17574581838774 +"874",0,1,9.18325567181697 +"875",1,7,9.19189149021753 +"876",1,7,9.19821243540035 +"877",1,1,9.22576128067905 +"878",0,4,9.26723323308178 +"879",0,2,9.29506858034534 +"880",1,4,9.31979224831214 +"881",1,4,9.35575416555767 +"882",1,7,9.36474756295889 +"883",1,7,9.37011989284837 +"884",1,5,9.38731487698061 +"885",1,1,9.4068545898335 +"886",1,8,9.43361295220371 +"887",1,7,9.4392533983935 +"888",0,2,9.45114846090673 +"889",0,1,9.4678811739086 +"890",1,4,9.4918931696066 +"891",0,3,9.50049630269975 +"892",1,1,9.50151338700854 +"893",1,9,9.5159290230517 +"894",1,4,9.53755662519851 +"895",1,3,9.53765222100442 +"896",1,10,9.55905850112537 +"897",1,9,9.56893221415636 +"898",1,10,9.57688587835249 +"899",0,2,9.59062742667421 +"900",0,8,9.61938493663426 +"901",1,1,9.65043530482793 +"902",0,3,9.65560504808446 +"903",1,7,9.65752741710915 +"904",0,7,9.66097610443097 +"905",1,1,9.66403159831382 +"906",1,3,9.67550307267967 +"907",1,5,9.67615987705222 +"908",1,4,9.68887278152743 +"909",1,7,9.70440249672189 +"910",1,7,9.71140695401455 +"911",1,4,9.71896299670396 +"912",1,1,9.72042907098254 +"913",0,1,9.72498445094509 +"914",0,9,9.72701336723698 +"915",1,9,9.73122632501635 +"916",1,2,9.75069152094328 +"917",1,7,9.75352854653166 +"918",0,4,9.76223946453866 +"919",0,3,9.76784596049555 +"920",1,3,9.77305446222753 +"921",0,8,9.79819145543703 +"922",1,1,9.80480159938537 +"923",1,10,9.8153209762105 +"924",1,8,9.81810006885018 +"925",1,6,9.83626556178631 +"926",0,9,9.84640568409149 +"927",1,8,9.85896251406512 +"928",1,9,9.86044981235246 +"929",0,5,9.8608584066674 +"930",1,4,9.86776858321464 +"931",1,8,9.87061278298211 +"932",1,10,9.88927142504233 +"933",0,2,9.89533092475126 +"934",1,10,9.90555509947839 +"935",1,9,9.9295750849555 +"936",1,3,9.96890765022749 +"937",1,3,10.007791426648 +"938",1,6,10.0090428372292 +"939",0,8,10.0093533862276 +"940",1,3,10.0128024030556 +"941",1,6,10.0136952140307 +"942",1,2,10.0219068452743 +"943",1,10,10.030473405479 +"944",1,7,10.037067654503 +"945",1,3,10.0528349269522 +"946",1,7,10.0571804667465 +"947",1,3,10.0651593441992 +"948",1,6,10.080455521485 +"949",1,7,10.0809592146568 +"950",1,9,10.0851293660097 +"951",0,8,10.0968603472456 +"952",1,3,10.1013688859434 +"953",1,4,10.1033175175481 +"954",1,1,10.1049650241092 +"955",1,4,10.10883281892 +"956",1,2,10.1127286610084 +"957",1,5,10.1644177063043 +"958",1,5,10.1790396966079 +"959",0,4,10.1937908904462 +"960",1,10,10.2074489778111 +"961",0,2,10.2501976287202 +"962",1,5,10.2544440464745 +"963",0,9,10.2573297115345 +"964",1,8,10.2575774971787 +"965",1,2,10.2627113002278 +"966",1,6,10.2848001217755 +"967",0,5,10.2885358053175 +"968",1,8,10.292939995286 +"969",1,2,10.2978037796131 +"970",0,5,10.3036457540649 +"971",0,4,10.3037603079745 +"972",1,8,10.3236747295089 +"973",1,5,10.3351399873897 +"974",1,4,10.350238087201 +"975",1,8,10.3523402127209 +"976",1,9,10.3547170156654 +"977",1,3,10.3781587761855 +"978",1,8,10.3812371285101 +"979",1,3,10.3922689508764 +"980",1,5,10.4040469284096 +"981",0,4,10.4059109058949 +"982",1,1,10.4076424896696 +"983",1,6,10.4147728065648 +"984",0,5,10.4291132478216 +"985",1,5,10.4335531722787 +"986",1,8,10.4396352795311 +"987",0,8,10.4541859170772 +"988",1,10,10.4754108621704 +"989",0,2,10.4763470037582 +"990",1,3,10.4765305347247 +"991",1,7,10.517204408948 +"992",1,4,10.524695387617 +"993",1,2,10.5305990086955 +"994",1,5,10.5391673866347 +"995",0,6,10.561788420179 +"996",1,8,10.5644669071679 +"997",0,2,10.565648077191 +"998",1,2,10.5822130709179 +"999",1,9,10.6063243762615 +"1000",1,6,10.6151280138103 diff --git a/test/runtests.jl b/test/runtests.jl index d8bab6f..67b2626 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using StableRNGs using Random using Distributions using DataFrames +using CSV using Turing using SafetySignalDetection diff --git a/test/small_historic_dataset.csv b/test/small_historic_dataset.csv new file mode 100644 index 0000000..bf51bee --- /dev/null +++ b/test/small_historic_dataset.csv @@ -0,0 +1,101 @@ +"","y","trial","time" +"1",0,5,0.508060037142824 +"2",0,4,0.5208363096442 +"3",0,2,0.556331578741567 +"4",0,4,0.563791380921675 +"5",1,2,0.571754800768225 +"6",0,5,0.590052150745623 +"7",0,5,0.591040918682293 +"8",0,2,0.595221235355885 +"9",0,4,0.597784269735739 +"10",1,2,0.607071721277385 +"11",0,1,0.611330809521675 +"12",0,4,0.613947846980951 +"13",0,4,0.618876835716143 +"14",0,2,0.640354126490284 +"15",0,3,0.645504716144507 +"16",0,3,0.648776558221881 +"17",0,4,0.653263780100834 +"18",0,1,0.664869826653364 +"19",0,4,0.671171271026457 +"20",1,2,0.674994900783847 +"21",0,3,0.679323147659513 +"22",0,3,0.681681905904048 +"23",0,2,0.693983536523541 +"24",0,4,0.713401719711471 +"25",0,3,0.715330274447549 +"26",0,5,0.715554267597332 +"27",0,5,0.724320719842258 +"28",0,3,0.743751407473829 +"29",0,1,0.744621157756557 +"30",0,5,0.749986377628643 +"31",0,2,0.756438894458238 +"32",0,3,0.769038580022455 +"33",0,5,0.781380031033633 +"34",0,3,0.782826549365607 +"35",0,3,0.78411009277818 +"36",0,2,0.791492904824573 +"37",1,4,0.801637802299991 +"38",0,3,0.811791839104605 +"39",0,5,0.837362781844232 +"40",0,3,0.841242652393076 +"41",0,5,0.862928944447668 +"42",1,2,0.869636503880502 +"43",0,2,0.886621511367525 +"44",0,1,0.889937450386028 +"45",0,2,0.89081983407345 +"46",0,4,0.894065672679249 +"47",0,2,0.895545855571809 +"48",0,4,0.909240301239177 +"49",0,3,0.919524707524235 +"50",0,4,0.924650129098574 +"51",0,1,0.945304040451132 +"52",0,5,0.946978089384413 +"53",1,3,0.951859002214535 +"54",0,4,0.957347076003118 +"55",0,1,0.960843806894175 +"56",0,3,0.965779375707409 +"57",1,5,0.966733577783565 +"58",0,1,0.96968774609861 +"59",0,1,0.972657807600261 +"60",1,4,0.979699684335007 +"61",0,4,0.979834431354125 +"62",0,1,0.990429732491942 +"63",0,4,0.993043940746905 +"64",0,5,0.9954099067624 +"65",0,3,0.997159593415083 +"66",0,5,1.00392896736284 +"67",0,3,1.01355110456085 +"68",0,2,1.04904199831044 +"69",1,3,1.06484850264858 +"70",0,4,1.06794773782012 +"71",0,3,1.0712293948771 +"72",0,2,1.08236479130322 +"73",0,4,1.0874704878659 +"74",0,3,1.08815656886374 +"75",0,4,1.09446561418411 +"76",0,4,1.11639969562923 +"77",1,5,1.13149712565995 +"78",0,4,1.15531456740156 +"79",0,2,1.17462593214085 +"80",0,1,1.18276276532017 +"81",1,1,1.18776627686728 +"82",0,4,1.20208448878763 +"83",0,1,1.21019494540132 +"84",0,2,1.21124169832635 +"85",0,1,1.2113759552136 +"86",1,2,1.21897941682307 +"87",0,4,1.22070786681776 +"88",0,2,1.24966304313974 +"89",1,4,1.2527988999829 +"90",0,5,1.25518082340978 +"91",0,3,1.25891242750253 +"92",0,5,1.26357294799558 +"93",0,2,1.26469715973161 +"94",0,4,1.27593270482697 +"95",0,5,1.28273087209758 +"96",1,2,1.29175769180613 +"97",0,3,1.31240415124692 +"98",1,5,1.31259288685395 +"99",0,4,1.33657055292152 +"100",0,2,1.34760702875274 diff --git a/test/test_fit_beta_mixture.jl b/test/test_fit_beta_mixture.jl index 137add2..a68c100 100644 --- a/test/test_fit_beta_mixture.jl +++ b/test/test_fit_beta_mixture.jl @@ -47,3 +47,57 @@ end @test mix_fit.prior.p[1] ≈ 1 - π rtol = 0.01 @test mix_fit.prior.p[2] ≈ π rtol = 0.01 end + +@testset "Reconcile fit_beta_mixture on known datasets" begin + # Create mixtures from MAP priors on known datasets + + seed = 2024 + prior_a = Beta(1 / 3, 1 / 3) + prior_b = Beta(5, 5) + prob_ctrl = 0.333 + + # Small historical dataset + df_small = DataFrame(CSV.File("small_historic_dataset.csv")) + df_small.y = map(x -> x == 1, df_small.y) # Cast to Vector{Bool} + rng = StableRNG(123) + map_small = sample( + rng, + meta_analytic(df_small.y, df_small.time, df_small.trial, + prior_a, prior_b), + HMC(0.05, 10), + 10_000 + ) + df_map_small = DataFrame(map_small) + pi_star_alpha = df_map_small.a .* df_map_small.b * nrow(df_small) + pi_star_beta = (1 .- df_map_small.a) .* df_map_small.b * nrow(df_small) + pi_star = [rand(Beta(alpha, beta), 1) + for (alpha, beta) in zip(pi_star_alpha, pi_star_beta)] + pi_star = collect(Iterators.flatten(pi_star)) # Flatten + mix_small = fit_beta_mixture(pi_star, 2) + @test mean(mix_small) ≈ mean(pi_star) atol = 0.01 + @test std(mix_small) ≈ std(pi_star) atol = 0.01 + + + # Large historical dataset + df_large = DataFrame(CSV.File("large_historic_dataset.csv")) + df_large.y = map(x -> x == 1, df_large.y) # Cast to Vector{Bool} + + rng = StableRNG(123) + map_large = sample( + rng, + meta_analytic(df_large.y, df_large.time, df_large.trial, + prior_a, prior_b), + HMC(0.05, 10), + 10_000 + ) + df_map_large = DataFrame(map_large) + pi_star_alpha = df_map_small.a .* df_map_small.b * nrow(df_large) + pi_star_beta = (1 .- df_map_small.a) .* df_map_small.b * nrow(df_large) + pi_star = [rand(Beta(alpha, beta), 1) + for (alpha, beta) in zip(pi_star_alpha, pi_star_beta)] + pi_star = collect(Iterators.flatten(pi_star)) # Flatten + mix_large = fit_beta_mixture(pi_star, 2) + @test mean(mix_large) ≈ mean(pi_star) atol = 0.01 + @test std(mix_large) ≈ std(pi_star) atol = 0.01 + +end \ No newline at end of file diff --git a/test/test_meta_analytic.jl b/test/test_meta_analytic.jl index 4674f86..cfe9703 100644 --- a/test/test_meta_analytic.jl +++ b/test/test_meta_analytic.jl @@ -19,3 +19,42 @@ check_numerical(chain, [:a], [0.223], rtol=0.001) check_numerical(chain, [:b], [0.485], rtol=0.001) end + +@testset "Reconcile meta_analytic on known datasets with rstan" begin + # Create MAP priors on known datasets and compare to rstan output + + seed = 2024 + prior_a = Beta(1 / 3, 1 / 3) + prior_b = Beta(5, 5) + prob_ctrl = 0.333 + + # Small historical dataset + df_small = DataFrame(CSV.File("small_historic_dataset.csv")) + df_small.y = map(x -> x == 1, df_small.y) # Cast to Vector{Bool} + rng = StableRNG(123) + map_small = sample( + rng, + meta_analytic(df_small.y, df_small.time, df_small.trial, + prior_a, prior_b), + HMC(0.05, 10), + 10_000 + ) + check_numerical(map_small, [:a], [0.17], atol = 0.01) + check_numerical(map_small, [:b], [0.51], atol = 0.01) + + + # Large historical dataset + df_large = DataFrame(CSV.File("large_historic_dataset.csv")) + df_large.y = map(x -> x == 1, df_large.y) # Cast to Vector{Bool} + rng = StableRNG(123) + map_large = sample( + rng, + meta_analytic(df_large.y, df_large.time, df_large.trial, + prior_a, prior_b), + HMC(0.05, 10), + 10_000 + ) + check_numerical(map_large, [:a], [0.13], atol = 0.01) + check_numerical(map_large, [:b], [0.55], atol = 0.01) + +end \ No newline at end of file From c783e54178ca987d872464fb4730574b6919d78e Mon Sep 17 00:00:00 2001 From: Kristian Brock Date: Thu, 11 Jan 2024 11:20:00 +0000 Subject: [PATCH 4/8] Split similar tests, and rolled back Julia version in Pluto nb --- design/ReconcileWithR.jl | 2 +- test/test_fit_beta_mixture.jl | 21 +++++++++++++-------- test/test_meta_analytic.jl | 17 +++++++++++------ 3 files changed, 25 insertions(+), 15 deletions(-) diff --git a/design/ReconcileWithR.jl b/design/ReconcileWithR.jl index e599846..43bbf9f 100644 --- a/design/ReconcileWithR.jl +++ b/design/ReconcileWithR.jl @@ -335,7 +335,7 @@ Turing = "~0.29.3" PLUTO_MANIFEST_TOML_CONTENTS = """ # This file is machine-generated - editing it directly is not advised -julia_version = "1.9.1" +julia_version = "1.9.3" manifest_format = "2.0" project_hash = "4359e7ffd037ffcea031434e2d83a0afe89b7541" diff --git a/test/test_fit_beta_mixture.jl b/test/test_fit_beta_mixture.jl index a68c100..779f5b0 100644 --- a/test/test_fit_beta_mixture.jl +++ b/test/test_fit_beta_mixture.jl @@ -48,15 +48,13 @@ end @test mix_fit.prior.p[2] ≈ π rtol = 0.01 end -@testset "Reconcile fit_beta_mixture on known datasets" begin - # Create mixtures from MAP priors on known datasets +@testset "Reconcile fit_beta_mixture on small historical dataset" begin + # Create mixture from MAP prior on known small historical dataset - seed = 2024 prior_a = Beta(1 / 3, 1 / 3) prior_b = Beta(5, 5) prob_ctrl = 0.333 - # Small historical dataset df_small = DataFrame(CSV.File("small_historic_dataset.csv")) df_small.y = map(x -> x == 1, df_small.y) # Cast to Vector{Bool} rng = StableRNG(123) @@ -76,9 +74,16 @@ end mix_small = fit_beta_mixture(pi_star, 2) @test mean(mix_small) ≈ mean(pi_star) atol = 0.01 @test std(mix_small) ≈ std(pi_star) atol = 0.01 + +end + +@testset "Reconcile fit_beta_mixture on large historical dataset" begin + # Create mixture from MAP prior on known large historical dataset + prior_a = Beta(1 / 3, 1 / 3) + prior_b = Beta(5, 5) + prob_ctrl = 0.333 - # Large historical dataset df_large = DataFrame(CSV.File("large_historic_dataset.csv")) df_large.y = map(x -> x == 1, df_large.y) # Cast to Vector{Bool} @@ -91,8 +96,8 @@ end 10_000 ) df_map_large = DataFrame(map_large) - pi_star_alpha = df_map_small.a .* df_map_small.b * nrow(df_large) - pi_star_beta = (1 .- df_map_small.a) .* df_map_small.b * nrow(df_large) + pi_star_alpha = df_map_large.a .* df_map_large.b * nrow(df_large) + pi_star_beta = (1 .- df_map_large.a) .* df_map_large.b * nrow(df_large) pi_star = [rand(Beta(alpha, beta), 1) for (alpha, beta) in zip(pi_star_alpha, pi_star_beta)] pi_star = collect(Iterators.flatten(pi_star)) # Flatten @@ -100,4 +105,4 @@ end @test mean(mix_large) ≈ mean(pi_star) atol = 0.01 @test std(mix_large) ≈ std(pi_star) atol = 0.01 -end \ No newline at end of file +end diff --git a/test/test_meta_analytic.jl b/test/test_meta_analytic.jl index cfe9703..d4a098c 100644 --- a/test/test_meta_analytic.jl +++ b/test/test_meta_analytic.jl @@ -20,15 +20,13 @@ check_numerical(chain, [:b], [0.485], rtol=0.001) end -@testset "Reconcile meta_analytic on known datasets with rstan" begin - # Create MAP priors on known datasets and compare to rstan output +@testset "Reconcile meta_analytic with rstan on small historical dataset" begin + # Create MAP priors on known small historical dataset and compare to rstan - seed = 2024 prior_a = Beta(1 / 3, 1 / 3) prior_b = Beta(5, 5) prob_ctrl = 0.333 - # Small historical dataset df_small = DataFrame(CSV.File("small_historic_dataset.csv")) df_small.y = map(x -> x == 1, df_small.y) # Cast to Vector{Bool} rng = StableRNG(123) @@ -41,9 +39,16 @@ end ) check_numerical(map_small, [:a], [0.17], atol = 0.01) check_numerical(map_small, [:b], [0.51], atol = 0.01) - - # Large historical dataset +end + +@testset "Reconcile meta_analytic with rstan on large historical dataset" begin + # Create MAP priors on known large historical dataset and compare to rstan + + prior_a = Beta(1 / 3, 1 / 3) + prior_b = Beta(5, 5) + prob_ctrl = 0.333 + df_large = DataFrame(CSV.File("large_historic_dataset.csv")) df_large.y = map(x -> x == 1, df_large.y) # Cast to Vector{Bool} rng = StableRNG(123) From 0517179f5e51228828d73bcb16d4a9fb771699ac Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 11 Jan 2024 14:52:05 +0100 Subject: [PATCH 5/8] use NUTS for better samples in tests, and rtol to be stricter --- test/test_fit_beta_mixture.jl | 12 ++++++------ test/test_meta_analytic.jl | 12 ++++++------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/test/test_fit_beta_mixture.jl b/test/test_fit_beta_mixture.jl index 779f5b0..ae6f230 100644 --- a/test/test_fit_beta_mixture.jl +++ b/test/test_fit_beta_mixture.jl @@ -62,7 +62,7 @@ end rng, meta_analytic(df_small.y, df_small.time, df_small.trial, prior_a, prior_b), - HMC(0.05, 10), + NUTS(0.65), 10_000 ) df_map_small = DataFrame(map_small) @@ -72,8 +72,8 @@ end for (alpha, beta) in zip(pi_star_alpha, pi_star_beta)] pi_star = collect(Iterators.flatten(pi_star)) # Flatten mix_small = fit_beta_mixture(pi_star, 2) - @test mean(mix_small) ≈ mean(pi_star) atol = 0.01 - @test std(mix_small) ≈ std(pi_star) atol = 0.01 + @test mean(mix_small) ≈ mean(pi_star) rtol = 0.01 + @test std(mix_small) ≈ std(pi_star) rtol = 0.01 end @@ -92,7 +92,7 @@ end rng, meta_analytic(df_large.y, df_large.time, df_large.trial, prior_a, prior_b), - HMC(0.05, 10), + NUTS(0.65), 10_000 ) df_map_large = DataFrame(map_large) @@ -102,7 +102,7 @@ end for (alpha, beta) in zip(pi_star_alpha, pi_star_beta)] pi_star = collect(Iterators.flatten(pi_star)) # Flatten mix_large = fit_beta_mixture(pi_star, 2) - @test mean(mix_large) ≈ mean(pi_star) atol = 0.01 - @test std(mix_large) ≈ std(pi_star) atol = 0.01 + @test mean(mix_large) ≈ mean(pi_star) rtol = 0.01 + @test std(mix_large) ≈ std(pi_star) rtol = 0.01 end diff --git a/test/test_meta_analytic.jl b/test/test_meta_analytic.jl index d4a098c..28d1d63 100644 --- a/test/test_meta_analytic.jl +++ b/test/test_meta_analytic.jl @@ -34,11 +34,11 @@ end rng, meta_analytic(df_small.y, df_small.time, df_small.trial, prior_a, prior_b), - HMC(0.05, 10), + NUTS(0.65), 10_000 ) - check_numerical(map_small, [:a], [0.17], atol = 0.01) - check_numerical(map_small, [:b], [0.51], atol = 0.01) + check_numerical(map_small, [:a], [0.17], rtol = 0.01) + check_numerical(map_small, [:b], [0.51], rtol = 0.01) end @@ -56,10 +56,10 @@ end rng, meta_analytic(df_large.y, df_large.time, df_large.trial, prior_a, prior_b), - HMC(0.05, 10), + NUTS(0.65), 10_000 ) - check_numerical(map_large, [:a], [0.13], atol = 0.01) - check_numerical(map_large, [:b], [0.55], atol = 0.01) + check_numerical(map_large, [:a], [0.13], rtol = 0.01) + check_numerical(map_large, [:b], [0.55], rtol = 0.01) end \ No newline at end of file From 2591b528486440af00c419f7116df7075ac4d4b8 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 11 Jan 2024 15:11:57 +0100 Subject: [PATCH 6/8] ok that was maybe a bit too strict :) --- test/test_fit_beta_mixture.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_fit_beta_mixture.jl b/test/test_fit_beta_mixture.jl index ae6f230..827853a 100644 --- a/test/test_fit_beta_mixture.jl +++ b/test/test_fit_beta_mixture.jl @@ -72,8 +72,8 @@ end for (alpha, beta) in zip(pi_star_alpha, pi_star_beta)] pi_star = collect(Iterators.flatten(pi_star)) # Flatten mix_small = fit_beta_mixture(pi_star, 2) - @test mean(mix_small) ≈ mean(pi_star) rtol = 0.01 - @test std(mix_small) ≈ std(pi_star) rtol = 0.01 + @test mean(mix_small) ≈ mean(pi_star) rtol = 0.02 + @test std(mix_small) ≈ std(pi_star) rtol = 0.02 end @@ -102,7 +102,7 @@ end for (alpha, beta) in zip(pi_star_alpha, pi_star_beta)] pi_star = collect(Iterators.flatten(pi_star)) # Flatten mix_large = fit_beta_mixture(pi_star, 2) - @test mean(mix_large) ≈ mean(pi_star) rtol = 0.01 - @test std(mix_large) ≈ std(pi_star) rtol = 0.01 + @test mean(mix_large) ≈ mean(pi_star) rtol = 0.02 + @test std(mix_large) ≈ std(pi_star) rtol = 0.02 end From 376589b76356e2052009920be4d97554a1253da6 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 2 Feb 2024 22:09:38 +0100 Subject: [PATCH 7/8] simplify tests further --- test/test_fit_beta_mixture.jl | 66 +++++------------------------------ 1 file changed, 9 insertions(+), 57 deletions(-) diff --git a/test/test_fit_beta_mixture.jl b/test/test_fit_beta_mixture.jl index 827853a..d59ee0b 100644 --- a/test/test_fit_beta_mixture.jl +++ b/test/test_fit_beta_mixture.jl @@ -48,61 +48,13 @@ end @test mix_fit.prior.p[2] ≈ π rtol = 0.01 end -@testset "Reconcile fit_beta_mixture on small historical dataset" begin - # Create mixture from MAP prior on known small historical dataset - - prior_a = Beta(1 / 3, 1 / 3) - prior_b = Beta(5, 5) - prob_ctrl = 0.333 - - df_small = DataFrame(CSV.File("small_historic_dataset.csv")) - df_small.y = map(x -> x == 1, df_small.y) # Cast to Vector{Bool} - rng = StableRNG(123) - map_small = sample( - rng, - meta_analytic(df_small.y, df_small.time, df_small.trial, - prior_a, prior_b), - NUTS(0.65), - 10_000 - ) - df_map_small = DataFrame(map_small) - pi_star_alpha = df_map_small.a .* df_map_small.b * nrow(df_small) - pi_star_beta = (1 .- df_map_small.a) .* df_map_small.b * nrow(df_small) - pi_star = [rand(Beta(alpha, beta), 1) - for (alpha, beta) in zip(pi_star_alpha, pi_star_beta)] - pi_star = collect(Iterators.flatten(pi_star)) # Flatten - mix_small = fit_beta_mixture(pi_star, 2) - @test mean(mix_small) ≈ mean(pi_star) rtol = 0.02 - @test std(mix_small) ≈ std(pi_star) rtol = 0.02 - -end - -@testset "Reconcile fit_beta_mixture on large historical dataset" begin - # Create mixture from MAP prior on known large historical dataset - - prior_a = Beta(1 / 3, 1 / 3) - prior_b = Beta(5, 5) - prob_ctrl = 0.333 - - df_large = DataFrame(CSV.File("large_historic_dataset.csv")) - df_large.y = map(x -> x == 1, df_large.y) # Cast to Vector{Bool} - - rng = StableRNG(123) - map_large = sample( - rng, - meta_analytic(df_large.y, df_large.time, df_large.trial, - prior_a, prior_b), - NUTS(0.65), - 10_000 - ) - df_map_large = DataFrame(map_large) - pi_star_alpha = df_map_large.a .* df_map_large.b * nrow(df_large) - pi_star_beta = (1 .- df_map_large.a) .* df_map_large.b * nrow(df_large) - pi_star = [rand(Beta(alpha, beta), 1) - for (alpha, beta) in zip(pi_star_alpha, pi_star_beta)] - pi_star = collect(Iterators.flatten(pi_star)) # Flatten - mix_large = fit_beta_mixture(pi_star, 2) - @test mean(mix_large) ≈ mean(pi_star) rtol = 0.02 - @test std(mix_large) ≈ std(pi_star) rtol = 0.02 - +@testset "Check that fit_beta_mixture can fit the mean and standard deviation of samples well" begin + Random.seed!(123) + N = 100 + normal_dist = Normal(0, 5) + x = rand(normal_dist, N) + y = exp.(x) / (1.0 .+ exp.(x)) + fit_y = fit_beta_mixture(y, 2) + @test mean(fit_y) ≈ mean(y) rtol = 0.02 + @test std(fit_y) ≈ std(y) rtol = 0.02 end From 68670fb7dc2f0c785622e0894dba8d7b96573154 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 2 Feb 2024 22:17:57 +0100 Subject: [PATCH 8/8] and even simpler --- test/test_fit_beta_mixture.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/test_fit_beta_mixture.jl b/test/test_fit_beta_mixture.jl index d59ee0b..b0dff31 100644 --- a/test/test_fit_beta_mixture.jl +++ b/test/test_fit_beta_mixture.jl @@ -51,10 +51,9 @@ end @testset "Check that fit_beta_mixture can fit the mean and standard deviation of samples well" begin Random.seed!(123) N = 100 - normal_dist = Normal(0, 5) - x = rand(normal_dist, N) - y = exp.(x) / (1.0 .+ exp.(x)) + normal_dist = Normal(0.5, 0.05) + y = rand(normal_dist, N) fit_y = fit_beta_mixture(y, 2) @test mean(fit_y) ≈ mean(y) rtol = 0.02 - @test std(fit_y) ≈ std(y) rtol = 0.02 + @test std(fit_y) ≈ std(y) rtol = 0.1 end