Skip to content

Commit

Permalink
Add member variables for current proposed and current accepted stats
Browse files Browse the repository at this point in the history
  • Loading branch information
@apulsipher
apulsipher committed Dec 16, 2024
1 parent f94a8c2 commit 1c86fd0
Show file tree
Hide file tree
Showing 3 changed files with 36 additions and 18 deletions.
29 changes: 20 additions & 9 deletions epiworld.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1206,6 +1206,8 @@ class LFMCMC {
std::vector< epiworld_double > m_initial_params; ///< Initial parameters
std::vector< epiworld_double > m_current_proposed_params; ///< Proposed parameters for the next sample
std::vector< epiworld_double > m_current_accepted_params; ///< Most recently accepted parameters (current state of MCMC)
std::vector< epiworld_double > m_current_proposed_stats; ///< Statistics from simulation run with proposed parameters
std::vector< epiworld_double > m_current_accepted_stats; ///< Statistics from simulation run with most recently accepted params

std::vector< epiworld_double > m_observed_stats; ///< Observed statistics

Expand Down Expand Up @@ -1302,6 +1304,8 @@ class LFMCMC {
const std::vector< epiworld_double > & get_initial_params() {return m_initial_params;};
const std::vector< epiworld_double > & get_current_proposed_params() {return m_current_proposed_params;};
const std::vector< epiworld_double > & get_current_accepted_params() {return m_current_accepted_params;};
const std::vector< epiworld_double > & get_current_proposed_stats() {return m_current_proposed_stats;};
const std::vector< epiworld_double > & get_current_accepted_stats() {return m_current_accepted_stats;};

const std::vector< epiworld_double > & get_observed_stats() {return m_observed_stats;};

Expand Down Expand Up @@ -1504,6 +1508,8 @@ class LFMCMC {
std::vector< epiworld_double > m_initial_params; ///< Initial parameters
std::vector< epiworld_double > m_current_proposed_params; ///< Proposed parameters for the next sample
std::vector< epiworld_double > m_current_accepted_params; ///< Most recently accepted parameters (current state of MCMC)
std::vector< epiworld_double > m_current_proposed_stats; ///< Statistics from simulation run with proposed parameters
std::vector< epiworld_double > m_current_accepted_stats; ///< Statistics from simulation run with most recently accepted params

std::vector< epiworld_double > m_observed_stats; ///< Observed statistics

Expand Down Expand Up @@ -1600,6 +1606,8 @@ class LFMCMC {
const std::vector< epiworld_double > & get_initial_params() {return m_initial_params;};
const std::vector< epiworld_double > & get_current_proposed_params() {return m_current_proposed_params;};
const std::vector< epiworld_double > & get_current_accepted_params() {return m_current_accepted_params;};
const std::vector< epiworld_double > & get_current_proposed_stats() {return m_current_proposed_stats;};
const std::vector< epiworld_double > & get_current_accepted_stats() {return m_current_accepted_stats;};

const std::vector< epiworld_double > & get_observed_stats() {return m_observed_stats;};

Expand Down Expand Up @@ -1872,6 +1880,8 @@ inline void LFMCMC<TData>::run(
m_n_stats = m_observed_stats.size();

// Reserving size
m_current_proposed_stats.resize(m_n_stats);
m_current_accepted_stats.resize(m_n_stats);
m_all_sample_drawn_prob.resize(m_n_samples);
m_all_sample_acceptance.resize(m_n_samples, false);
m_all_sample_stats.resize(m_n_samples * m_n_stats);
Expand All @@ -1883,15 +1893,16 @@ inline void LFMCMC<TData>::run(

TData data_i = m_simulation_fun(m_initial_params, this);

std::vector< epiworld_double > proposed_stats_i;
m_summary_fun(proposed_stats_i, data_i, this);
m_summary_fun(m_current_proposed_stats, data_i, this);
m_all_accepted_kernel_scores[0u] = m_kernel_fun(
proposed_stats_i, m_observed_stats, m_epsilon, this
m_current_proposed_stats, m_observed_stats, m_epsilon, this
);

// Recording statistics
for (size_t i = 0u; i < m_n_stats; ++i)
m_all_sample_stats[i] = proposed_stats_i[i];
m_all_sample_stats[i] = m_current_proposed_stats[i];

m_current_accepted_stats = m_current_proposed_stats;

for (size_t k = 0u; k < m_n_params; ++k)
m_all_accepted_params[k] = m_initial_params[k];
Expand All @@ -1916,18 +1927,18 @@ inline void LFMCMC<TData>::run(
m_simulated_data->operator[](i) = data_i;

// Step 3: Generate the summary statistics of the data
m_summary_fun(proposed_stats_i, data_i, this);
m_summary_fun(m_current_proposed_stats, data_i, this);

// Step 4: Compute the hastings ratio using the kernel function
epiworld_double hr = m_kernel_fun(
proposed_stats_i, m_observed_stats, m_epsilon, this
m_current_proposed_stats, m_observed_stats, m_epsilon, this
);

m_all_sample_kernel_scores[i] = hr;

// Storing data
for (size_t k = 0u; k < m_n_stats; ++k)
m_all_sample_stats[i * m_n_stats + k] = proposed_stats_i[k];
m_all_sample_stats[i * m_n_stats + k] = m_current_proposed_stats[k];

// Running Hastings ratio
epiworld_double r = runif();
Expand All @@ -1941,10 +1952,10 @@ inline void LFMCMC<TData>::run(

for (size_t k = 0u; k < m_n_stats; ++k)
m_all_accepted_stats[i * m_n_stats + k] =
proposed_stats_i[k];
m_current_proposed_stats[k];

m_current_accepted_params = m_current_proposed_params;

m_current_accepted_stats = m_current_proposed_stats;
} else
{

Expand Down
4 changes: 4 additions & 0 deletions include/epiworld/math/lfmcmc/lfmcmc-bones.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,6 +144,8 @@ class LFMCMC {
std::vector< epiworld_double > m_initial_params; ///< Initial parameters
std::vector< epiworld_double > m_current_proposed_params; ///< Proposed parameters for the next sample
std::vector< epiworld_double > m_current_accepted_params; ///< Most recently accepted parameters (current state of MCMC)
std::vector< epiworld_double > m_current_proposed_stats; ///< Statistics from simulation run with proposed parameters
std::vector< epiworld_double > m_current_accepted_stats; ///< Statistics from simulation run with most recently accepted params

std::vector< epiworld_double > m_observed_stats; ///< Observed statistics

Expand Down Expand Up @@ -240,6 +242,8 @@ class LFMCMC {
const std::vector< epiworld_double > & get_initial_params() {return m_initial_params;};
const std::vector< epiworld_double > & get_current_proposed_params() {return m_current_proposed_params;};
const std::vector< epiworld_double > & get_current_accepted_params() {return m_current_accepted_params;};
const std::vector< epiworld_double > & get_current_proposed_stats() {return m_current_proposed_stats;};
const std::vector< epiworld_double > & get_current_accepted_stats() {return m_current_accepted_stats;};

const std::vector< epiworld_double > & get_observed_stats() {return m_observed_stats;};

Expand Down
21 changes: 12 additions & 9 deletions include/epiworld/math/lfmcmc/lfmcmc-meat.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,6 +239,8 @@ inline void LFMCMC<TData>::run(
m_n_stats = m_observed_stats.size();

// Reserving size
m_current_proposed_stats.resize(m_n_stats);
m_current_accepted_stats.resize(m_n_stats);
m_all_sample_drawn_prob.resize(m_n_samples);
m_all_sample_acceptance.resize(m_n_samples, false);
m_all_sample_stats.resize(m_n_samples * m_n_stats);
Expand All @@ -250,15 +252,16 @@ inline void LFMCMC<TData>::run(

TData data_i = m_simulation_fun(m_initial_params, this);

std::vector< epiworld_double > proposed_stats_i;
m_summary_fun(proposed_stats_i, data_i, this);
m_summary_fun(m_current_proposed_stats, data_i, this);
m_all_accepted_kernel_scores[0u] = m_kernel_fun(
proposed_stats_i, m_observed_stats, m_epsilon, this
m_current_proposed_stats, m_observed_stats, m_epsilon, this
);

// Recording statistics
for (size_t i = 0u; i < m_n_stats; ++i)
m_all_sample_stats[i] = proposed_stats_i[i];
m_all_sample_stats[i] = m_current_proposed_stats[i];

m_current_accepted_stats = m_current_proposed_stats;

for (size_t k = 0u; k < m_n_params; ++k)
m_all_accepted_params[k] = m_initial_params[k];
Expand All @@ -283,18 +286,18 @@ inline void LFMCMC<TData>::run(
m_simulated_data->operator[](i) = data_i;

// Step 3: Generate the summary statistics of the data
m_summary_fun(proposed_stats_i, data_i, this);
m_summary_fun(m_current_proposed_stats, data_i, this);

// Step 4: Compute the hastings ratio using the kernel function
epiworld_double hr = m_kernel_fun(
proposed_stats_i, m_observed_stats, m_epsilon, this
m_current_proposed_stats, m_observed_stats, m_epsilon, this
);

m_all_sample_kernel_scores[i] = hr;

// Storing data
for (size_t k = 0u; k < m_n_stats; ++k)
m_all_sample_stats[i * m_n_stats + k] = proposed_stats_i[k];
m_all_sample_stats[i * m_n_stats + k] = m_current_proposed_stats[k];

// Running Hastings ratio
epiworld_double r = runif();
Expand All @@ -308,10 +311,10 @@ inline void LFMCMC<TData>::run(

for (size_t k = 0u; k < m_n_stats; ++k)
m_all_accepted_stats[i * m_n_stats + k] =
proposed_stats_i[k];
m_current_proposed_stats[k];

m_current_accepted_params = m_current_proposed_params;

m_current_accepted_stats = m_current_proposed_stats;
} else
{

Expand Down

0 comments on commit 1c86fd0

Please sign in to comment.