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

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

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 s!).

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

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

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

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

Источник

5 ответов

Решение

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

do {
    h = r();
} while (h == 0.0);
Источник

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

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

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

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

Источник

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

Источник

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

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

Источник

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

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++. Стоит также отметить, что в этом случае деление на ноль происходит ОЧЕНЬ дорого.

Источник
Другие вопросы по тегам