diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..2f97691 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,14 @@ +name: CI +on: [push, pull_request] +permissions: + contents: read +jobs: + regression: + strategy: + matrix: + os: [ubuntu-latest, macos-latest] + runs-on: ${{ matrix.os }} + steps: + - uses: actions/checkout@v4 + - run: c++ -std=c++17 -O2 -Wall -Wextra -Wpedantic -pthread -Ibest tests/regression.cpp -o regression && ./regression + - run: c++ -std=c++17 -O1 -g -fsanitize=address,undefined -fno-omit-frame-pointer -pthread -Ibest tests/regression.cpp -o sanitized && ./sanitized diff --git a/.gitignore b/.gitignore index 005815c..f55635e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,67 +1,4 @@ -# Compiled Object files -*.slo -*.lo -*.o -*.obj - -# Precompiled Headers -*.gch -*.pch - -# Compiled Dynamic libraries -*.so -*.dylib -*.dll - -# Fortran module files -*.mod -*.smod - -# Compiled Static libraries -*.lai -*.la -*.a -*.lib - -# Executables -*.exe -*.out -*.app - -# OSX / XCode +build*/ +install/ +bazel-* .DS_Store -*.pbxuser -xcuserdata/ - -# Visual Studio user-specific files -*.suo -*.user -*.userosscache -*.sln.docstates -*.vcxproj.filters - -# Visual Studio cache/options directory -.vs/ - -# Visual C++ cache files -ipch/ -*.aps -*.ncb -*.opendb -*.opensdf -*.sdf -*.cachefile -*.VC.db -*.VC.VC.opendb - -# Build results -[Dd]ebug/ -[Dd]ebugPublic/ -[Rr]elease/ -[Rr]eleases/ -x64/ -x86/ -bld/ -[Bb]in/ -[Oo]bj/ -[Ll]og/ diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..fe71088 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,13 @@ +# Changelog + +## Unreleased + +- Replace speculative parallel proposals with sequential adaptive Metropolis-within-Gibbs. Remove shared random-generator races and proposal-selection bias. Remove the thread-count argument from `AMWG::Init`; callers now pass only the starting state and log-posterior callback. +- `Sample(n)` records exactly n completed sweeps; `Burn(n)` advances without retaining samples. Adaptation tracks sweeps independently of chain storage. Reinitialization resets the seeded generator and adaptation. Sampling no longer prints progress. +- Validate initialization, observations, batch sizes and statistical parameters. Callback exceptions propagate on the calling thread; a partially completed sweep may have advanced the state, but is not added to the chain. +- Enforce both sigma bounds and nu >= 1, evaluate the model in log space, retain double precision, and reuse Student-t normalization across observations. +- Own model observations to prevent dangling references. BEST cannot be copied or moved because its callback refers to the model. +- Define the empirical 95% interval as the shortest window containing ceil(0.95*n) observations, choosing the first window on ties. Reject empty/nonfinite inputs. +- Compute population standard deviation without a temporary allocation. Add numerical, lifecycle, seeded-distribution and independent-chain regression checks. + +These corrections change seeded sequences, inference results and invalid-input behavior. Old results should be recomputed. Reproducibility is within the same implementation and standard library; C++ random-distribution algorithms are not portable bit-for-bit. diff --git a/best/amwg.hpp b/best/amwg.hpp index c572882..0e35372 100644 --- a/best/amwg.hpp +++ b/best/amwg.hpp @@ -1,100 +1,79 @@ #pragma once - +#include #include -#include -#include +#include +#include +#include #include +#include #include -#include +#include #include +#include #include static const size_t kDefaultBatchSize = 50; -/** - * Implementation of an Adaptive Metropolis within Gibbs sampler. - */ +// Sequential adaptive Metropolis-within-Gibbs; one proposal per coordinate. +// Instances are not safe for concurrent mutation. Independent chains may run +// on separate threads with independent seeds and posterior state. template class AMWG { + static_assert(std::is_floating_point::value, "floating point required"); + static_assert(NumParams > 0, "at least one parameter required"); public: - using ParamArray = std::array; using PosteriorFunc = std::function; - - AMWG(size_t batchSize = kDefaultBatchSize, uint32_t seed = std::mt19937::default_seed) - : rng_(seed) - , rand_(0.0, 1.0) - , currentPosteriorDensity_(0.0) - , batchSize_(batchSize) - , batchCount_(0) - , logSD_() - , acceptanceCount_() - , running_(false) - , epoch_(0) - , completedCount_(0) - , curParam_(0) - {}; - - ~AMWG() { - { - std::lock_guard lock(mutex_); - running_ = false; - epoch_ = SIZE_MAX; - } - - cv_.notify_all(); - - for (auto&& t : threads_) { - if (t.joinable()) t.join(); - } + explicit AMWG(size_t batchSize = kDefaultBatchSize, + uint32_t seed = std::mt19937::default_seed) + : rng_(seed), seed_(seed), batchSize_(batchSize) { + if (!batchSize) throw std::invalid_argument("batch size must be positive"); } - - void Init(std::array startValues, PosteriorFunc posteriorFunc, uint32_t threads = std::thread::hardware_concurrency()) { - state_ = startValues; - posteriorFunc_ = posteriorFunc; + void Init(ParamArray start, PosteriorFunc posterior) { + if (!posterior) throw std::invalid_argument("posterior is required"); + for (auto x : start) + if (!std::isfinite(x)) throw std::invalid_argument("nonfinite initial state"); + RealType density = posterior(start); + if (!std::isfinite(density)) throw std::invalid_argument("initial log density must be finite"); + posterior_ = std::move(posterior); + state_ = start; + density_ = density; chain_.clear(); - currentPosteriorDensity_ = posteriorFunc_(state_); - - running_ = true; - - threads_.reserve(threads); - threadStates_.resize(threads); - - for (uint32_t threadId = 0; threadId < threads; threadId++) { - threads_.emplace_back(std::thread(&AMWG::threadWorker, this, threadId)); - } - }; - - size_t NextSample(); - void Sample(size_t n); - void Burn(size_t n); - RealType posterior_density() { return currentPosteriorDensity_; }; - std::vector& chain() { return chain_; }; - + rng_.seed(seed_); + logSD_.fill(0); + accepted_.fill(0); + batchCount_ = withinBatch_ = 0; + initialized_ = true; + } + size_t NextSample() { + requireInit(); + chain_.reserve(chain_.size() + 1); + step(); + chain_.push_back(state_); + return NumParams; + } + void Sample(size_t n) { + requireInit(); + if (n > chain_.max_size() - chain_.size()) throw std::length_error("chain too large"); + chain_.reserve(chain_.size() + n); + for (size_t i=0; i& chain() { return chain_; } + const std::vector& chain() const { return chain_; } private: - - void threadWorker(uint32_t threadId); - void createProposal(uint32_t threadId); - - PosteriorFunc posteriorFunc_; + void requireInit() const { if (!initialized_) throw std::logic_error("Init must be called first"); } + void step(); std::mt19937 rng_; - std::uniform_real_distribution rand_; + uint32_t seed_; + size_t batchSize_, batchCount_=0, withinBatch_=0; + bool initialized_=false; + PosteriorFunc posterior_; + ParamArray state_{}, logSD_{}; + std::array accepted_{}; + RealType density_=0; std::vector chain_; - RealType currentPosteriorDensity_; - size_t batchSize_; - size_t batchCount_; - ParamArray state_; - ParamArray logSD_; - std::array acceptanceCount_; - - std::atomic_bool running_; - std::atomic epoch_; - std::atomic completedCount_; - std::atomic curParam_; - std::vector threads_; - std::vector> threadStates_; - std::mutex mutex_; - std::condition_variable cv_; }; - #include "amwg.inl" diff --git a/best/amwg.inl b/best/amwg.inl index b195132..b5c61cd 100644 --- a/best/amwg.inl +++ b/best/amwg.inl @@ -1,127 +1,31 @@ -#include -#include -#include - +#pragma once template -void AMWG::threadWorker(uint32_t threadId) { - int curEpoch = 0; - - while (running_) { - { - std::unique_lock lock(mutex_); - while (curEpoch >= epoch_) cv_.wait(lock); +void AMWG::step() { + std::uniform_real_distribution uniform(0, 1); + for (size_t i=0; i normal(state_[i], std::exp(logSD_[i])); + proposal[i] = normal(rng_); + if (!std::isfinite(proposal[i])) continue; + RealType next = posterior_(proposal); + // Reject invalid/out-of-support proposals, including +infinity and NaN. + RealType u = uniform(rng_); + if (std::isfinite(next) && (next >= density_ || std::log(u) < next-density_)) { + state_ = proposal; + density_ = next; + ++accepted_[i]; } - - if (!running_) return; - - createProposal(threadId); - - ++curEpoch; - ++completedCount_; } -} - -template -void AMWG::createProposal(uint32_t threadId) { - ParamArray localState = state_; - - // Modify one parameter in the current state by sampling from a Gaussian distribution - // with mean and standard deviation of the current parameter - RealType stdev = exp(logSD_[curParam_]); - std::normal_distribution normal(state_[curParam_], stdev); - localState[curParam_] = normal(rng_); - - // Measure the posterior density of this proposal - RealType proposedPosteriorDensity = posteriorFunc_(localState); - bool accept; - if (!isfinite(proposedPosteriorDensity)) { - proposedPosteriorDensity = RealType(0.0); - localState[curParam_] = state_[curParam_]; - accept = false; - } else { - // If this proposal is better than the previous, we always accept. - // Otherwise, we accept proportional to the likelihood ratio - RealType acceptProb = exp(proposedPosteriorDensity - currentPosteriorDensity_); - RealType rand = rand_(rng_); - accept = (acceptProb > rand); - } - - threadStates_[threadId] = std::make_tuple(localState[curParam_], proposedPosteriorDensity, accept); -} - -template -size_t AMWG::NextSample() { - static const RealType kTargetAcceptRate(0.44f); - - size_t threadCount = threads_.size(); - - chain_.emplace_back(state_); - - size_t steps = 0; - - for (curParam_ = 0; curParam_ < NumParams; curParam_++) { - { - std::lock_guard lock(mutex_); - //std::unique_lock lock(mutex_); - ++epoch_; - completedCount_ = 0; - } - - cv_.notify_all(); - - while (completedCount_ < threadCount) std::this_thread::yield(); - - for (auto&& result : threadStates_) { - RealType proposal = std::get<0>(result); - RealType proposedPosteriorDensity = std::get<1>(result); - bool accepted = std::get<2>(result); - - ++steps; - - // Iterate over each result, stopping on the first acceptance - if (accepted) { - acceptanceCount_[curParam_]++; - currentPosteriorDensity_ = proposedPosteriorDensity; - state_[curParam_] = proposal; - break; - } + if (++withinBatch_ == batchSize_) { + ++batchCount_; + RealType delta = std::min(RealType(0.01), RealType(1)/std::sqrt(RealType(batchCount_))); + for (size_t i=0; i= RealType(0.44) ? delta : -delta; + // Keep proposal scales representable over very long chains. + logSD_[i] = std::max(std::log(std::numeric_limits::min())/2, + std::min(std::log(std::numeric_limits::max())/2, logSD_[i])); + accepted_[i] = 0; } + withinBatch_ = 0; } - - // If this iteration completes a batch, update logSD and reset acceptanceCount - if (chain_.size() % batchSize_ == 0) { - batchCount_++; - RealType ooBatchCount = std::min(RealType(0.01), RealType(1.0) / std::sqrt(static_cast(batchCount_))); - - for (size_t i = 0; i < NumParams; i++) { - RealType pctOfBatch = static_cast(acceptanceCount_[i]) / static_cast(batchSize_); - logSD_[i] += (pctOfBatch >= kTargetAcceptRate) ? ooBatchCount : -ooBatchCount; - acceptanceCount_[i] = 0; - } - } - - return steps; -} - -template -void AMWG::Sample(size_t n) { - chain_.reserve(chain_.size() + n); - - size_t totalSteps = n * NumParams; - - - for (size_t i = 0, s = 0; i < n && s < totalSteps; i++) { - s += NextSample(); - if (i % 10 == 0) { - double pct = (static_cast(s) / static_cast(totalSteps)) * 100.0; - std::cout << std::fixed << std::setprecision(2) << pct << "%" << std::endl; - } - } -} - -template -void AMWG::Burn(size_t n) { - std::vector prevChain(std::move(chain_)); - Sample(n); - chain_ = std::move(prevChain); } diff --git a/best/best.hpp b/best/best.hpp index e8e929e..e31a684 100644 --- a/best/best.hpp +++ b/best/best.hpp @@ -1,114 +1,65 @@ #pragma once - -#include -#include -#include - #include "amwg.hpp" #include "stats.hpp" -/** - * Implementation of "Bayesian estimation supersedes the t test". Estimate the - * difference between two sets of real values. - */ -template -class BEST { - // Length of startValues array below - static const size_t kParamCount = 5; - +// Two-group Student-t BEST model. Owns observations; copying/moving is disabled +// because the sampler's posterior callback refers to this model. +template class BEST { + static const size_t kParamCount=5; public: - - using RealType = typename Container::value_type; - using ParamArray = typename AMWG::ParamArray; - - BEST(const Container& y1, const Container& y2, size_t batchSize = kDefaultBatchSize, uint32_t seed = std::mt19937::default_seed) - : sampler_(batchSize, seed) - , y1_(y1) - , y2_(y2) - { - size_t jointSize = y1.size() + y2.size(); - - // Find the mean of {y1,y2} - RealType sum = std::accumulate(y1.begin(), y1.end(), RealType(0.0)) + - std::accumulate(y2.begin(), y2.end(), RealType(0.0)); - RealType mean = RealType(sum / static_cast(jointSize)); - - // Find the standard deviation of {y1,y2} - std::vector diff(jointSize); - std::transform(y1.begin(), y1.end(), diff.begin(), [mean](RealType x) { return x - mean; }); - std::transform(y2.begin(), y2.end(), diff.begin() + y1.size(), [mean](RealType x) { return x - mean; }); - RealType sqSum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0f); - RealType stdev = std::sqrt(sqSum / jointSize); - - meanMu_ = mean; - RealType sdMu = stdev; - scaledSdMu_ = sdMu * RealType(1000000.0); - sigmaLow_ = sdMu / RealType(1000.0); - sigmaHigh_ = sdMu * RealType(1000.0); - - // Initial parameter values: mu_y1, mu_y2, sigma_y1, sigma_y2, dof - std::array startValues = { - stats::mean(y1), - stats::mean(y2), - stats::stdev(y1), - stats::stdev(y2), - RealType(5.0) - }; - - sampler_.Init(startValues, std::bind(&BEST::JointPosterior, this, std::placeholders::_1)); + using RealType=typename Container::value_type; + using ParamArray=typename AMWG::ParamArray; + BEST(const Container& y1,const Container& y2,size_t batchSize=kDefaultBatchSize, + uint32_t seed=std::mt19937::default_seed) + : y1_(y1), y2_(y2), sampler_(batchSize,seed) { + stats::detail::values(y1_); stats::detail::values(y2_); + std::vector joint(y1_.begin(),y1_.end()); + joint.insert(joint.end(),y2_.begin(),y2_.end()); + meanMu_=stats::mean(joint); + RealType sd=stats::stdev(joint); + if (!(sd>0) || !std::isfinite(sd)) throw std::invalid_argument("pooled variance must be positive and finite"); + scaledSdMu_=sd*RealType(1000000); + sigmaLow_=sd/RealType(1000); sigmaHigh_=sd*RealType(1000); + if (!(sigmaLow_>0) || !std::isfinite(scaledSdMu_) || !std::isfinite(sigmaHigh_)) + throw std::invalid_argument("data scale cannot represent prior bounds"); + ParamArray start={{stats::mean(y1_),stats::mean(y2_), + std::max(sigmaLow_,stats::stdev(y1_)),std::max(sigmaLow_,stats::stdev(y2_)),RealType(5)}}; + sampler_.Init(start,[this](const ParamArray& p){return LogPosterior(p);}); + } + BEST(const BEST&)=delete; + BEST& operator=(const BEST&)=delete; + BEST(BEST&&)=delete; + BEST& operator=(BEST&&)=delete; + void Burn(size_t n) { sampler_.Burn(n); } + void Sample(size_t n) { sampler_.Sample(n); } + const std::vector& chain() const { return sampler_.chain(); } + void ComputeStats(std::pair& hdi,RealType& mean) const { + std::vector diff; + diff.reserve(chain().size()); + for (const auto& p:chain()) diff.push_back(p[0]-p[1]); + hdi=stats::highestDensityInterval(diff); mean=stats::mean(diff); } - - void Burn(size_t n) { sampler_.Burn(n); }; - - void Sample(size_t n) { sampler_.Sample(n); }; - - void ComputeStats(std::pair& hdi, RealType& mean) { - std::vector& chain = sampler_.chain(); - - std::vector muDiff(chain.size()); - std::transform(chain.begin(), chain.end(), muDiff.begin(), [](ParamArray params) { return params[0] - params[1]; }); - hdi = stats::highestDensityInterval(muDiff); - mean = stats::mean(muDiff); + RealType LogPosterior(const ParamArray& p) const { + for(auto x:p) if(!std::isfinite(x)) return -std::numeric_limits::infinity(); + if(p[2]sigmaHigh_||p[3]sigmaHigh_||p[4]<1) + return -std::numeric_limits::infinity(); + return -std::log(RealType(29))-(p[4]-1)/RealType(29) + +posterior(p[2],p[0],p[4],y1_)+posterior(p[3],p[1],p[4],y2_); } - private: - - RealType Posterior(RealType sigma, RealType mu, RealType nu, const Container& data) { - RealType logP = log(stats::UniformPDF(sigmaLow_, sigmaHigh_)); - logP += log(stats::NormalPDF(mu, meanMu_, scaledSdMu_)); - - RealType ooSD = RealType(1.0) / sigma; - for (RealType val : data) { - logP += log(ooSD * stats::StudentTPDF((val - mu) * ooSD, nu)); + RealType posterior(RealType sigma,RealType mu,RealType nu,const Container& data) const { + RealType logP=-std::log(sigmaHigh_-sigmaLow_)+stats::NormalLogPDF(mu,meanMu_,scaledSdMu_); + // The Student-t normalization is constant across all observations. + RealType norm=stats::StudentTLogPDF(RealType(0),nu)-std::log(sigma); + for(auto x:data) { + RealType z=(x-mu)/sigma; + RealType q=2*std::log(std::abs(z))-std::log(nu); + RealType tail=q>0 ? q+std::log1p(std::exp(-q)) : std::log1p(std::exp(q)); + logP+=norm-(nu+1)/2*tail; } - return logP; - }; - - RealType JointPosterior(const ParamArray& params) { - // A trick to get an exponentially distributed prior on nu that starts at 1 - const RealType kOOTwentyNine(1.0f / 29.0f); - - RealType mu1 = params[0]; - RealType mu2 = params[1]; - RealType sigma1 = params[2]; - RealType sigma2 = params[3]; - RealType nu = params[4]; - - if (sigma1 < sigmaLow_ || sigma2 < sigmaLow_) return -std::numeric_limits::infinity(); - - RealType logP = log(stats::ExponentialPDF(nu - RealType(1.0), kOOTwentyNine)); - logP += Posterior(sigma1, mu1, nu, y1_); - logP += Posterior(sigma2, mu2, nu, y2_); - - return logP; - }; - - const Container& y1_; - const Container& y2_; - AMWG sampler_; - RealType meanMu_; - RealType scaledSdMu_; - RealType sigmaLow_; - RealType sigmaHigh_; + } + Container y1_,y2_; + AMWG sampler_; + RealType meanMu_,scaledSdMu_,sigmaLow_,sigmaHigh_; }; diff --git a/best/main.cpp b/best/main.cpp index 89529e9..35c29a4 100644 --- a/best/main.cpp +++ b/best/main.cpp @@ -1,54 +1,30 @@ +#include "best.hpp" #include -#include #include -#include +#include +#include #include -#include "amwg.hpp" -#include "best.hpp" - -std::vector readFile(const char* filename); - -int main(int argc, const char* argv[]) { - /*static auto y1 = std::vector{1.96f, 2.06f, 2.03f, 2.11f, 1.88f, 1.88f, 2.08f, 1.93f, 2.03f, 2.03f, 2.03f, 2.08f, 2.03f, 2.11f, 1.93f}; - static auto y2 = std::vector{1.83f, 1.93f, 1.88f, 1.85f, 1.85f, 1.91f, 1.91f, 1.85f, 1.78f, 1.91f, 1.93f, 1.80f, 1.80f, 1.85f, 1.93f, - 1.85f, 1.83f, 1.85f, 1.91f, 1.85f, 1.91f, 1.85f, 1.80f, 1.80f, 1.85f};*/ - - if (argc != 3) { - std::cout << "Usage: best " << std::endl; - return 1; - } - - auto y1 = readFile(argv[1]); - auto y2 = readFile(argv[2]); - - uint32_t now = static_cast(std::chrono::system_clock::now().time_since_epoch().count()); - - BEST> best(y1, y2, kDefaultBatchSize, now); - - std::cout << "Running burn-in" << std::endl; - best.Burn(5000); - std::cout << "Running sampler" << std::endl; - best.Sample(5000); - - std::pair hdi; - float mean; - best.ComputeStats(hdi, mean); - - std::cout.precision(std::numeric_limits::max_digits10); - std::cout << "hdi = " << hdi.first << "," << hdi.second << ", mean = " << mean << std::endl; - return 0; -} - -std::vector readFile(const char* filename) { - std::ifstream infile(filename); - - std::vector values; - float value; - - while (infile >> value) { - values.emplace_back(value); - } - +std::vector readFile(const char* filename) { + std::ifstream input(filename); + if(!input) throw std::runtime_error(std::string("cannot open ")+filename); + std::vector values; + double value; + while(input>>value) values.push_back(value); + if(!input.eof()) throw std::runtime_error(std::string("invalid numeric data in ")+filename); + if(values.empty()) throw std::runtime_error(std::string("empty input: ")+filename); return values; } +int main(int argc,const char* argv[]) { + if(argc!=3) {std::cerr<<"Usage: best \n";return 1;} + try { + auto y1=readFile(argv[1]), y2=readFile(argv[2]); + uint32_t seed=static_cast(std::chrono::system_clock::now().time_since_epoch().count()); + BEST> model(y1,y2,kDefaultBatchSize,seed); + model.Burn(5000);model.Sample(5000); + std::pair hdi;double mean; + model.ComputeStats(hdi,mean); + std::cout.precision(std::numeric_limits::max_digits10); + std::cout<<"hdi = "< - -#define M_LOG_2PI 1.837877066409345483560659472811235279 - -template -RealType stats::NormalPDF(RealType x, RealType mean, RealType std) { - return exp(RealType(-0.5) * RealType(M_LOG_2PI) - log(std) - pow(x - mean, 2) / (RealType(2.0) * std * std)); +#pragma once +namespace stats { +namespace detail { +template void positive(T x) { + if (!(x > 0) || !std::isfinite(x)) throw std::invalid_argument("finite positive parameter required"); } - -template -RealType stats::ExponentialPDF(RealType x, RealType rate) { - return std::max(RealType(0.0), rate * exp(-rate * x)); +template void values(const C& c) { + static_assert(std::is_floating_point::value, "floating point required"); + if (c.empty()) throw std::invalid_argument("empty data"); + for (auto x:c) if (!std::isfinite(x)) throw std::invalid_argument("nonfinite data"); } - -template -RealType stats::UniformPDF(RealType a, RealType b) { - return RealType(1.0) / (b - a); } - -template -RealType stats::StudentTPDF(RealType x, RealType dof) { - const RealType one = RealType(1.0); - const RealType half = RealType(0.5); - - return one / (std::sqrt(dof) * Beta(half, dof * half)) * pow(one + x * x / dof, -((dof + one) * half)); +template T NormalLogPDF(T x,T m,T sd) { + detail::positive(sd); + T z=(x-m)/sd; + return -T(0.91893853320467274178L)-std::log(sd)-T(0.5)*z*z; } - -template -RealType stats::Beta(RealType x, RealType y) { - return exp(std::lgamma(x) + std::lgamma(y) - std::lgamma(x + y)); +template T NormalPDF(T x,T m,T sd) { return std::exp(NormalLogPDF(x,m,sd)); } +template T ExponentialPDF(T x,T rate) { + detail::positive(rate); + return x<0 ? T(0) : rate*std::exp(-rate*x); } - -template -typename Container::value_type stats::mean(const Container& c) { - using RealType = typename Container::value_type; - - RealType sum = std::accumulate(c.begin(), c.end(), RealType(0.0)); - return sum / static_cast(c.size()); +template T UniformPDF(T a,T b) { + if (!std::isfinite(a)||!std::isfinite(b)||!(a -typename Container::value_type stats::stdev(const Container& c) { - using RealType = typename Container::value_type; - - RealType m = mean(c); - - std::vector diff(c.size()); - std::transform(c.begin(), c.end(), diff.begin(), [m](RealType x) { return x - m; }); - - RealType sqSum = std::inner_product(diff.begin(), diff.end(), diff.begin(), RealType(0.0)); - return std::sqrt(sqSum / static_cast(c.size())); +template T StudentTLogPDF(T x,T dof) { + detail::positive(dof); + // Evaluate the quadratic in log space to avoid overflow in x*x. + T q = T(2)*std::log(std::abs(x))-std::log(dof); + T tail = q>0 ? q+std::log1p(std::exp(-q)) : std::log1p(std::exp(q)); + // Asymptotic gamma ratio avoids catastrophic cancellation at large dof. + T normalizer = dof > T(1e6) + ? -T(0.91893853320467274178L)-T(0.25)/dof + : std::lgamma((dof+1)/2)-std::lgamma(dof/2) + -(std::log(dof)+T(1.14472988584940017414L))/2; + return normalizer-(dof+1)/2*tail; } - -template -typename std::pair stats::highestDensityInterval(const Container& c) { - using RealType = typename Container::value_type; - - const RealType p = 0.95f; - - // Build a sorted copy of the data - std::vector x(c); - std::sort(x.begin(), x.end()); - - // Choose a credible interval - size_t ciNumPoints = static_cast(std::floor(static_cast(x.size()) * p)); - std::pair minWidthCI = std::make_pair(x.front(), x.back()); - - for (size_t i = 0; i < x.size() - ciNumPoints; i++) { - RealType ciWidth = x[i + ciNumPoints] - x[i]; - if (ciWidth < minWidthCI.second - minWidthCI.first) { - minWidthCI = std::make_pair(x[i], x[i + ciNumPoints]); +template T StudentTPDF(T x,T dof) { return std::exp(StudentTLogPDF(x,dof)); } +template T Beta(T x,T y) { + detail::positive(x); detail::positive(y); + return std::exp(std::lgamma(x)+std::lgamma(y)-std::lgamma(x+y)); +} +template typename C::value_type mean(const C& c) { + detail::values(c); + long double sum=0; + for (auto x:c) sum+=static_cast(x)/c.size(); + return static_cast(sum); +} +template typename C::value_type stdev(const C& c) { + detail::values(c); + long double m=mean(c), scale=0, ss=1; + // Scaled sum of squares needs no temporary allocation and avoids overflow. + for (auto x:c) { + long double d=std::abs(static_cast(x)-m); + if (d!=0) { + if (scale(scale*std::sqrt(ss/c.size())); +} +template std::pair +highestDensityInterval(const C& c,double mass) { + detail::values(c); + if (!(mass>0 && mass<=1)) throw std::invalid_argument("mass must be in (0,1]"); + using T=typename C::value_type; + std::vector x(c.begin(),c.end()); + std::sort(x.begin(),x.end()); + size_t count=std::max(size_t(1),static_cast(std::ceil(mass*x.size()))); + size_t best=0; + for(size_t i=1;i+count<=x.size();++i) + if(static_cast(x[i+count-1])-x[i](x[best+count-1])-x[best]) best=i; + return {x[best],x[best+count-1]}; +} } diff --git a/tests/regression.cpp b/tests/regression.cpp new file mode 100644 index 0000000..a4ac942 --- /dev/null +++ b/tests/regression.cpp @@ -0,0 +1,92 @@ +#include "best.hpp" +#include +#include +#include +#include +#include + +int checks=0; +void check(bool ok,const char* message) { ++checks; if(!ok) throw std::runtime_error(message); } +void near(double x,double y,double tolerance,const char* message) { check(std::abs(x-y) void rejects(F f) { bool threw=false; try { f(); } catch(const std::exception&) { threw=true; } check(threw,"expected exception"); } +using V=std::vector; +using S=AMWG; +int main() { try { + near(stats::NormalPDF(0.,0.,1.),0.3989422804014327,1e-14,"normal reference"); + near(stats::StudentTPDF(0.,1.),1/3.141592653589793,1e-14,"Cauchy reference"); + near(stats::StudentTPDF(2.,1.),1/(5*3.141592653589793),1e-14,"Cauchy tail"); + near(stats::StudentTLogPDF(0.,1e100),stats::NormalLogPDF(0.,0.,1.),1e-12,"large dof normal limit"); + near(stats::Beta(2.,3.),1./12,1e-14,"beta reference"); + check(stats::ExponentialPDF(-1.,1.)==0,"exponential support"); + near(stats::NormalLogPDF(100.,0.,1.),-5000.918938533205,1e-10,"normal log tail"); + check(std::isfinite(stats::StudentTLogPDF(1e200,3.)),"Student log tail"); + near(stats::mean(V{1,2,3}),2,1e-14,"mean"); + near(stats::stdev(V{1,2,3}),std::sqrt(2./3),1e-14,"population sd"); + near(stats::stdev(V{1e100,-1e100})/1e100,1,1e-14,"large sd"); + check(stats::highestDensityInterval(std::list{4})==std::make_pair(4.,4.),"singleton HDI"); + V ordered; for(int i=0;i<20;++i) ordered.push_back(i); + check(stats::highestDensityInterval(ordered)==std::make_pair(0.,18.),"95 percent includes 19 of 20"); + rejects([]{stats::mean(V{});}); rejects([]{stats::stdev(V{NAN});}); + rejects([]{stats::highestDensityInterval(V{});}); rejects([]{stats::highestDensityInterval(V{1},0);}); + rejects([]{stats::StudentTPDF(0.,0.);}); rejects([]{S s(0);}); + S s(50,42); rejects([&]{s.Sample(1);}); rejects([&]{s.Burn(0);}); + rejects([&]{s.Init({{0}},{});}); rejects([&]{s.Init({{0}},[](const S::ParamArray&){return NAN;});}); + auto normal=[](const S::ParamArray& p){return -p[0]*p[0]/2;}; + { + s.Init({{0}},normal); s.Burn(1000); s.Sample(101); + check(s.chain().size()==101,"exact sample count"); + auto before=s.chain(); s.Burn(100); check(s.chain()==before,"burn preserves chain"); + s.Sample(0); check(s.chain()==before,"zero samples"); + check(s.NextSample()==1,"one proposal per coordinate"); + check(s.chain().back()==s.state(),"stores completed sweep"); + } + S a(50,12),b(50,12); a.Init({{0}},normal); b.Init({{0}},normal); + a.Burn(77); a.Sample(123); b.Burn(77); b.Sample(23); b.Sample(100); + check(a.chain()==b.chain(),"split-call reproducibility"); + a.Init({{0}},normal); a.Burn(77); a.Sample(123); check(a.chain()==b.chain(),"reinitialization resets adaptation and RNG"); + a.chain().clear(); a.Sample(50); b.Sample(50); + check(a.state()==b.state(),"adaptation independent of stored chain"); + S bounded; bounded.Init({{0}},[](const S::ParamArray& p){return p[0]==0 ? 0. : -INFINITY;}); bounded.Sample(100); + check(bounded.state()[0]==0,"reject out of support"); + S throwing; throwing.Init({{0}},[](const S::ParamArray& p)->double {if(p[0]!=0)throw std::runtime_error("callback");return 0;}); + rejects([&]{throwing.Sample(1);}); check(throwing.chain().empty(),"callback exception propagation"); + // Fixed seeds and generous tolerances: regression checks, not convergence proofs. + for(unsigned seed:{7u,42u,123u}) { + S n(50,seed); n.Init({{0}},normal); n.Burn(5000); n.Sample(60000); + V draws; for(auto p:n.chain()) draws.push_back(p[0]); + near(stats::mean(draws),0,0.06,"normal empirical mean"); near(stats::stdev(draws),1,0.06,"normal empirical sd"); + S e(50,seed); e.Init({{1}},[](const S::ParamArray& p){return p[0]<0?-INFINITY:-p[0];}); e.Burn(5000);e.Sample(60000); + draws.clear();for(auto p:e.chain())draws.push_back(p[0]); + near(stats::mean(draws),1,0.10,"exponential empirical mean"); + } + // Independent chains can be owned and run by separate threads. + V outcomes(2); std::thread t1([&]{S n;n.Init({{0}},normal);n.Sample(100);outcomes[0]=n.state()[0];}); + std::thread t2([&]{S n;n.Init({{0}},normal);n.Sample(100);outcomes[1]=n.state()[0];});t1.join();t2.join();check(outcomes[0]==outcomes[1],"independent chains"); + rejects([]{BEST m(V{},V{1});}); rejects([]{BEST m(V{1},V{1});}); + BEST model(V{1,2,3},V{3,4,5},50,42); + std::pair hdi;double mean; + rejects([&]{model.ComputeStats(hdi,mean);}); + check(!std::isfinite(model.LogPosterior({{2,4,1,1,0.5}})),"nu lower bound"); + check(!std::isfinite(model.LogPosterior({{2,4,1e10,1,5}})),"sigma upper bound"); + check(std::isfinite(model.LogPosterior({{2,4,1,1,5}})),"valid posterior"); + // Independent Cauchy (nu=1) formula for all factors of the joint model. + double sd=std::sqrt(5./3), priorSd=sd*1e6; + double expected=-std::log(29.)-2*std::log(sd*1000-sd/1000) + -2*std::log(priorSd*std::sqrt(2*3.141592653589793))-1/(priorSd*priorSd) + -6*std::log(3.141592653589793)-4*std::log(2.); + near(model.LogPosterior({{2,4,1,1,1}}),expected,1e-12,"joint posterior reference"); + BEST> floats({1,2,3},{3,4,5}); + near(floats.LogPosterior({{2,4,1,1,1}}),expected,1e-4,"float posterior precision"); + V original{1,2,3};BEST owned(original,V{3,4,5}); + auto density=owned.LogPosterior({{2,4,1,1,1}});original[0]=100; + check(owned.LogPosterior({{2,4,1,1,1}})==density,"model owns observations"); + // A well-identified fixture avoids treating a short diffuse chain as a + // convergence guarantee for two groups with only three observations. + V y1,y2; for(int i=0;i<40;++i) {double x=(i%5-2)*0.3; y1.push_back(x);y2.push_back(x+2);} + BEST identified(y1,y2,50,42); + identified.Burn(5000);identified.Sample(20000);identified.ComputeStats(hdi,mean); + check(identified.chain().size()==20000 && std::isfinite(mean) && hdi.first