Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rename member variables and add new ones #45

Merged
merged 4 commits into from
Dec 16, 2024
Merged

Rename member variables and add new ones #45

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f94a8c2
Rename member variables and associated getters
apulsipher Dec 16, 2024
1c86fd0
Add member variables for current proposed and current accepted stats
apulsipher Dec 16, 2024
a1a79ed
Populate m_all_sample_params
apulsipher Dec 16, 2024
f7f3ab5
Add const keyword to getters
apulsipher Dec 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
190 changes: 104 additions & 86 deletions epiworld.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1203,21 +1203,23 @@ class LFMCMC {

epiworld_double m_epsilon;

std::vector< epiworld_double > m_initial_params; ///< Initial parameters
std::vector< epiworld_double > m_current_params; ///< Parameters for the current sample
std::vector< epiworld_double > m_previous_params; ///< Parameters from the previous sample
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
std::vector< epiworld_double > m_observed_stats; ///< Observed statistics

std::vector< epiworld_double > m_sample_params; ///< Parameter samples
std::vector< epiworld_double > m_sample_stats; ///< Statistic samples
std::vector< bool > m_sample_acceptance; ///< Indicator if sample was accepted
std::vector< epiworld_double > m_sample_drawn_prob; ///< Drawn probabilities (runif()) for each sample
std::vector< epiworld_double > m_sample_kernel_scores; ///< Kernel scores for each sample
std::vector< epiworld_double > m_all_sample_params; ///< Parameter samples
std::vector< epiworld_double > m_all_sample_stats; ///< Statistic samples
std::vector< bool > m_all_sample_acceptance; ///< Indicator if sample was accepted
std::vector< epiworld_double > m_all_sample_drawn_prob; ///< Drawn probabilities (runif()) for each sample
std::vector< epiworld_double > m_all_sample_kernel_scores; ///< Kernel scores for each sample

std::vector< epiworld_double > m_accepted_params; ///< Posterior distribution of parameters from accepted samples
std::vector< epiworld_double > m_accepted_stats; ///< Posterior distribution of statistics from accepted samples
std::vector< epiworld_double > m_accepted_kernel_scores; ///< Kernel scores for each accepted sample
std::vector< epiworld_double > m_all_accepted_params; ///< Posterior distribution of parameters from accepted samples
std::vector< epiworld_double > m_all_accepted_stats; ///< Posterior distribution of statistics from accepted samples
std::vector< epiworld_double > m_all_accepted_kernel_scores; ///< Kernel scores for each accepted sample

// Functions
LFMCMCSimFun<TData> m_simulation_fun;
Expand Down Expand Up @@ -1300,20 +1302,22 @@ class LFMCMC {
epiworld_double get_epsilon() const {return m_epsilon;};

const std::vector< epiworld_double > & get_initial_params() {return m_initial_params;};
const std::vector< epiworld_double > & get_current_params() {return m_current_params;};
const std::vector< epiworld_double > & get_previous_params() {return m_previous_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;};

const std::vector< epiworld_double > & get_sample_params() {return m_sample_params;};
const std::vector< epiworld_double > & get_sample_stats() {return m_sample_stats;};
const std::vector< bool > & get_sample_acceptance() {return m_sample_acceptance;};
const std::vector< epiworld_double > & get_sample_drawn_prob() {return m_sample_drawn_prob;};
const std::vector< epiworld_double > & get_sample_kernel_scores() {return m_sample_kernel_scores;};
const std::vector< epiworld_double > & get_all_sample_params() {return m_all_sample_params;};
const std::vector< epiworld_double > & get_all_sample_stats() {return m_all_sample_stats;};
const std::vector< bool > & get_all_sample_acceptance() {return m_all_sample_acceptance;};
const std::vector< epiworld_double > & get_all_sample_drawn_prob() {return m_all_sample_drawn_prob;};
const std::vector< epiworld_double > & get_all_sample_kernel_scores() {return m_all_sample_kernel_scores;};

const std::vector< epiworld_double > & get_accepted_params() {return m_accepted_params;};
const std::vector< epiworld_double > & get_accepted_stats() {return m_accepted_stats;};
const std::vector< epiworld_double > & get_accepted_kernel_scores() {return m_accepted_kernel_scores;};
const std::vector< epiworld_double > & get_all_accepted_params() {return m_all_accepted_params;};
const std::vector< epiworld_double > & get_all_accepted_stats() {return m_all_accepted_stats;};
const std::vector< epiworld_double > & get_all_accepted_kernel_scores() {return m_all_accepted_kernel_scores;};

std::vector< TData > * get_simulated_data() {return m_simulated_data;};

Expand Down Expand Up @@ -1501,21 +1505,23 @@ class LFMCMC {

epiworld_double m_epsilon;

std::vector< epiworld_double > m_initial_params; ///< Initial parameters
std::vector< epiworld_double > m_current_params; ///< Parameters for the current sample
std::vector< epiworld_double > m_previous_params; ///< Parameters from the previous sample
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
std::vector< epiworld_double > m_observed_stats; ///< Observed statistics

std::vector< epiworld_double > m_sample_params; ///< Parameter samples
std::vector< epiworld_double > m_sample_stats; ///< Statistic samples
std::vector< bool > m_sample_acceptance; ///< Indicator if sample was accepted
std::vector< epiworld_double > m_sample_drawn_prob; ///< Drawn probabilities (runif()) for each sample
std::vector< epiworld_double > m_sample_kernel_scores; ///< Kernel scores for each sample
std::vector< epiworld_double > m_all_sample_params; ///< Parameter samples
std::vector< epiworld_double > m_all_sample_stats; ///< Statistic samples
std::vector< bool > m_all_sample_acceptance; ///< Indicator if sample was accepted
std::vector< epiworld_double > m_all_sample_drawn_prob; ///< Drawn probabilities (runif()) for each sample
std::vector< epiworld_double > m_all_sample_kernel_scores; ///< Kernel scores for each sample

std::vector< epiworld_double > m_accepted_params; ///< Posterior distribution of parameters from accepted samples
std::vector< epiworld_double > m_accepted_stats; ///< Posterior distribution of statistics from accepted samples
std::vector< epiworld_double > m_accepted_kernel_scores; ///< Kernel scores for each accepted sample
std::vector< epiworld_double > m_all_accepted_params; ///< Posterior distribution of parameters from accepted samples
std::vector< epiworld_double > m_all_accepted_stats; ///< Posterior distribution of statistics from accepted samples
std::vector< epiworld_double > m_all_accepted_kernel_scores; ///< Kernel scores for each accepted sample

// Functions
LFMCMCSimFun<TData> m_simulation_fun;
Expand Down Expand Up @@ -1598,20 +1604,22 @@ class LFMCMC {
epiworld_double get_epsilon() const {return m_epsilon;};

const std::vector< epiworld_double > & get_initial_params() {return m_initial_params;};
const std::vector< epiworld_double > & get_current_params() {return m_current_params;};
const std::vector< epiworld_double > & get_previous_params() {return m_previous_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;};

const std::vector< epiworld_double > & get_sample_params() {return m_sample_params;};
const std::vector< epiworld_double > & get_sample_stats() {return m_sample_stats;};
const std::vector< bool > & get_sample_acceptance() {return m_sample_acceptance;};
const std::vector< epiworld_double > & get_sample_drawn_prob() {return m_sample_drawn_prob;};
const std::vector< epiworld_double > & get_sample_kernel_scores() {return m_sample_kernel_scores;};
const std::vector< epiworld_double > & get_all_sample_params() {return m_all_sample_params;};
const std::vector< epiworld_double > & get_all_sample_stats() {return m_all_sample_stats;};
const std::vector< bool > & get_all_sample_acceptance() {return m_all_sample_acceptance;};
const std::vector< epiworld_double > & get_all_sample_drawn_prob() {return m_all_sample_drawn_prob;};
const std::vector< epiworld_double > & get_all_sample_kernel_scores() {return m_all_sample_kernel_scores;};

const std::vector< epiworld_double > & get_accepted_params() {return m_accepted_params;};
const std::vector< epiworld_double > & get_accepted_stats() {return m_accepted_stats;};
const std::vector< epiworld_double > & get_accepted_kernel_scores() {return m_accepted_kernel_scores;};
const std::vector< epiworld_double > & get_all_accepted_params() {return m_all_accepted_params;};
const std::vector< epiworld_double > & get_all_accepted_stats() {return m_all_accepted_stats;};
const std::vector< epiworld_double > & get_all_accepted_kernel_scores() {return m_all_accepted_kernel_scores;};

std::vector< TData > * get_simulated_data() {return m_simulated_data;};

Expand Down Expand Up @@ -1858,43 +1866,50 @@ inline void LFMCMC<TData>::run(
if (seed >= 0)
this->seed(seed);

m_current_params.resize(m_n_params);
m_previous_params.resize(m_n_params);
m_current_proposed_params.resize(m_n_params);
m_current_accepted_params.resize(m_n_params);

if (m_simulated_data != nullptr)
m_simulated_data->resize(m_n_samples);

m_previous_params = m_initial_params;
m_current_params = m_initial_params;
m_current_accepted_params = m_initial_params;
m_current_proposed_params = m_initial_params;

// Computing the baseline sufficient statistics
m_summary_fun(m_observed_stats, m_observed_data, this);
m_n_stats = m_observed_stats.size();

// Reserving size
m_sample_drawn_prob.resize(m_n_samples);
m_sample_acceptance.resize(m_n_samples, false);
m_sample_stats.resize(m_n_samples * m_n_stats);
m_sample_kernel_scores.resize(m_n_samples);

m_accepted_params.resize(m_n_samples * m_n_params);
m_accepted_stats.resize(m_n_samples * m_n_stats);
m_accepted_kernel_scores.resize(m_n_samples);
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_params.resize(m_n_samples * m_n_params);
m_all_sample_stats.resize(m_n_samples * m_n_stats);
m_all_sample_kernel_scores.resize(m_n_samples);

m_all_accepted_params.resize(m_n_samples * m_n_params);
m_all_accepted_stats.resize(m_n_samples * m_n_stats);
m_all_accepted_kernel_scores.resize(m_n_samples);

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_accepted_kernel_scores[0u] = m_kernel_fun(
proposed_stats_i, m_observed_stats, m_epsilon, this
m_summary_fun(m_current_proposed_stats, data_i, this);
m_all_accepted_kernel_scores[0u] = m_kernel_fun(
m_current_proposed_stats, m_observed_stats, m_epsilon, this
);

// Recording statistics
for (size_t i = 0u; i < m_n_stats; ++i)
m_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_accepted_params[k] = m_initial_params[k];
m_all_accepted_params[k] = m_initial_params[k];

for (size_t k = 0u; k < m_n_params; ++k)
m_all_sample_params[k] = m_initial_params[k];

// Init progress bar
progress_bar = Progress(m_n_samples, 80);
Expand All @@ -1905,59 +1920,62 @@ inline void LFMCMC<TData>::run(
// Run LFMCMC
for (size_t i = 1u; i < m_n_samples; ++i)
{
// Step 1: Generate a proposal and store it in m_current_params
m_proposal_fun(m_current_params, m_previous_params, this);
// Step 1: Generate a proposal and store it in m_current_proposed_params
m_proposal_fun(m_current_proposed_params, m_current_accepted_params, this);

// Step 2: Using m_current_params, simulate data
TData data_i = m_simulation_fun(m_current_params, this);
// Step 2: Using m_current_proposed_params, simulate data
TData data_i = m_simulation_fun(m_current_proposed_params, this);

// Are we storing the data?
if (m_simulated_data != nullptr)
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_sample_kernel_scores[i] = hr;
m_all_sample_kernel_scores[i] = hr;

// Storing data
for (size_t k = 0u; k < m_n_params; ++k)
m_all_sample_params[i * m_n_params + k] = m_current_proposed_params[k];

for (size_t k = 0u; k < m_n_stats; ++k)
m_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();
m_sample_drawn_prob[i] = r;
m_all_sample_drawn_prob[i] = r;

// Step 5: Update if likely
if (r < std::min(static_cast<epiworld_double>(1.0), hr / m_accepted_kernel_scores[i - 1u]))
if (r < std::min(static_cast<epiworld_double>(1.0), hr / m_all_accepted_kernel_scores[i - 1u]))
{
m_accepted_kernel_scores[i] = hr;
m_sample_acceptance[i] = true;
m_all_accepted_kernel_scores[i] = hr;
m_all_sample_acceptance[i] = true;

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

m_previous_params = m_current_params;
m_all_accepted_stats[i * m_n_stats + k] =
m_current_proposed_stats[k];

m_current_accepted_params = m_current_proposed_params;
m_current_accepted_stats = m_current_proposed_stats;
} else
{

for (size_t k = 0u; k < m_n_stats; ++k)
m_accepted_stats[i * m_n_stats + k] =
m_accepted_stats[(i - 1) * m_n_stats + k];
m_all_accepted_stats[i * m_n_stats + k] =
m_all_accepted_stats[(i - 1) * m_n_stats + k];

m_accepted_kernel_scores[i] = m_accepted_kernel_scores[i - 1u];
m_all_accepted_kernel_scores[i] = m_all_accepted_kernel_scores[i - 1u];
}


for (size_t k = 0u; k < m_n_params; ++k)
m_accepted_params[i * m_n_params + k] = m_previous_params[k];
m_all_accepted_params[i * m_n_params + k] = m_current_accepted_params[k];

if (verbose) {
progress_bar.next();
Expand Down Expand Up @@ -2167,7 +2185,7 @@ inline void LFMCMC<TData>::print(size_t burnin) const
std::vector< epiworld_double > par_i(n_samples_print);
for (size_t i = burnin; i < m_n_samples; ++i)
{
par_i[i-burnin] = m_accepted_params[i * m_n_params + k];
par_i[i-burnin] = m_all_accepted_params[i * m_n_params + k];
summ_params[k * 3] += par_i[i-burnin]/n_samples_dbl;
}

Expand All @@ -2189,7 +2207,7 @@ inline void LFMCMC<TData>::print(size_t burnin) const
std::vector< epiworld_double > stat_k(n_samples_print);
for (size_t i = burnin; i < m_n_samples; ++i)
{
stat_k[i-burnin] = m_accepted_stats[i * m_n_stats + k];
stat_k[i-burnin] = m_all_accepted_stats[i * m_n_stats + k];
summ_stats[k * 3] += stat_k[i-burnin]/n_samples_dbl;
}

Expand Down Expand Up @@ -2412,7 +2430,7 @@ inline std::vector< epiworld_double > LFMCMC<TData>::get_mean_params()
for (size_t k = 0u; k < m_n_params; ++k)
{
for (size_t i = 0u; i < m_n_samples; ++i)
res[k] += (this->m_accepted_params[k + m_n_params * i])/
res[k] += (this->m_all_accepted_params[k + m_n_params * i])/
static_cast< epiworld_double >(m_n_samples);
}

Expand All @@ -2428,7 +2446,7 @@ inline std::vector< epiworld_double > LFMCMC<TData>::get_mean_stats()
for (size_t k = 0u; k < m_n_stats; ++k)
{
for (size_t i = 0u; i < m_n_samples; ++i)
res[k] += (this->m_accepted_stats[k + m_n_stats * i])/
res[k] += (this->m_all_accepted_stats[k + m_n_stats * i])/
static_cast< epiworld_double >(m_n_samples);
}

Expand Down
Loading
Loading