Allowing convolutions of differences as well as additions

pythonScripts/convolveLN.py

@@ -65,6 +65,17 @@ def getLTransformGrid(sigma,fz, mu=0):
     fn=len(fa)
     return numpy.array([getLTransform(z,fa,fb,mu) for z in fz])
 
+def getComplexConjugatedLTransformAtMinusComplexConjugatedZGrid(sigma,fz, mu=0):
+   #evaluate the Laplace transform of a gamma-set convolution 
+   #that approximates a log-normal distribution with parameters
+   #sigma_s=sigma, mu_s=mu at an array of complex numbers fz
+    fa,fb=readCoef(sigma)
+    fn=len(fa)
+    zPrime=-numpy.conj(fz)
+    return numpy.conj(numpy.array([getLTransform(z,fa,fb,mu) for z in zPrime]))
+
+
+
 def getLTransformPart(z, fa, fb, i):
    #get a Laplace transform of a single gamma distribution, specified by index
    #into a set of alpha and beta arrays fa and fb, respectively
@@ -83,13 +94,11 @@ def addExpA(fz,fq,x):
     #apply the exp(c*x)*exp(Re(a)*u*x term, see paper for reasoning
     return fq*numpy.exp(numpy.real(fz)*x)
 
-def convolveLN(fx,sigma,mu):
-   #convolve a set of LN with parameters sigma and mu as numpy arrays and
-   #return the pdf of convolution evaluated at the set of points fx
-    
+
+def getZArray(n=2000):
     #internal computation parameters
     #number of grid points 1M suggested in paper, good fit found for few thousand
-    n=2000
+    #n=2000
     #integral range, 250 to 500 suggested in paper
     u0=0
     u1=300
@@ -104,23 +113,38 @@ def convolveLN(fx,sigma,mu):
     h=(u1-u0)/2/n
     #curve segmentation
     qz=c+fa/numpy.absolute(fa)*qh
-    #function values at grid points
-    qPhi=numpy.ones(qh.size,dtype=complex)
-    for i in range(len(sigma)):
-        m=mu[i]
-        s=sigma[i]
-        qPhi*=getLTransformGrid(sigma[i],qz,mu[i])
-        
+    return qz,fa,u0,h
+
+def inverseL(fx,qz,qPhi,fa,u0,h,n):
+   #Filion quadratures by parts of complex numbers
+   fI=numpy.zeros(len(fx),dtype=complex)
+   for i in range(len(fx)):
+      x=fx[i]
+      qPhi1=addExpA(qz,qPhi,x)
+      fRc=filionQuadrature.filionQuadrature(numpy.real(qPhi1),x,u0,h,n,'cos')
+      fIc=filionQuadrature.filionQuadrature(numpy.imag(qPhi1),x,u0,h,n,'cos')
+      fRs=filionQuadrature.filionQuadrature(numpy.real(qPhi1),x,u0,h,n,'sin')
+      fIs=filionQuadrature.filionQuadrature(numpy.imag(qPhi1),x,u0,h,n,'sin')
+      fI[i]=numpy.complex(fRc-fIs,fIc+fRs)
+   fI1=fa*fI
+   return numpy.imag(fI1)
+
+def convolveLN(fx,sigma,mu,sign):
+   #convolve a set of LN with parameters sigma and mu as numpy arrays and
+   #return the pdf of convolution evaluated at the set of points fx
     
-    #Filion quadratures by parts of complex numbers
-    fI=numpy.zeros(len(fx),dtype=complex)
-    for i in range(len(fx)):
-        x=fx[i]
-        qPhi1=addExpA(qz,qPhi,x)
-        fRc=filionQuadrature.filionQuadrature(numpy.real(qPhi1),x,u0,h,n,'cos')
-        fIc=filionQuadrature.filionQuadrature(numpy.imag(qPhi1),x,u0,h,n,'cos')
-        fRs=filionQuadrature.filionQuadrature(numpy.real(qPhi1),x,u0,h,n,'sin')
-        fIs=filionQuadrature.filionQuadrature(numpy.imag(qPhi1),x,u0,h,n,'sin')
-        fI[i]=numpy.complex(fRc-fIs,fIc+fRs)
-    fI1=fa*fI
-    return numpy.imag(fI1)
+   #zArray
+   #number of grid points 1M suggested in paper, good fit found for few thousand
+   n=2000
+   qz,fa,u0,h=getZArray(n)
+   #function values at grid points
+   qPhi=numpy.ones(qz.size,dtype=complex)
+   for i in range(len(sigma)):
+      m=mu[i]
+      s=sigma[i]
+      if sign[i]:
+         qPhi*=getLTransformGrid(sigma[i],qz,mu[i])
+         continue
+      qPhi*=getComplexConjugatedLTransformAtMinusComplexConjugatedZGrid(s,qz,m)
+        
+   return inverseL(fx,qz,qPhi,fa,u0,h,n)