Сетка над векторами вероятности - PullRequest
Новая
       5

Сетка над векторами вероятности

0 голосов
/ 29 мая 2018

Я пытаюсь получить «сетку» из n-мерных вероятностных векторов - векторов, в которых каждая запись находится в диапазоне от 0 до 1, и все записи составляют в сумме 1. Я хочу, чтобы все возможные векторы, в которых координатыможет принимать любое из числа v равномерно распределенных значений между 0 и 1.

Чтобы проиллюстрировать это, то, что следует, является ужасно неэффективной реализацией, для n = 3 и v = 3: 

from itertools import product
grid_redundant = product([0, .5, 1], repeat=3)
grid = [point for point in grid_redundant if sum(point)==1]

сейчас grid содержит [(0, 0, 1), (0, 0.5, 0.5), (0, 1, 0), (0.5, 0, 0.5), (0.5, 0.5, 0), (1, 0, 0)].

Эта "реализация" просто ужасна для сеток с большими размерами и более мелкозернистыми.Есть ли хороший способ сделать это, возможно, используя numpy?

Я мог бы добавить момент о мотивации: я был бы совершенно счастлив, если бы простая выборка из случайного распределения дала мне достаточно экстремальные очки, но это не так.См. этот вопрос .«Решетка», которую я ищу, не случайна, а систематически охватывает симплекс (пространство векторов вероятности).

Ответы [ 2 ]

0 голосов
/ 05 июня 2018

Делать это в полной общности для больших размерных векторов и даже с умным решением в принятом ответе довольно неуправляемо.В моем собственном случае стоит рассчитать соответствующее подмножество всех значений.Например, следующая функция вычисляет все dimension -мерные вероятностные векторы только с n ненулевыми равновероятными записями: 

import itertools as it
import numpy as np

def equip_n(dimension, n):
"""
Calculate all possible <dimension>-dimensional probability vectors with n nonzero,
equiprobable entries
"""
combinations  = np.array([comb for comb in it.combinations(range(dimension), n)])
vectors = np.zeros((combinations.shape[0], dimension))
for line, comb in zip(vectors, combinations):
    line[comb] = 1/n
return vectors 

print(equip_n(6, 3))

Возвращает 

[[ 0.3333  0.3333  0.3333  0.      0.      0.    ]
 [ 0.3333  0.3333  0.      0.3333  0.      0.    ] 
 [ 0.3333  0.3333  0.      0.      0.3333  0.    ]
 [ 0.3333  0.3333  0.      0.      0.      0.3333]
 [ 0.3333  0.      0.3333  0.3333  0.      0.    ]
 [ 0.3333  0.      0.3333  0.      0.3333  0.    ]
 [ 0.3333  0.      0.3333  0.      0.      0.3333]
 [ 0.3333  0.      0.      0.3333  0.3333  0.    ]
 [ 0.3333  0.      0.      0.3333  0.      0.3333]
 [ 0.3333  0.      0.      0.      0.3333  0.3333]
 [ 0.      0.3333  0.3333  0.3333  0.      0.    ]
 [ 0.      0.3333  0.3333  0.      0.3333  0.    ]
 [ 0.      0.3333  0.3333  0.      0.      0.3333]
 [ 0.      0.3333  0.      0.3333  0.3333  0.    ]
 [ 0.      0.3333  0.      0.3333  0.      0.3333]
 [ 0.      0.3333  0.      0.      0.3333  0.3333]
 [ 0.      0.      0.3333  0.3333  0.3333  0.    ]
 [ 0.      0.      0.3333  0.3333  0.      0.3333]
 [ 0.      0.      0.3333  0.      0.3333  0.3333]
 [ 0.      0.      0.      0.3333  0.3333  0.3333]]

Это очень быстро.%timeit equip_n(6, 3) возвращает 

15.1 µs ± 74.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
0 голосов
/ 29 мая 2018

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

import math
from itertools import permutations

def probability_grid(values, n):
    values = set(values)
    # Check if we can extend the probability distribution with zeros
    with_zero = 0. in values
    values.discard(0.)
    if not values:
        raise StopIteration
    values = list(values)
    for p in _probability_grid_rec(values, n, [], 0.):
        if with_zero:
            # Add necessary zeros
            p += (0.,) * (n - len(p))
        if len(p) == n:
            yield from set(permutations(p))  # faster: more_itertools.distinct_permutations(p)

def _probability_grid_rec(values, n, current, current_sum, eps=1e-10):
    if not values or n <= 0:
        if abs(current_sum - 1.) <= eps:
            yield tuple(current)
    else:
        value, *values = values
        inv = 1. / value
        # Skip this value
        yield from _probability_grid_rec(
            values, n, current, current_sum, eps)
        # Add copies of this value
        precision = round(-math.log10(eps))
        adds = int(round((1. - current_sum) / value, precision))
        for i in range(adds):
            current.append(value)
            current_sum += value
            n -= 1
            yield from _probability_grid_rec(
                values, n, current, current_sum, eps)
        # Remove copies of this value
        if adds > 0:
            del current[-adds:]

print(list(probability_grid([0, 0.5, 1.], 3)))

Вывод: 

[(1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0), (0.5, 0.5, 0.0), (0.0, 0.5, 0.5), (0.5, 0.0, 0.5)]

Быстрое сравнение с опубликованным методом: 

from itertools import product

def probability_grid_basic(values, n):
    grid_redundant = product(values, repeat=n)
    return [point for point in grid_redundant if sum(point)==1]

values = [0, 0.25, 1./3., .5, 1]
n = 6
%timeit list(probability_grid(values, n))
1.61 ms ± 20.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit probability_grid_basic(values, n)
6.27 ms ± 186 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
...