hops
PotentialScaleReductionFactor.hpp
Go to the documentation of this file.
1 #ifndef HOPS_POTENTIALSCALEREDUCTIONFACTOR_HPP
2 #define HOPS_POTENTIALSCALEREDUCTIONFACTOR_HPP
3 
4 #include <string>
5 #include <stdexcept>
6 #include <vector>
7 #include <cmath>
8 #include <cassert>
9 #include <memory>
10 
11 #include "IsConstantChain.hpp"
12 
13 namespace hops {
14  /*
15  * states.size() is the number of chains, states[i].size() is the number of draws, states[i][j].rows() is the dimensionality
16  *
17  *
18  */
19  template<typename StateType>
20  double computePotentialScaleReductionFactor (const std::vector<const std::vector<StateType>*>& chains,
21  unsigned long dimension,
22  unsigned long numUnseen,
23  std::vector<double>& sampleVariancesSeen,
24  std::vector<double>& intraChainExpectationsSeen,
25  double& interChainExpectationSeen,
26  unsigned long& numSeen) {
27  using Scalar = typename StateType::Scalar;
28 
29  unsigned long d = dimension;
30  unsigned long numChains = chains.size();
31  assert(numChains >= 2);
32 
33  assert(numChains == intraChainExpectationsSeen.size());
34  assert(numChains == sampleVariancesSeen.size());
35 
36  unsigned long numDraws = chains[0]->size();
37  assert(numDraws >= numUnseen);
38 #ifndef NDEBUG // prevent warning
39  for (auto& draws : chains) {
40  assert(numDraws == draws->size());
41  }
42 #endif // NDEBUG
43 
44  if (isConstantChain(chains, dimension)) {
45  return std::numeric_limits<double>::quiet_NaN();
46  }
47 
48  auto intraChainExpectationsUnseen = std::vector<double>(numChains, 0),
49  intraChainExpectations = std::vector<double>(numChains, 0);
50  double interChainExpectation,
51  interChainExpectationUnseen = 0,
52  eta = numSeen / double((numUnseen + numSeen));
53  for (unsigned long m = 0; m < numChains; ++m) {
54  for (unsigned long n = numDraws - numUnseen; n < numDraws; ++n) {
55  intraChainExpectationsUnseen[m] += (*chains[m])[n](d) / double(numUnseen);
56  }
57  intraChainExpectations[m] = eta * intraChainExpectationsSeen[m] + (1 - eta) * intraChainExpectationsUnseen[m];
58  interChainExpectationUnseen += intraChainExpectationsUnseen[m] / numChains;
59  }
60 
61  interChainExpectation = eta * interChainExpectationSeen + (1 - eta) * interChainExpectationUnseen;
62 
63  auto sampleVariances = std::vector<double>(numChains, 0);
64  Scalar betweenChainVariance = 0,
65  withinChainVariance = 0,
66  biasedVariance = 0,
67  kappa = (numUnseen + numSeen - 1);
68 
69  for (unsigned long m = 0; m < numChains; ++m) {
70  betweenChainVariance += (numSeen + numUnseen) * std::pow(intraChainExpectations[m] - interChainExpectation, 2) / (numChains - 1);
71  biasedVariance = 0;
72  for (unsigned long n = numDraws - numUnseen; n < numDraws; ++n) {
73  biasedVariance += std::pow((*chains[m])[n](d) - intraChainExpectations[m], 2);
74  }
75  sampleVariances[m] =
76  ( (numSeen - 1) * sampleVariancesSeen[m] +
77  numSeen * std::pow(intraChainExpectationsSeen[m] - intraChainExpectations[m], 2) +
78  biasedVariance ) / kappa ;
79 
80  withinChainVariance += sampleVariances[m] / numChains;
81  }
82 
83  Scalar rHat = std::sqrt(((numUnseen + numSeen - 1) + betweenChainVariance / withinChainVariance) / (numSeen + numUnseen));
84 
85  // update seen values
86  sampleVariancesSeen = sampleVariances;
87  intraChainExpectationsSeen = intraChainExpectations;
88  interChainExpectationSeen = interChainExpectation;
89  numSeen = numSeen + numUnseen;
90 
91  return rHat;
92  }
93 
94  /*
95  * states.size() is the number of chains, states[i].size() is the number of draws, states[i][j].rows() is the dimensionality
96  *
97  *
98  */
99  template<typename StateType>
100  double computePotentialScaleReductionFactor (const std::vector<std::vector<StateType>>& chains,
101  unsigned long dimension,
102  unsigned long numUnseen,
103  std::vector<double>& sampleVariancesSeen,
104  std::vector<double>& intraChainExpectationsSeen,
105  double& interChainExpectationSeen,
106  unsigned long& numSeen) {
107  std::vector<const std::vector<StateType>*> chainsPtrArray;
108  for (auto& chain : chains) {
109  chainsPtrArray.push_back(&chain);
110  }
111  return computePotentialScaleReductionFactor(chainsPtrArray,
112  dimension,
113  numUnseen,
114  sampleVariancesSeen,
115  intraChainExpectationsSeen,
116  interChainExpectationSeen,
117  numSeen);
118  }
119 
120  /*
121  * Compute PSRF for all dimensions at once.
122  *
123  *
124  */
125  template<typename StateType>
126  std::vector<double> computePotentialScaleReductionFactor (const std::vector<const std::vector<StateType>*>& chains,
127  unsigned long numUnseen,
128  std::vector<std::vector<double>>& sampleVariancesSeen,
129  std::vector<std::vector<double>>& intraChainExpectationsSeen,
130  std::vector<double>& interChainExpectationSeen,
131  unsigned long& numSeen) {
132  unsigned long dimensions = (*chains[0])[0].rows();
133  unsigned long numChains = chains.size();
134 
135  if (numChains <= 1) {
136  throw std::runtime_error(std::string("PSRF is only defined for at least two chains, your data consists of only ") + std::to_string(numChains));
137  }
138 
139  // initialize empty vectors correctly
140  if (sampleVariancesSeen.size() == 0) {
141  sampleVariancesSeen = std::vector<std::vector<double>>(dimensions,
142  std::vector<double>(numChains, 0));
143  }
144 
145  if (intraChainExpectationsSeen.size() == 0) {
146  intraChainExpectationsSeen = std::vector<std::vector<double>>(dimensions,
147  std::vector<double>(numChains, 0));
148  }
149 
150  if (interChainExpectationSeen.size() == 0) {
151  interChainExpectationSeen = std::vector<double>(dimensions, 0);
152  }
153 
154  assert(dimensions == sampleVariancesSeen.size());
155  assert(dimensions == intraChainExpectationsSeen.size());
156  assert(dimensions == interChainExpectationSeen.size());
157 
158  unsigned long numSeenUpdated;
159 
160  std::vector<double> rhats = std::vector<double>(dimensions);
161  for (unsigned long dimension = 0; dimension < dimensions; ++dimension) {
162  unsigned long _numSeen = numSeen;
163  rhats[dimension] = computePotentialScaleReductionFactor(chains,
164  dimension,
165  numUnseen,
166  sampleVariancesSeen[dimension],
167  intraChainExpectationsSeen[dimension],
168  interChainExpectationSeen[dimension],
169  _numSeen);
170  numSeenUpdated = _numSeen;
171  }
172 
173  numSeen = numSeenUpdated;
174 
175  return rhats;
176  }
177 
178  /*
179  * Compute PSRF for all dimensions at once.
180  *
181  *
182  */
183  template<typename StateType>
184  std::vector<double> computePotentialScaleReductionFactor (const std::vector<std::vector<StateType>>& chains,
185  unsigned long numUnseen,
186  std::vector<std::vector<double>>& sampleVariancesSeen,
187  std::vector<std::vector<double>>& intraChainExpectationsSeen,
188  std::vector<double>& interChainExpectationSeen,
189  unsigned long& numSeen) {
190  std::vector<const std::vector<StateType>*> chainsPtrArray;
191  for (auto& chain : chains) {
192  chainsPtrArray.push_back(&chain);
193  }
194  return computePotentialScaleReductionFactor(chainsPtrArray,
195  numUnseen,
196  sampleVariancesSeen,
197  intraChainExpectationsSeen,
198  interChainExpectationSeen,
199  numSeen);
200  }
201 
202  template<typename StateType>
203  double computePotentialScaleReductionFactor (const std::vector<const std::vector<StateType>*>& chains,
204  unsigned long dimension) {
205  std::vector<double> sampleVariances(chains.size());
206  std::vector<double> intraChainExpectations(chains.size());
207  double interChainExpectation = 0;
208  unsigned long numSeen = 0;
209  return computePotentialScaleReductionFactor(chains,
210  dimension,
211  chains[0]->size(),
212  sampleVariances,
213  intraChainExpectations,
214  interChainExpectation,
215  numSeen);
216  }
217 
218  template<typename StateType>
219  double computePotentialScaleReductionFactor (const std::vector<std::vector<StateType>>& chains,
220  unsigned long dimension) {
221  std::vector<const std::vector<StateType>*> chainsPtrArray;
222  for (auto& chain : chains) {
223  chainsPtrArray.push_back(&chain);
224  }
225  return computePotentialScaleReductionFactor(chainsPtrArray, dimension);
226  }
227 
228  template<typename StateType>
229  std::vector<double> computePotentialScaleReductionFactor (const std::vector<const std::vector<StateType>*>& chains) {
230  unsigned long dimensions = (*chains[0])[0].rows();
231  std::vector<double> rhats = std::vector<double>(dimensions);
232  for (unsigned long dimension = 0; dimension < dimensions; ++dimension) {
233  rhats[dimension] = computePotentialScaleReductionFactor(chains, dimension);
234  }
235  return rhats;
236  }
237 
238  template<typename StateType>
239  std::vector<double> computePotentialScaleReductionFactor (const std::vector<std::vector<StateType>>& chains) {
240  std::vector<const std::vector<StateType>*> chainsPtrArray;
241  for (auto& chain : chains) {
242  chainsPtrArray.push_back(&chain);
243  }
244  return computePotentialScaleReductionFactor(chainsPtrArray);
245  }
246 }
247 
248 
249 #endif //HOPS_POTENTIALSCALEREDUCTIONFACTOR_HPP
250 
IsConstantChain.hpp
hops
Definition: CsvReader.hpp:8
string
NAME string(REPLACE ".cpp" "_bin" example_name ${example_filename}) if($
Definition: hops/Third-party/HighFive/src/examples/CMakeLists.txt:6
hops::computePotentialScaleReductionFactor
double computePotentialScaleReductionFactor(const std::vector< const std::vector< StateType > * > &chains, unsigned long dimension, unsigned long numUnseen, std::vector< double > &sampleVariancesSeen, std::vector< double > &intraChainExpectationsSeen, double &interChainExpectationSeen, unsigned long &numSeen)
Definition: PotentialScaleReductionFactor.hpp:20
hops::isConstantChain
bool isConstantChain(const std::vector< const std::vector< StateType > * > &chains, long dimension)
Definition: IsConstantChain.hpp:25
Generated by   doxygen 1.8.17