Гауссово случайное распределение в 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();
Я искал способ генерировать числа в соответствии с гауссовым распределением и впервые нашел этот пост. Вот почему я делюсь тем, что нашел сразу после:
Существует, по крайней мере, 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