Как наиболее эффективно предотвратить то, чтобы моя нормально распределенная случайная величина была нулевой? - PullRequest
Новая
       73

Как наиболее эффективно предотвратить то, чтобы моя нормально распределенная случайная величина была нулевой?

2 голосов
/ 19 июня 2011

Я пишу алгоритм Монте-Карло, в котором в какой-то момент мне нужно разделить случайную величину. Точнее: случайная величина используется в качестве ширины шага для разностного отношения, поэтому я фактически сначала умножаю что-то на переменную, а затем снова делю ее на некоторую локально линейную функцию этого выражения. Как 

double f(double);

std::tr1::variate_generator<std::tr1::mt19937, std::tr1::normal_distribution<> >
  r( std::tr1::mt19937(time(NULL)),
     std::tr1::normal_distribution<>(0) );

double h = r();
double a = ( f(x+h) - f(x) ) / h;

Это прекрасно работает большую часть времени, но не работает, когда h=0. Математически, это не проблема, потому что при любом конечном (или даже счетном) выборе нормально распределенных случайных величин все они будут отличны от нуля с вероятностью 1. Но в цифровой реализации я столкнусь с h==0 каждые ≈ Вызовы функций 2³² (независимо от того, имеет ли мерсенновый твистер период больше, чем вселенная, он все равно выдает обычные long с!).

Довольно просто избежать этой проблемы, в данный момент я делаю 

double h = r();
while (h==0) h=r();

но я не считаю это особенно элегантным. Есть ли лучший способ?

Функция, которую я оцениваю, на самом деле не просто ℝ-> ℝ, как f, а ℝᵐxℝⁿ -> which, в которой я вычисляю градиент по ℝᵐ переменным, в то же время численно интегрируя по ℝⁿ переменным. Вся функция наложена на непредсказуемый (но «когерентный») шум, иногда с определенными (но неизвестными) выдающимися частотами, вот что доставляет мне неприятности, когда я пробую его с фиксированными значениями для h.

Ответы [ 5 ]

3 голосов
/ 19 июня 2011

Отношение двух нормально распределенных случайных величин является распределением Коши.Распределение Коши является одним из тех неприятных распределений с бесконечной дисперсией.Очень противный действительно.Распределение Коши запутает ваш эксперимент Монте-Карло.

Во многих случаях, когда вычисляется отношение двух случайных величин, знаменатель не является нормальным.Люди часто используют нормальное распределение для аппроксимации этой ненормально распределенной случайной величины, потому что с

  • нормальными распределениями обычно так легко работать,
  • обычно имеют такие хорошие математические свойства,
  • нормальное предположение кажется более или менее правильным, а
  • реальным распределением является медведь.

Предположим, вы делите на расстояние.Расстояние является полу-положительно определенным по определению и часто положительно определенным как случайная величина.Таким образом, расстояние от биты никогда не может быть нормально распределено.Тем не менее, люди часто принимают нормальное распределение для расстояния в тех случаях, когда среднее значение намного, намного больше, чем стандартное отклонение.Когда это нормальное предположение сделано, вам нужно защищаться от этих нереальных ценностей.Одно простое решение - усеченная нормаль.

3 голосов
/ 19 июня 2011

ваш путь кажется достаточно элегантным, может быть, немного другим: 

do {
    h = r();
} while (h == 0.0);
0 голосов
/ 19 июня 2011

Функция (f (x + h) -f (x)) / h имеет предел как h-> 0, и поэтому, если вы встретите h == 0, вы должны использовать этот предел.Пределом будет f '(x), поэтому, если вы знаете производную, вы можете использовать ее.

Если то, что вы на самом деле делаете, - это создание числа дискретных точек, хотя это приблизительно соответствует нормальному распределению, и этого достаточнодля вашего дистрибутива создайте его так, чтобы ни один из них не имел значения 0.

0 голосов
/ 19 июня 2011

В зависимости от того, что вы пытаетесь вычислить, возможно, что-то вроде этого будет работать: 

double h = r();
double a;
if (h != 0)
    a = ( f(x+h) - f(x) ) / h;
else
    a = 0;

Если f - линейная функция, это должно (я думаю?) Оставаться непрерывным при h =0.

Возможно, вы захотите вместо этого отловить исключения деления на ноль, чтобы избежать затрат на ветку. Обратите внимание, что это может иметь или не иметь пагубного влияния на производительность. - тест в обоих направлениях!

В Linux вам нужно будет построить файл, содержащий ваше потенциальное деление на ноль, с -fnon-call-exceptions и установите обработчик SIGFPE: 

struct fp_exception { };

void sigfpe(int) {
  signal(SIGFPE, sigfpe);
  throw fp_exception();
}

void setup() {
  signal(SIGFPE, sigfpe);
}

// Later...
    try {
        run_one_monte_carlo_trial();
    } catch (fp_exception &) {
        // skip this trial
    }

В Windows используйте SEH: 

__try 
{ 
    run_one_monte_carlo_trial();
} 
__except(GetExceptionCode() == EXCEPTION_INT_DIVIDE_BY_ZERO ? 
         EXCEPTION_EXECUTE_HANDLER : EXCEPTION_CONTINUE_SEARCH)
{ 
    // skip this trial
}

Это дает преимущество, заключающееся в том, что оно потенциально меньше влияет на быстрый путь.Ветви нет, хотя может быть некоторая корректировка записей обработчика исключений.В Linux может быть небольшое снижение производительности из-за компилятора, генерирующего более консервативный код для -fnon-call-exceptionsС меньшей вероятностью это станет проблемой, если код, скомпилированный в -fnon-call-exceptions, не выделяет никаких автоматических (стековых) объектов C ++.Стоит также отметить, что в этом случае деление на ноль делает ОЧЕНЬ дорогим.

0 голосов
/ 19 июня 2011

Если вы хотите сохранить нормальное распределение, вы должны либо исключить 0, либо присвоить 0 новому ранее не встречающемуся значению. Поскольку второе, скорее всего, невозможно в ограниченных пределах информатики, первое - наш единственный вариант.

...