Оптимизировать метод отклонения для генерации переменных

Question

У меня проблема с оптимизацией метода отклонения генерации непрерывных случайных величин.У меня плотность: f(x) = 3/2 (1-x^2).Вот мой код:

import random
import matplotlib.pyplot as plt
import numpy  as np
import time
import scipy.stats as ss

a=0   # xmin
b=1   # xmax

m=3/2 # ymax
variables = [] #list for variables

def f(x):
    return 3/2 * (1 - x**2)  #probability density function

reject = 0   # number of rejections
start = time.time()
while len(variables) < 100000:  #I want to generate 100 000 variables
    u1 = random.uniform(a,b)
    u2 = random.uniform(0,m)

    if u2 <= f(u1):
        variables.append(u1)
    else:
        reject +=1
end = time.time()

print("Time: ", end-start)
print("Rejection: ", reject)
x = np.linspace(a,b,1000)
plt.hist(variables,50, density=1)
plt.plot(x, f(x))
plt.show()

ss.probplot(variables, plot=plt)
plt.show()

Мой первый вопрос: правильно ли построен мой график вероятности?И второе, что есть в названии.Как оптимизировать этот метод?Я хотел бы получить несколько советов по оптимизации кода.Теперь этот код занимает около 0,5 секунд, и есть около 50 000 отклонений.Можно ли сократить время и количество отказов?При необходимости я могу оптимизировать, используя другой метод генерации переменных.

Answer 1

Что касается вашего первого вопроса, scipy.stats.probplot сравнивает вашу выборку с квантилями нормального распределения. Если вы хотите сравнить его с квантилями вашего f(x) дистрибутива, проверьте параметр dist в probplot.

С точки зрения ускорения этой процедуры отбора проб, как правило, следует избегать циклов. Замена кода между start = ... и end = ... на следующее привела к ускорению> 20x.

n_before_accept_reject = 150000
u1 = np.random.uniform(a, b, size=n_before_accept_reject)
u2 = np.random.uniform(0, m, size=n_before_accept_reject)
variables = u1[u2 <= f(u1)]
reject = n_before_accept_reject - len(variables)

Обратите внимание, что при каждом запуске вы будете получать приблизительно 100000 принятых образцов. Вы можете немного увеличить значение n_before_accept_reject, чтобы эффективно гарантировать, что variables всегда будет иметь> 100000 принятых значений, а затем просто ограничить размер переменных, чтобы при необходимости получить ровно 100000.

Answer 2

Другие говорили о вероятностном построении, я собираюсь обсудить эффективность алгоритма отклонения.

Схемы принятия / отклонения основаны на m (x), «мажорирующей функции».Мажорирующая функция должна иметь два свойства: 1) m (x) ≥f (x) ∀ x;и 2) m (x), когда масштабируется как распределение, должно легко генерировать значения из.Вы использовали постоянную функцию m = 3/2, которая удовлетворяет обоим требованиям, но не очень тесно связывает f (x).Интегрируется от нуля до единицы, имеет площадь 3/2.Ваша функция f (x), будучи действительной функцией плотности, имеет площадь 1. Следовательно, ∫f (x)) / ∫m (x)) = 1 / (3/2) = 2/3.Другими словами, принимаются 2/3 значений, которые вы генерируете из мажорирующей функции, и вы отклоняете 1/3 попыток.

Вам нужен m (x), который обеспечивает более жесткую границу для f(Икс).Я пошел с линией, которая касается f (x) в x = 1/2.С небольшим исчислением, чтобы получить наклон, я получил m(x) = 15/8 - 3x/2.

Этот выбор m (x) имеет площадь 9/ 8, поэтому только 1/9 значений будет отклонено.Чуть больше исчисления дало генератор обратного преобразования для x, основанный на этом m (x): x = (5 - sqrt(25 - 24U)) / 4, где U - это равномерная (0,1) случайная переменная.

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

import random
import matplotlib.pyplot as plt
import numpy  as np
import time
import math
import scipy.stats as ss

a = 0   # xmin
b = 1   # xmax

reject = 0   # number of rejections

def f(x):
    return 3.0 / 2.0 * (1.0 - x**2)  #probability density function

def m(x):
    return 1.875 - 1.5 * x

def generate_x():
    global reject
    while True:
        x = (5.0 - math.sqrt(25.0 - random.uniform(0.0, 24.0))) / 4.0
        u = random.uniform(0, m(x))
        if u <= f(x):
            return x 
        reject += 1    

start = time.time()
variables = [generate_x() for _ in range(100000)]
end = time.time()

print("Time: ", end-start)
print("Rejection: ", reject)
x = np.linspace(a,b,1000)
plt.hist(variables,50, density=1)
plt.plot(x, f(x))
plt.show()

Answer 3

Мой первый вопрос: правильно ли построен мой график вероятности?

Нет.Это сделано по сравнению с нормальным распределением по умолчанию.Вы должны упаковать свою функцию f(x) в класс, полученный из stats.rv_continuous, превратить ее в метод _pdf и передать ее в probplot

И второе, что находится в заголовке.Как оптимизировать этот метод?Можно ли сократить время и количество отклонений?

Конечно, у вас есть сила векторных способностей NumPy в ваших руках.Никогда не пишите явные циклы - векторизация, векторизация и векторизация!

Посмотрите на приведенный ниже модифицированный код, а не на один цикл, все делается через векторы NumPy.Время на моем компьютере сократилось на 100000 образцов (Xeon, Win10 x64, Anaconda Python 3.7) с 0,19 до 0,003.

import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import time

a = 0.  # xmin
b = 1.  # xmax

m = 3.0/2.0 # ymax

def f(x):
    return 1.5 * (1.0 - x*x)  # probability density function

start  = time.time()

N = 100000
u1 = np.random.uniform(a, b, N)
u2 = np.random.uniform(0.0, m, N)

negs = np.empty(N)
negs.fill(-1)
variables = np.where(u2 <= f(u1), u1, negs) # accepted samples are positive or 0, rejected are -1

end = time.time()

accept = np.extract(variables>=0.0, variables)
reject = N - len(accept)

print("Time: ", end-start)
print("Rejection: ", reject)

x = np.linspace(a, b, 1000)
plt.hist(accept, 50, density=True)
plt.plot(x, f(x))
plt.show()

ss.probplot(accept, plot=plt) # against normal distribution
plt.show()

Что касается уменьшения количества отклонений, вы можете выполнить выборку с 0 отклонениями, используя метод обратного анализа,является кубическим уравнением, поэтому он может работать с легкими

ОБНОВЛЕНИЕ

Вот код, который нужно использовать для probplot:

class my_pdf(ss.rv_continuous):
    def _pdf(self, x):
        return 1.5 * (1.0 - x*x)

ss.probplot(accept, dist=my_pdf(a=a, b=b, name='my_pdf'), plot=plt)

, и вы должны получить что-то вроде