Я заинтересован в вычислении моментов высшего порядка некоторых наблюдаемых в программе на C ++. В настоящее время у меня есть рабочий метод, но он очень медленный, и моя реализация, вероятно, была наивной. Под «моментом» я подразумеваю только то, что обозначено здесь . У меня есть следующая функция, которая вычисляет момент до для усреднения (что я делаю «онлайн»),
vector<double> ExpVec(const vector<double>& Data, const double& momentOrder) {
vector<double> ExpedVec(Data.size());
for (size_t i = 0; i < Data.size(); ++i)
ExpedVec[i] = std::pow(Data[i], momentOrder);
return ExpedVec;
}
Каждый элемент ExpedVec затем передается в объект, который вычисляет скользящее среднее. Таким образом, Data.size() - это количество выборок, и каждый раз, когда обнаруживается новый Data, мы пропускаем его через вышеупомянутое и добавляем новый материал к скользящему среднему. Процедура усреднения в основном то, что можно найти здесь здесь . Таким образом, используя класс RunningStats, мы имеем что-то вроде
int main() {
uint nMoments = 5;
double stepSize = 0.5; // Consider half integer moments too
uint nSamples = 1000;
vector<vector<RunningStats>> RS(nMoments, vector<RunningStats>(nSamples));
for (uint n = 0; n < nSamples; ++n) {
vector<double> Data = getData(n); // A vector of length M say
for (size_t j = 0; j < nMoments; ++j) {
vector<double> Exped = ExpVec(Data, (j+1) * stepSize);
for (size_t k = 0; k < Exped.size(); ++k) {
RS[j][k].Push(Exped[j]); // Details aren't super important, it just stores the running average for each Data element
}
}
}
}
Следует отметить, что да, схема RunningStats класса здесь уже вычисляет до момента порядка 4. Однако, если ваш общий размер выборки nSamples очень большой и вы можете вычислять подмножества общего количества выборок параллельно, у вас есть проблема с «объединением» ваших отклонений ». С другой стороны, когда мы вычисляем средние значения параллельно, мы имеем среднее значение средних значений, и нет необходимости объединять дисперсию , перекос и c.
в любом В случае, если детали main, я думаю, не так важны. Кажется, главная проблема с производительностью - это (используя ответ, обозначенный здесь ) pow функция. Это обсуждалось здесь и здесь и кажется, что если требуются нецелые степени, std::pow - это путь к go.
Если есть более подходящие способы решения этой проблемы go, и если у людей есть явный опыт решения проблемы, тогда я приветствую любые комментарии. С другой стороны, с чисто программной точки зрения c, что делать? Так много вложенных циклов, что я не уверен, как их избежать.
РЕДАКТИРОВАТЬ: Работоспособный пример ниже, с помощью комментариев:
#include <vector>
#include <iostream>
#include <math.h>
#include <algorithm>
#include <random>
#include <iterator>
#include <functional>
#include <chrono>
#include <iomanip>
using namespace std;
using namespace std::chrono;
struct RunningAverage {
public:
RunningAverage() : nElems(0), mean(0.0) {}
void Push(const double& x) {
nElems++;
mean += (x - mean) / nElems;
}
double getMean() const { return mean; }
private:
double mean;
int nElems;
};
void ExpVec(vector<double>& Vec, const double& momentOrder) {
for (size_t i = 0; i < Vec.size(); ++i)
Vec[i] = std::pow(Vec[i], momentOrder);
}
void MultVecbyVec(vector<double>& Vec1, const vector<double>& Vec2) {
for (size_t i = 0; i < Vec1.size(); ++i)
Vec1[i] = Vec1[i] * Vec2[i];
}
int main() {
uint nMoments = 20;
uint nSamples = 1000000; // No. times to collect data
uint DataSize = 1; // Not an actual use case
double momentStep = 0.5; // nMoments/momentStep = largest moment
uniform_real_distribution<double> unif(0.0, 1.0); // Used for generating dummy data
default_random_engine re;
// With recycling
high_resolution_clock::time_point t0 = high_resolution_clock::now();
vector<vector<RunningAverage>> Moments1(nMoments, vector<RunningAverage>(DataSize));
for (uint i = 0; i < nSamples; ++i) {
vector<double> Data(DataSize);
for (uint i = 0; i < DataSize; ++i)
Data[i] = unif(re);
vector<double> Accu = Data;
ExpVec(Accu, momentStep);
vector<double> Base = Accu;
for (uint j = 0; j < nMoments; ++j) {
for (uint k = 0; k < DataSize; ++k)
Moments1[j][k].Push(Accu[k]);
MultVecbyVec(Accu, Base);
}
}
high_resolution_clock::time_point t1 = high_resolution_clock::now();
auto duration1 = duration_cast<milliseconds>(t1 - t0).count();
// Without recycling
high_resolution_clock::time_point t2 = high_resolution_clock::now();
vector<vector<RunningAverage>> Moments2(nMoments, vector<RunningAverage>(DataSize));
for (uint i = 0; i < nSamples; ++i) {
vector<double> Data(DataSize);
for (uint i = 0; i < DataSize; ++i)
Data[i] = unif(re); // Dummy data
for (size_t j = 0; j < nMoments; ++j) {
vector<double> Exped = Data;
ExpVec(Exped, (j + 1.) * momentStep);
for (size_t k = 0; k < Exped.size(); ++k)
Moments2[j][k].Push(Exped[k]);
}
}
high_resolution_clock::time_point t3 = high_resolution_clock::now();
auto duration2 = duration_cast<milliseconds>(t3 - t2).count();
// Summary
for (uint j = 0; j < nMoments; ++j)
for (uint k = 0; k < DataSize; ++k)
cout << "Moment order (" << fixed << setprecision(3) << (j + 1) * momentStep << "): With recycling = " << Moments1[j][k].getMean()
<< ". Without recycling = " << Moments2[j][k].getMean()
<< ". Expected value = " << 1./(momentStep * (j + 1.) + 1.) << endl;
cout << endl << "Run time with recycling: " << duration1 << endl;
cout << "Run time without recycling: " << duration2 << endl;
return 0;
}
Что дает выход:
Moment order (0.500): With recycling = 0.667. Without recycling = 0.667. Expected value = 0.667
Moment order (1.000): With recycling = 0.500. Without recycling = 0.500. Expected value = 0.500
Moment order (1.500): With recycling = 0.400. Without recycling = 0.400. Expected value = 0.400
Moment order (2.000): With recycling = 0.334. Without recycling = 0.333. Expected value = 0.333
.
.
.
Moment order (9.500): With recycling = 0.095. Without recycling = 0.095. Expected value = 0.095
Moment order (10.000): With recycling = 0.091. Without recycling = 0.091. Expected value = 0.091
Run time with recycling: 1056
Run time without recycling: 4124