Как ускорить генерацию случайных массивов в python - PullRequest
0 голосов
/ 10 июня 2018

Я хочу использовать метод интеграции Монте-Карло, чтобы найти объем n-сферы и сравнить его с аналитическим результатом.Я хочу сделать это для 10 ^ 6 баллов и 10 ^ 9 баллов, и хотя это работает для 10 ^ 6 баллов (занимает около минуты для n = 2 (круг), n = 3 (сфера) и n = 12),это ужасно медленно для 10 ^ 9 баллов.

Краткое объяснение MC-метода: чтобы найти Объем n-сферы с радиусом r = 1, я представляю себе легко известный Объем (скажем, n-куб со сторонами длины 2 * r), который полностьюсодержит н-сферу.Затем я выбираю из равномерного распределения точек в n-кубе и проверяю, лежит ли точка в сфере.Я считаю все сгенерированные точки внутри n-сферы.Отношение V_sphere / V_cube может быть аппроксимировано как N_inside / N_total и, следовательно, V_sphere = V_cube * N_inside / N_total

Вот функция:

def hyp_sphere_mc(d,samples):    

    inside = 0                     #number of points inside sphere                                             
    sum = 0                        #sum of squared components

    for j in range(0,samples):        

        x2 = np.random.uniform(0,1,d)**2     #generate random point in d-dimensions                           
        sum = np.sum(x2)                     #sum its components

        if(sum < 1.0):                                              
            inside += 1                      #count points inside sphere

    V = ((2)**d)*(float(inside)/samples)     #V = V_cube * N_inside/N_total           

    V_true = float(math.pi**(float(d)/2))/math.gamma(float(d)/2 + 1) #analytical result 

    ERR = (float(abs(V_true-V))/V_true)*100        #relative Error

    print "Numerical:", V, "\t" , "Exact: ", V_true, "\t", "Error: ", ERR

Я думаю, проблема в том, чтодля каждой итерации я генерирую новый случайный массив, и это занимает много времени, особенно если у меня есть 10 ^ 9 итераций.Есть ли способ ускорить это?

1 Ответ

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

Вы можете заменить цикл следующим:

inside = np.sum(np.sum(np.random.rand(samples,d)**2,1)<1)

При использовании numpy вы должны стараться избегать циклов.Идея состоит в том, что вы можете просто сгенерировать все образцы одновременно в матрице, а затем векторизовать все последующие операции.

...