From 17a150ed909262b870351489aea023fd7b61bddd Mon Sep 17 00:00:00 2001 From: Steven Alfonso Date: Sun, 18 Feb 2024 21:32:40 -0500 Subject: [PATCH] Remove dependencies and add sampler --- LICENSE | 2 +- Project.toml | 14 +---- README.md | 7 ++- src/McmcHermes.jl | 147 ++++++++++++++++++++++++++++++++++++---------- 4 files changed, 124 insertions(+), 46 deletions(-) diff --git a/LICENSE b/LICENSE index 7cb768a..c1750cd 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2023 stevenalfonso and contributors +Copyright (c) 2023 stevenalfonso and contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/Project.toml b/Project.toml index 6b6fbd5..966bde7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,22 +1,14 @@ name = "McmcHermes" uuid = "1dcdc7e9-a8b1-4892-8a0a-cdab5a60ee98" -authors = ["stevenalfonso "] -version = "1.0.0" +authors = ["stevenalfonso "] +version = "1.0.1" [deps] Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" -DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" -LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] -DataFrames = "1.5" Distributions = "0.25" -LaTeXStrings = "1.3" -Plots = "1.38" -ProgressMeter = "1.7" -julia = "1.8" +julia = "1.10" Compat = "4.6" \ No newline at end of file diff --git a/README.md b/README.md index 667eef5..a269381 100644 --- a/README.md +++ b/README.md @@ -19,15 +19,18 @@ run_mcmc: run multiple chains with a specific number of walkers. get\_flat\_chain: get the stored chain of MCMC samples. -get\_gelman\_rubin: get the Gelman Rubin convergence diagnostic of the chains. +gelman\_rubin\_diagnostic: get the Gelman Rubin convergence diagnostic of the chains. -## Pkg Registry +## Installation ```julia using Pkg Pkg.add("McmcHermes") ``` +## Example + + See here the [documentation](https://stevenalfonso.github.io/McmcHermes.jl/dev/). \ No newline at end of file diff --git a/src/McmcHermes.jl b/src/McmcHermes.jl index 6ad0456..fe1df8c 100644 --- a/src/McmcHermes.jl +++ b/src/McmcHermes.jl @@ -1,85 +1,119 @@ module McmcHermes -using ProgressMeter using Distributions -export one_mcmc, run_mcmc, get_flat_chain, get_gelman_rubin +export one_mcmc, run_mcmc, get_flat_chain, gelman_rubin_diagnostic, sampler + + """ - one_mcmc(log_prob, data, parameters, n_iter, a) +one_mcmc(log_prob, data, parameters, n_iter, a) -Returns one chain with the dimension of (n_iter, 1, n_dim) +Returns one chain with dimension of (n_iter, 1, n_dim) """ -function one_mcmc(log_prob::Function, data::Vector, parameters::Vector, n_iter::Int64, a::Float64) - + + +function one_mcmc(log_prob::Function, + data::Vector, + parameters::Vector, + n_iter::Int64, + a::Float64) + x = ones(( n_iter, length(parameters) )) x[1,:] = [parameters][1] - space = zeros(length(parameters)) + new_parameters = zeros(length(parameters)) for i in 2:n_iter - + + # Propose new parameters for j in 1:length(parameters) - space[j] = x[i-1,j] + a * rand(Uniform(-1, 1), 1)[1] + new_parameters[j] = x[i-1,j] + a * randn(1)[1] end present = log_prob(data, parameters) - future = log_prob(data, space) - - alpha = min( exp(future - present), 1.0) - g = rand() + future = log_prob(data, new_parameters) - if alpha > g + # Accept or reject the proposed parameters + alpha = min( exp(future - present), 1.0) + if rand() < alpha for k in 1:length(parameters) - x[i,k] = space[k] + x[i,k] = new_parameters[k] end - else - x[i,:] = x[i-1,:] + x[i,:] = x[i-1,:] end + end return x[:,:,:] end + """ - run_mcmc(log_prob, data, parameters, n_iter, n_walkers, n_dim, a) +run_mcmc(log_prob, data, seed, n_iter, n_walkers, n_dim, a) -Returns chains with the dimension of (n_iter, n_walkers, n_dim) +Returns chains with dimension of (n_iter, n_walkers, n_dim) """ -function run_mcmc(log_prob::Function, data::Vector, parameters::Vector, n_iter::Int64, n_walkers::Int64, - n_dim::Int64; a::Float64=0.1) + + +function run_mcmc(log_prob::Function, + data::Vector, + seed::Matrix, + n_iter::Int64, + n_walkers::Int64, + n_dim::Int64; + a::Float64=0.1) + + println("Running chains...") mcmc_chains = ones(( n_iter, n_dim, n_walkers )) - seed = rand(n_walkers, n_dim) * 1e-2 .+ transpose(parameters) - @showprogress "Running chains..." for walkers in 1:n_walkers + start_time = time_ns() + + Threads.@threads for walkers in 1:n_walkers one_chain = one_mcmc(log_prob, data, seed[walkers,:], n_iter, a) mcmc_chains[:,:,walkers] = one_chain - sleep(0.1) + end + + elapsed_time = (time_ns() - start_time) / 1_000_000_000 # Convert to seconds + + if elapsed_time < 60 + println("Time elapsed running mcmc: $(elapsed_time) seconds.") + + elseif (elapsed_time >= 60) && (elapsed_time < 3600) + println("Time elapsed running mcmc: $(elapsed_time / 60) minutes.") + + else + println("Time elapsed running mcmc: $(elapsed_time / 3600) hours.") + end chains = permutedims(mcmc_chains, [1, 3, 2]) return chains + end + """ - get_flat_chain(array, burn_in, thin) +get_flat_chain(array, burn_in, thin) Returns the stored chain of MCMC samples. """ + + function get_flat_chain(array::Array; burn_in::Int=1, thin::Int=1) l = ones(( size([i for i in 1:size(array)[1] if i%thin==0])[1], size(array)[2], size(array)[3] )) @@ -101,23 +135,72 @@ end """ - get_gelman_rubin(chains) +gelman_rubin_diagnostic(chains) Get the Gelman Rubin convergence diagnostic of the chains. Returns the Gelman-Rubin number. """ -function get_gelman_rubin(chains) - s_j = var(chains, dims=2) + +function gelman_rubin_diagnostic(chains) + + s_j = var(chains, dims=2) # within chain variance W = mean(s_j) - theta = mean(chains, dims=2) - theta_j = mean(theta, dims=1) + + theta = mean(chains, dims=2) # chain mean + theta_j = mean(theta, dims=1) # grand mean M, N = size(chains)[1], size(chains)[2] - B = N / (M - 1) * sum((theta_j .- theta).^2) + B = N / (M - 1) * sum((theta_j .- theta).^2) # between chain variance + var_theta = (N - 1) / N * W + B / N - R_hat = sqrt(var_theta / W) + R_hat = sqrt(var_theta / W) # Gelma-Rubin statistic + return R_hat + end + + +""" +sampler(pdf, n_samples, intervals, params) + +Get n_samples samples from a 1D Probability Density Distribution +given some parameters. The interval works for the range of samples. + +""" + + + +function sampler(pdf::Function, + n_samples::Number, + interval::Vector, + params::Vector) + + + #states = [] + states = ones(n_samples) + current = rand(Uniform(interval[1], interval[2]), 1)[1] + + for i in 1:n_samples + + #push!(states, current) + states[i] = current + movement = rand(Uniform(interval[1], interval[2]), 1)[1] + + current_prob = pdf(current, params) + movement_prob = pdf(movement, params) + + alpha = min(movement_prob / current_prob, 1.0) + + if rand() < alpha + current = movement + end + end + + return states + +end + + end # module McmcHermes