Таким образом, мне кажется, что для генерации заданного распределения вероятности мне нужна квантильная функция , обратная
совокупная функция распределения , как говорит @dmckee.
Возникает вопрос: как лучше всего сгенерировать и сохранить квантильную функцию, описывающую данную непрерывную гистограмму? У меня есть ощущение, что ответ будет во многом зависеть от формы ввода - если он следует какой-либо схеме, должны быть упрощения в наиболее общем случае. Я буду обновлять здесь, как я иду.
Edit:
У меня был разговор на этой неделе, который напомнил мне об этой проблеме. Если я перестану описывать гистограмму как уравнение и просто сохраню таблицу, могу ли я делать выборки за O (1) время? Оказывается, вы можете без потери точности за счет времени строительства O (N lgN).
Создать массив из N элементов. Равномерный случайный выбор в массиве найдет элемент с вероятностью 1 / N. Для каждого элемента сохраните долю попаданий, для которой этот элемент должен быть фактически выбран, и индекс другого элемента, который будет выбран, если этого нет.
Взвешенная случайная выборка, реализация C:
//data structure
typedef struct wrs_data {
double share;
int pair;
int idx;
} wrs_t;
//sort helper
int wrs_sharecmp(const void* a, const void* b) {
double delta = ((wrs_t*)a)->share - ((wrs_t*)b)->share;
return (delta<0) ? -1 : (delta>0);
}
//Initialize the data structure
wrs_t* wrs_create(int* weights, size_t N) {
wrs_t* data = malloc(sizeof(wrs_t));
double sum = 0;
int i;
for (i=0;i<N;i++) { sum+=weights[i]; }
for (i=0;i<N;i++) {
//what percent of the ideal distribution is in this bucket?
data[i].share = weights[i]/(sum/N);
data[i].pair = N;
data[i].idx = i;
}
//sort ascending by size
qsort(data,N, sizeof(wrs_t),wrs_sharecmp);
int j=N-1; //the biggest bucket
for (i=0;i<j;i++) {
int check = i;
double excess = 1.0 - data[check].share;
while (excess>0 && i<j) {
//If this bucket has less samples than a flat distribution,
//it will be hit more frequently than it should be.
//So send excess hits to a bucket which has too many samples.
data[check].pair=j;
// Account for the fact that the paired bucket will be hit more often,
data[j].share -= excess;
excess = 1.0 - data[j].share;
// If paired bucket now has excess hits, send to new largest bucket at j-1
if (excess >= 0) { check=j--;}
}
}
return data;
}
int wrs_pick(wrs_t* collection, size_t N)
//O(1) weighted random sampling (after preparing the collection).
//Randomly select a bucket, and a percentage.
//If the percentage is greater than that bucket's share of hits,
// use it's paired bucket.
{
int idx = rand_in_range(0,N);
double pct = rand_percent();
if (pct > collection[idx].share) { idx = collection[idx].pair; }
return collection[idx].idx;
}
Редактировать 2:
После небольшого исследования я обнаружил, что даже возможно сделать конструкцию за O (N) время. Благодаря тщательному отслеживанию вам не нужно сортировать массив, чтобы найти большие и маленькие ячейки. Обновленная реализация здесь