Как генерировать случайные числа распределенные по...

yurimedvedev

степенному закону?
Я посмотрел в книжках и пришел к выводу, что нужно делать так:
1. Вводим ф-цию распределения y = x^a
2. Находим обратную ф-цию x = y^(1/a)
3. Генерируем число r, равномерно распределенное на [0:1]
4. Утверждаем, что x = r^(1/a) - есть случайное число, распределенное по нужному нам закону.
Так ли это? Проблема в том, что я написал прогу, которая должна генерить случайное число E, распределенное по закону f(E) ~ E^(-2). Однако, когда я сгенерил массив, нашел зависимость вероятности от E, она оказалась степенной с другим показателем (примерно 2.5). Где ошибка?
 
#!c:/perl/bin/perl -w
my ($N);
my ($r, $i, $E, $alpha);

$N = 100000;
$alpha = -2;

for($i=0; $i<$N; $i++){
$r = rand;
if($r==0){
next;
}
$E = exp( log($r) * 1.0/$alpha );
print($E, "\n");
}

yurimedvedev

Спасибо, но это не совсем то, что мне нужно.
Сейчас проблема носит характер академический. Мне нужно уметь генерировать любые распределения. Например, степенное. Т.е. проблема ближе к теорверу чем к прикладному программированию. Хотя результатом работы должна являться программа.

parfum74

alpha = -2
Вводим ф-цию распределения y = x^a
Ты уверен, что это можно назвать функцией распределения? (Напомню, что любая ф.р. F(x) при x->infty стремится к 1, при x->-infty стремится к 0.)

yurimedvedev

http://www.xycoon.com/pow_distribution.htm
С определениями трудности...

parfum74

Обрати внимание, что a>0 и попробуй протестить свою прогу снова

yurimedvedev

Хм... Я неправильно обрабатывал получившийся массив.
Чтобы узнать вид распределения, я разбивал множество значений величины E на отрезки и смотрел сколько раз у меня выпадало E в каждом отрезке, потом считал среднее значение. В моем случае Е принимало значения от 1 до ~500. Если я разбивал это на отрезки постоянной длины, то получал неправильный наклон (например 2.5 или 5). Если разбивал на отрезки в логарифмическом масштабе, то получал наклон ~2.04, ~2.01.
Что-то я не понял почему...

#!c:/perl/bin/perl -w

my ($E, @flux, $N, $str);
my ($minE, $maxE, $NE);
my ($first);

my @Scale;

@Scale = (1, 3, 5, 10, 30, 50, 100, 300, 500);
#@Scale = (0, 50, 100, 150, 200, 250, 300, 350, 400);

$NE = @Scale-1;
$N = 0;
$minE = 1;
$maxE = 100;

open(IN, "r.txt"); #r.txt - файл со случайными числами
while($E = <IN>){
for($i=0; $i<$NE; $i++){
if($E>= $Scale[$i] && $E<$Scale[$i+1]){
$flux[$i]++;
$N++;
}
}
}
for($i=0; $i<$NE; $i++){
$flux[$i] /= $N;
print$Scale[$i+1]+$Scale[$i])/2, "\t", $flux[$i], "\n");
}
close(IN);

Anastasia85

У нас на войне это было - доказывалось, что такой способ действительно дает нужное распределение. Писали проги, которые это делали.
Называется "метод обратной функции".
Только правильно тут сказали: не перепутаны ли тут функция распределения и плотность с.в.?
Ботать войну тут:
http://old.cmc-msu.ru/files05.html

ereyzer

Все так, только для того, чтобы твоя "функция распределения" была настоящей ф-ей распределения, нужно еще сделать нормировку: N=\int_Xmin^Xmax f(x)dx
Сначала для общего случая:
Пусть мы хотим сгенерировать случайную величину XX с плотностью вероятности y=f(x)/N на промежутке от Xmin до Xmax. Тогда ф-я распределения будет такой:
F(x) = \int_Xmin^x f(x)dx / N, а XX = F^{-1}(r 0<r<=1.
Если все проинтегрировать, то для степенной функции y=x^a получим:
XX = (r*(Xmax^{a+1}-Xmin^{a+1}) + Xmin^{a+1})^{1/{a+1}}
Когда ты будешь рисовать распределение по log10(x то наклон из интегрального превратится в дифференциальный. Как это показать я сходу не помню - можем вечером вместе подумать...
Оставить комментарий
Имя или ник:
Комментарий: