Гауссово случайное распределение в Postgresql

У меня есть таблица, скажем, 250 URL:

create table url (
  id serial,
  url varchar(64)
)

Эти URL-адреса соответствуют каждому веб-сайту. Каждый из сайтов имеет разную популярность. Допустим, что id=125 (тот, который сосредоточен на гауссиане) является самым популярным, те в id=1 или же id=250 наименее популярны.

Я хочу заполнить таблицу "log", подобную следующей, значением url среди тех, что указаны в таблице "url", но с учетом того, что различные URL-адреса могут появляться чаще (например, для URL-адреса, идентификатор которого равен 125, будет самый популярный).

create table log (
  id serial,
  url_id integer
)

Я хочу избежать использования random() так как оно однородно и не очень "реально".

Как этого можно достичь с помощью Postgresql?

6 ответов

Решение

Сумма из 12 равномерных распределений в диапазоне [0, 1) является хорошим приближением к гауссову распределению, ограниченному в диапазоне [0, 12). Затем его можно легко перемасштабировать, умножив на константу, а затем добавив / вычтя константу.

select
    random() + 
    random() + 
    random() +
    random() + 
    random() + 
    random() +
    random() + 
    random() + 
    random() +
    random() + 
    random() + 
    random();

http://books.google.com/books?id=EKA-yeX2GVgC&pg=PA185&lpg=PA185&dq=%22sum+of+12+uniform+random+variables%22&source=bl&ots=YfwwE0fBB3&sig=HX9J9Oe6x316kVL8uamDU_GOsn4&hl=en&sa=X&ei=bJLZUur1GozaqwGHm4DQDQ&ved=0CEUQ6AEwAw

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

Существует, по крайней мере, PostGreSQL 8.4, дополнительный модуль под названием tablefunc ( http://www.postgresql.org/docs/9.2/static/tablefunc.html).

Он предлагает функцию normal_rand(n, mean, stddev), генерирующую n псевдослучайных чисел с использованием гауссовского распределения (поэтому эта функция возвращает набор значений, обычно используемых в предложении FROM). Однако, если вы установите n равным 1, его можно использовать как функцию, возвращающую значение, а не набор значений.

Рассматривая таблицу nb10, содержащую 10 записей, два следующих запроса возвращают набор из 10 псевдослучайных чисел после стандартного распределения Гаусса (среднее = 0, stddev = 1)

SELECT normal_rand(1, 0, 1) FROM nb10;

а также

SELECT * from normal_rand(10, 0, 1);

Я надеюсь, что это может помочь любому в будущем...:-)

Чтобы конкретно ответить на ваш вопрос, вы можете использовать что-то вроде:

SELECT floor(random_rand(1, 0, 1) * 250 + 125);

