Szybsze przekształcanie funkcji gęstości prawdopodobieństwa w Pythonie

Question

Załóżmy, że splot ogólnej liczby dyskretnych funkcji gęstości prawdopodobieństwa musi zostać obliczony. Dla przykładu poniżej istnieją cztery dystrybucje, które odbywają się na wartościach 0,1,2 z określonymi prawdopodobieństwami:

import numpy as np 
pdfs = np.array([[0.6,0.3,0.1],[0.5,0.4,0.1],[0.3,0.7,0.0],[1.0,0.0,0.0]]) 

splotu można znaleźć tak:

pdf = pdfs[0]   
for i in range(1,pdfs.shape[0]): 
    pdf = np.convolve(pdfs[i], pdf) 

Prawdopodobieństwa widząc 0, 1, ..., 8 są następnie przez

array([ 0.09 , 0.327, 0.342, 0.182, 0.052, 0.007, 0. , 0. , 0. ]) 

Ta część jest wąskim gardłem w moim kodu i wydaje się, że coś musi być dostępny do wektoryzacji tę operację. Czy ktoś ma sugestię, aby przyspieszyć działanie?

Alternatywnie, roztwór, gdzie można użyć

pdf1 = np.array([[0.6,0.3,0.1],[0.5,0.4,0.1]]) 
pdf2 = np.array([[0.3,0.7,0.0],[1.0,0.0,0.0]]) 
convolve(pd1,pd2) 

i uzyskać parami zwoje

array([[ 0.18, 0.51, 0.24, 0.07, 0. ], 
     [ 0.5, 0.4, 0.1, 0. , 0. ]]) 

pomogłoby także ogromnie.

Answer 1

Można efektywnie obliczyć splot wszystkich plików PDF za pomocą szybkich transformatów Fouriera (FFT): kluczowym faktem jest to, że FFT of the convolution jest produktem funkcji FFT poszczególnych funkcji gęstości prawdopodobieństwa. Więc przekształć każdy plik PDF, pomnóż przekształcone pliki PDF razem, a następnie wykonaj przekształcenie odwrotne. Będziesz musiał dopełnić każdy wejściowy plik PDF zerami do odpowiedniej długości, aby uniknąć efektów zawijania.

ten powinien być w miarę wydajny: jeśli masz m PDF, każda zawierająca n wpisy, to czas, aby obliczyć splot przy użyciu tej metody powinny rosnąć jak (m^2)n log(mn). Czas jest zdominowany przez FFT, a my efektywnie obliczamy niezależne transformaty FFT (m transformaty do przodu i jedną transformatę odwrotną), każda z tablic o długości nie większej niż mn. Ale jak zawsze, jeśli chcesz mieć prawdziwe czasy, powinieneś profilować.

Oto niektóre kodu:

import numpy.fft 

def convolve_many(arrays): 
    """ 
    Convolve a list of 1d float arrays together, using FFTs. 
    The arrays need not have the same length, but each array should 
    have length at least 1. 

    """ 
    result_length = 1 + sum((len(array) - 1) for array in arrays) 

    # Copy each array into a 2d array of the appropriate shape. 
    rows = numpy.zeros((len(arrays), result_length)) 
    for i, array in enumerate(arrays): 
     rows[i, :len(array)] = array 

    # Transform, take the product, and do the inverse transform 
    # to get the convolution. 
    fft_of_rows = numpy.fft.fft(rows) 
    fft_of_convolution = fft_of_rows.prod(axis=0) 
    convolution = numpy.fft.ifft(fft_of_convolution) 

    # Assuming real inputs, the imaginary part of the output can 
    # be ignored. 
    return convolution.real 

Stosując to do Twojego przykładu, oto co mam:

>>> convolve_many([[0.6, 0.3, 0.1], [0.5, 0.4, 0.1], [0.3, 0.7], [1.0]]) 
array([ 0.09 , 0.327, 0.342, 0.182, 0.052, 0.007]) 

to podstawowa idea. Jeśli chcesz to poprawić, możesz również spojrzeć na numpy.fft.rfft (i jego odwrotność, numpy.fft.irfft), które wykorzystują fakt, że dane wejściowe są prawdziwe, aby utworzyć bardziej zwarte transformowane tablice. Możesz również uzyskać pewną prędkość, wypełniając tablicę rows zerami, aby całkowita liczba kolumn była optymalna do wykonywania FFT. Definicja "optymalnej" tutaj zależałaby od implementacji FFT, ale uprawnienia dwóch mogłyby być na przykład dobrym celem. Wreszcie, istnieją pewne oczywiste uproszczenia, które można wprowadzić podczas tworzenia rows, jeśli wszystkie tablice wejściowe mają tę samą długość. Ale zostawię ci te potencjalne ulepszenia.