Общие функции для стандартного отклонения в монте-карло C++

Я должен вычислить функцию стандартного отклонения в некоторых симуляциях Монте-Карло. Формула такая: введите здесь описание изображения

Я думаю, что мои результаты далеки от того, какими они должны быть. Моя функция использует кортежи из библиотеки boost и выглядит так:

double add_square(double prev_sum, double new_val)
{
  return prev_sum + new_val*new_val;
}

template <typename V>
double vec_add_squares(const V<double>& v)
{
  return std::accumulate(v.begin(), v.end(), 0.0, add_square);
}

    template <class T> 
    boost::tuple<double,double> get_std_dev_and_error(const vector<T>& input, double r, double N)
{
 double M = double(input.size());

 double sum = std::accumulate(input.begin(),input.end(),0.0);
 double Squared_sum = vec_add_squares(input);

 std::cout << "sum " << Squared_sum << endl;

 // Calls Sum
 double term1 = Squared_sum - (sum/M)*sum;

 double SD = (sqrt(term1) * exp(-2.0 * r *N))/(M-1) ;
 double SE = SD/sqrt(M);
 std::cout << "SD = " << SD << endl;
 std::cout << "SE = " << SE << endl;

 return boost::tuple<double,double>(SD, SE) ;
 }
  1. Может ли кто-нибудь увидеть здесь какие-либо ошибки?
  2. Кроме того, в библиотеке STL есть функция «накопить» - существует ли накопление в квадрате (члены контейнера)?

c++ generic-programming montecarlo
person Mathias    schedule 15.11.2012    source источник
comment
Для накопления вы можете написать свой собственный функтор   -  person PSIAlt    schedule 15.11.2012
comment
Спасибо за ответ PSIA. Я написал свой собственный, как в приведенном выше коде, но я получаю неверные результаты, поэтому он должен быть неправильным или неправильно реализованным в программе.   -  person Mathias    schedule 15.11.2012
comment
К вашему сведению, собственная ошибка индикатора ошибки может быть довольно большой. Получение доверительных интервалов для моделирования методом Монте-Карло сложнее, чем обычно считается.   -  person Alexandre C.    schedule 16.11.2012
comment
Да, он довольно большой... причина, по которой я знаю, что это неправильно, в том, что мне сообщили результаты.   -  person Mathias    schedule 16.11.2012
comment
Может быть, stats.stackexchange.com — хорошее место, чтобы спросить о формуле расчета?   -  person Riga    schedule 17.11.2012

Ответы (1)


arrow_upward
4
arrow_downward

Просто используйте Boost.Accumulators (поскольку вы уже используете boost ):

#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/range/algorithm.hpp>
#include <iostream>
#include <ostream>

using namespace boost;
using namespace boost::accumulators;
using namespace std;

int main()
{
    accumulator_set<double, stats<tag::sum , tag::variance, tag::mean > > acc;
    double data[] = {1., 2., 3.};
    acc = for_each(data, acc);
    cout << "sum = " << sum(acc) << endl;
    cout << "variance = " << variance(acc) << endl;
    cout << "sqrt(variance()) = " << sqrt(variance(acc)) << endl;
    cout << "mean = " << mean(acc) << endl;
}

Выход:

sum = 6
variance = 0.666667
sqrt(variance()) = 0.816497
mean = 2
person Evgeny Panasyuk    schedule 15.11.2012
comment
Мне очень нравится Boost.Accumulators - я слишком поторопился со стандартным отклонением, это не часть Accumulators (к сожалению)... Я добавлю изображение формулы - person Mathias; 16.11.2012
comment
@Mathias: у вас есть дисперсия, доступная в boost::accumulators. Просто возьмите его квадратный корень (не забудьте дополнительный множитель exp (-2rT) / (M - 1)). В любом случае, не используйте опубликованную вами формулу - во многих случаях она ужасно нестабильна (boost использует стабильную формулу on-line для дисперсии). - person Alexandre C.; 16.11.2012
comment
Когда я пытаюсь переписать уравнения из boost.org/doc/libs/1_52_0/doc/html/boost/accumulators/impl/ Я все еще не могу получить ту же формулу, что и выше (за исключением exp(-2Tr)) - person Mathias; 16.11.2012
comment
@АлександрК. извини забыл отметить тебя - person Mathias; 16.11.2012
comment
@Mathias: Ваша формула гласит прибл. SD = sqrt (дисперсия) * exp (-2rT) / sqrt (M - 1). Перед первой суммой должно быть M /(M - 1), а знаменатель должен быть M вместо M - 1, чтобы быть точно стандартной оценкой stdev, которую вы получили бы из boost::accumulators. Однако это не должно сильно изменить результат; - person Alexandre C.; 16.11.2012
comment
@АлександрК. Я знаю, но меня попросили закодировать именно эту формулу (иначе я бы использовал библиотеку boost)... Я немного изменил свой код - это делает его более стабильным: - person Mathias; 16.11.2012
comment
@АлександрК. Я внес некоторые изменения в свой код - теперь он стал лучше. Я должен использовать точную формулу (как указано), поэтому я не могу использовать функции повышения. - person Mathias; 16.11.2012
comment
@Mathias Матиас, значит, boost дает вам правильный ответ, но ваш код неверен? Если проблема связана с нестабильностью вычислений, попробуйте другую статистику с обеими формулами — вашей и форсированной. Если ответ будет таким же, то ваша формула конкретно плоха на исходном наборе данных. - person Riga; 17.11.2012
comment
@Riga, это хороший совет - спасибо. Я сравнил свои симуляции с результатами других людей, и оказалось, что в некоторых случаях я получаю одни и те же результаты, а в других — совершенно другие — это очень раздражает. - person Mathias; 17.11.2012