К сожалению, с помощью этого запроса можно получить ответ не в диапазоне [0, 249]. Вы могли бы, например:

  • использовать рекурсивный запрос, который я нахожу немного излишним, для отбрасывания значений не в диапазоне [0, 249], или
  • сделайте свой выбор в цикле на вашем главном языке, принимая значение только в том случае, если оно находится в диапазоне [0, 249], или
  • используйте оператор по модулю, чтобы остаться в диапазоне [0, 250[, я думаю, что это лучшее решение, хотя оно слегка изменяет гауссову кривую. Вот последний запрос, который я предлагаю вам использовать (трюки по модулю /+/modulo заключаются в том, что -x по модулю y с положительным числом x дает отрицательное число в PostGreSQL, что неплохо:p):

    SELECT ((floor(normal_rand(1,0,1)*250 + 125)::int % 250) + 250) % 250 as v;
    

Простой факт заключается в том, что вы хотите создать свою собственную функцию, которая упаковывает rand() во что-то, что обеспечивает гауссово распределение неявно или явно.

У меня нет статистических данных, чтобы рассказать вам, как преобразовать равномерное распределение в гауссовское, но вам придется написать конвертер. Примерно так, как указано на http://www.perlmonks.org/?node_id=26889 (если вам не нравится Perl, вы можете переписать это в pl/pgsql или даже в простом SQL).

CREATE OR REPLACE FUNCTION gaussian_rand() RETURNS numeric LANGUAGE PLPERL VOLATILE AS
$$
    my ($u1, $u2);  # uniformly distributed random numbers
    my $w;          # variance, then a weight
    my ($g1, $g2);  # gaussian-distributed numbers

    do {
        $u1 = 2 * rand() - 1;
        $u2 = 2 * rand() - 1;
        $w = $u1*$u1 + $u2*$u2;
    } while ( $w >= 1 );

    $w = sqrt( (-2 * log($w))  / $w );
    $g2 = $u1 * $w;
    $g1 = $u2 * $w;
    # return both if wanted, else just one
    return $g1;

$$;

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

SELECT normal_rand(1, 0, 1); -- generates 1 single value with mean 0 and a standard deviation of 1

Приведенный выше запрос должен генерировать одно значение в нормальном распределении

Если он не установлен, попробуйте это:

CREATE EXTENSION "tablefunc";

В противном случае вам нужно будет войти как суперпользователь и установить модуль.

Вы также можете реализовать это непосредственно на встроенном языке PL/PgSQL.

create or replace
function random_gauss( avg real = 0, stddev real = 1 )
returns real language plpgsql as $$
declare x1 real; x2 real; w real;
begin
  loop
    x1 = 2.0 * random() - 1.0;
    x2 = 2.0 * random() - 1.0;
    w = x1*x1 + x2*x2;
    exit when w < 1.0;
  end loop;
  return avg + x1 * sqrt(-2.0*ln(w)/w) * stddev;
end; $$;

with data as (
  select t, random_gauss(100,15)::integer score from generate_series(1,1000000) t
)
select
  score,
  sum(1),
  repeat('=',sum(1)::integer/500) bar
from data
where score between 60 and 140
group by score
order by 1;

rollback;

Это дает нам что-то подобное из выборки из 1 миллиона чисел со средним значением 100 и стандартным отклонением 15.

 score |  sum  |                          bar
-------+-------+--------------------------------------------------------
    60 |   764 | =
    61 |   893 | =
    62 |  1059 | ==
    63 |  1269 | ==
    64 |  1524 | ===
    65 |  1740 | ===
    66 |  1990 | ===
    67 |  2346 | ====
    68 |  2741 | =====
    69 |  3160 | ======
    70 |  3546 | =======
    71 |  4109 | ========
    72 |  4633 | =========
    73 |  5252 | ==========
    74 |  5952 | ===========
    75 |  6536 | =============
    76 |  7429 | ==============
    77 |  8140 | ================
    78 |  9061 | ==================
    79 | 10063 | ====================
    80 | 10844 | =====================
    81 | 11911 | =======================
    82 | 13180 | ==========================
    83 | 13880 | ===========================
    84 | 15111 | ==============================
    85 | 16016 | ================================
    86 | 17310 | ==================================
    87 | 18262 | ====================================
    88 | 19615 | =======================================
    89 | 20400 | ========================================
    90 | 21186 | ==========================================
    91 | 22190 | ============================================
    92 | 23103 | ==============================================
    93 | 24150 | ================================================
    94 | 24327 | ================================================
    95 | 24992 | =================================================
    96 | 25505 | ===================================================
    97 | 25868 | ===================================================
    98 | 26146 | ====================================================
    99 | 26574 | =====================================================
   100 | 27104 | ======================================================
   101 | 26599 | =====================================================
   102 | 26345 | ====================================================
   103 | 25940 | ===================================================
   104 | 25485 | ==================================================
   105 | 25157 | ==================================================
   106 | 24827 | =================================================
   107 | 23844 | ===============================================
   108 | 23262 | ==============================================
   109 | 22211 | ============================================
   110 | 21326 | ==========================================
   111 | 20315 | ========================================
   112 | 19496 | ======================================
   113 | 18026 | ====================================
   114 | 17182 | ==================================
   115 | 16026 | ================================
   116 | 14979 | =============================
   117 | 13959 | ===========================
   118 | 12840 | =========================
   119 | 11718 | =======================
   120 | 11169 | ======================
   121 | 10037 | ====================
   122 |  9273 | ==================
   123 |  8041 | ================
   124 |  7402 | ==============
   125 |  6761 | =============
   126 |  5827 | ===========
   127 |  5257 | ==========
   128 |  4736 | =========
   129 |  4153 | ========
   130 |  3494 | ======
   131 |  3103 | ======
   132 |  2731 | =====
   133 |  2379 | ====
   134 |  2064 | ====
   135 |  1696 | ===
   136 |  1481 | ==
   137 |  1246 | ==
   138 |  1024 | ==
   139 |   910 | =
   140 |   788 | =

Для таких людей, как я, которые пришли сюда, чтобы попробовать истинное непрерывное распределение Гаусса в PostgreSQL, вот как это сделать:

      SELECT cos(2*pi()*random()) * sqrt(-2*ln(random())) AS my_standard_normal_random_variable

Объяснение: это состоит из выборки полярных координат двумерной гауссианы (угол следует равномерному распределению по[0, 2π], а CDF радиуса r определяется выражениемr -> 1 - exp(-0.5 * r^2)соответствующий обратному CDFq -> sqrt(-2 * ln(1 - q))), затем проецируется на координату x с помощью косинуса, который работает, потому что координата x стандартного гауссова 2D является стандартной гауссовой 1D. Эта стратегия известна как преобразование Бокса-Мюллера .

Проверить: давайте посмотрим на некоторую сводную статистику:

      WITH my_gaussian_samples AS (
  SELECT cos(2*pi()*random()) * sqrt(-2*ln(random())) AS x
  FROM generate_series(1,1000000)
) SELECT 
  AVG(x) AS mean,
  AVG(x^2) AS standard_deviation, 
  COUNT((x > -1.96 AND x < 1.96) OR NULL)::float/COUNT(*) AS in_95_interval
FROM my_gaussian_samples;

          mean          | standard_deviation | in_95_interval 
------------------------+--------------------+----------------
 -0.0012901731683899918 | 0.9988458311780355 |       0.950221
Другие вопросы по тегам