Schnellere Faltung von Wahrscheinlichkeitsdichtefunktionen in Python

Question

Angenommen, die Faltung einer allgemeinen Anzahl von diskreten Wahrscheinlichkeitsdichtefunktionen muss berechnet werden. Für das Beispiel unten gibt es vier Verteilungen, die auf Werten 0,1,2 mit den angegebenen Wahrscheinlichkeiten nehmen:

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]]) 

Die Faltung kann wie folgt gefunden werden:

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

Die Wahrscheinlichkeiten des Sehens 0, 1, ..., 8 sind dann gegeben durch

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

Dieser Teil der Engpass in meinem Code ist, und es scheint, dass es verfügbar etwas sein muss, um diesen Vorgang vektorisieren. Hat jemand einen Vorschlag, es schneller zu machen?

Alternativ kann eine Lösung, wo Sie

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) 

und erhalten die paarweise Faltungen

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

auch enorm helfen würde, nutzen könnten.

Answer 1

Sie können die Faltung all Ihrer PDFs effizient mit schnellen Fourier-Transformationen (FFTs) berechnen: Die wichtigste Tatsache ist, dass die FFT of the convolution das Produkt der FFTs der einzelnen Wahrscheinlichkeitsdichtefunktionen ist. Also transformiere jedes PDF, multipliziere die transformierten PDFs und führe dann die inverse Transformation durch. Sie müssen jede eingegebene PDF-Datei mit Nullen auf die entsprechende Länge auffüllen, um Effekte durch Umbrechen zu vermeiden.

Dies sollte einigermaßen effizient sein: wenn Sie m PDFs, die jeweils n Einträge, dann die Zeit die Faltung mit dieser Methode zu berechnen, sollte als (m^2)n log(mn) wachsen. Die Zeit wird von den FFTs dominiert, und wir berechnen effektiv m + 1 unabhängige FFTs (m Vorwärtstransformationen und eine inverse Transformation), jedes von einem Array mit einer Länge von nicht größer als mn. Aber wie immer, wenn Sie echte Timings wollen, sollten Sie profilieren.

Hier einige Code:

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 

Angewandt auf Ihrem Beispiel, hier ist was ich bekommen:

>>> 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]) 

, dass die grundlegende Idee. Wenn Sie dies optimieren möchten, können Sie auch numpy.fft.rfft (und seine inverse, numpy.fft.irfft) betrachten, die die Tatsache ausnutzen, dass die Eingabe real ist, um kompaktere transformierte Arrays zu erzeugen. Sie können möglicherweise auch etwas schneller werden, indem Sie das Array rows mit Nullen auffüllen, sodass die Gesamtzahl der Spalten optimal für die Durchführung einer FFT ist. Die Definition von "optimal" würde hier von der FFT-Implementierung abhängen, aber Zweierpotenzen wären zum Beispiel gute Ziele. Schließlich gibt es einige offensichtliche Vereinfachungen, die vorgenommen werden können, wenn rows erstellt wird, wenn alle Eingabearrays dieselbe Länge haben. Aber ich werde diese potenziellen Verbesserungen Ihnen überlassen.