Changing propagateErrorLN to separatly convolve positive and negative derivatives and perform cross-correlation at the end

pythonScripts/convolveLN.py

@@ -41,9 +41,22 @@ def readCoef(sigma):
     fbeta=[]
     for i in range(fn):
         f2a = scipy.interpolate.interp1d(fx, fa[:,i], kind='cubic')
-        falpha.append(f2a(sigma))
+        try:
+         falpha.append(f2a(sigma))
+        except ValueError:
+         #interpolation failed
+         if sigma<fx[0]:
+            falpha.append(fa[0,i])
+            print('Using alpha[{}]={} at sigma={}'.format(fx[0],fa[0,i],sigma))
+
         f2b = scipy.interpolate.interp1d(fx, fb[:,i], kind='cubic')
-        fbeta.append(f2b(sigma))
+        try:
+         fbeta.append(f2b(sigma))
+        except ValueError:
+         if sigma<fx[0]:
+            fbeta.append(fb[0,i])
+            print('Using beta[{}]={} at sigma={}'.format(fx[0],fb[0,i],sigma))
+
     return numpy.array(falpha),numpy.array(fbeta)
 
 def getLTransform(z,fa,fb, mu=0):
@@ -61,6 +74,7 @@ def getLTransformGrid(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
+    print('Reading for sigma={}'.format(sigma))
     fa,fb=readCoef(sigma)
     fn=len(fa)
     return numpy.array([getLTransform(z,fa,fb,mu) for z in fz])
@@ -69,6 +83,7 @@ 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
+    print('Reading for sigma={}'.format(sigma))
     fa,fb=readCoef(sigma)
     fn=len(fa)
     zPrime=-numpy.conj(fz)
@@ -129,7 +144,7 @@ def inverseL(fx,qz,qPhi,fa,u0,h,n):
    fI1=fa*fI
    return numpy.imag(fI1)
 
-def convolveLN(fx,sigma,mu,sign):
+def convolveLN(fx,sigma,mu,sign=numpy.array([])):
    #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
     
@@ -142,9 +157,15 @@ def convolveLN(fx,sigma,mu,sign):
    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)
+      try:
+         if not sign[i]:
+            qPhi*=getComplexConjugatedLTransformAtMinusComplexConjugatedZGrid(s,qz,m)
+            continue
+      except IndexError:
+         pass
+      
+      qPhi*=getLTransformGrid(sigma[i],qz,mu[i])
         
-   return inverseL(fx,qz,qPhi,fa,u0,h,n)    
+   return inverseL(fx,qz,qPhi,fa,u0,h,n)   
+
+   

pythonScripts/propagateErrorLN.py

@@ -1,5 +1,8 @@
 import numpy
+import scipy.interpolate
+import importlib
 import convolveLN
+importlib.reload(convolveLN)
 
 def generateSample(mu,cv):
    return numpy.array([convolveLN.ln(mu,cv) for i in range(10000)])
@@ -23,17 +26,90 @@ def calculateDistribution(fy,muTarget,mu,cv,dydp):
    #muTarget is the target mean of y distribution
    #fy is the set of values y where distribution is evaluated at
    #Parameters mu, cv and dydp are given as vectors/lists
-    
+   
+
+   #parameters qa should be around 1
    qa=mu*dydp
-   A=numpy.sum(qa)
-   qa*=muTarget/A
-   cvPrime=cv*A/muTarget
-   sigmaS=numpy.sqrt(numpy.log(1+cvPrime*cvPrime))
+   print('qa={}'.format(qa))
    sign=qa>=0
+   Ax=numpy.sum(qa[sign])
+   Ay=numpy.sum(-qa[~sign])
+   print('A={},{}'.format(Ax,Ay))
+   if Ax>0:
+      qa[sign]/=Ax
+   if Ay>0:
+      qa[~sign]/=Ay
+   cvPrime=cv
+   print('cvPrime={}'.format(cvPrime))
+
+   #convolve qa multiplied (1,cv) distributions 
+   sigmaS=numpy.sqrt(numpy.log(1+cvPrime*cvPrime))
+   print('sigmaS={}'.format(sigmaS))
    muS=numpy.zeros(qa.shape)
    muS[sign]=numpy.log(qa[sign]/numpy.sqrt(1+cvPrime[sign]*cvPrime[sign]))
    muS[~sign]=numpy.log(-qa[~sign]/numpy.sqrt(1+cvPrime[~sign]*cvPrime[~sign]))
-   y=convolveLN.convolveLN(fy,sigmaS,muS,sign)
+   print('muS={}'.format(muS))
+
+   #do the convolution
+   nbins=100
+   umax=5
+   fw=numpy.linspace(0,umax,nbins)
+   wx=numpy.array([])
+   wy=numpy.array([])
+   
+   if Ax>0:
+      wx=convolveLN.convolveLN(fw,sigmaS[sign],muS[sign])
+   
+   if Ay>0:
+      wy=convolveLN.convolveLN(fw,sigmaS[~sign],muS[~sign])
+
+   #for correlate, binning of wx and wy should be equal, but one has scale A1, the other A2
+   #how do you deal with that?
+   if Ax>0:
+      fwx=Ax*fw
+   if Ay>0:
+      fwy=Ay*fw
+
+   fw=fwx
+   if not Ay>Ax:
+      if wy.size>0:
+      #re-sample wy
+         wy = scipy.interpolate.CubicSpline(fwy, wy)(fw)
+
+   if Ay>Ax:
+      fw=fwy
+      #resample wx
+      if wx.size>0:
+         wx = scipy.interpolate.CubicSpline(fwx, wx)(fw)
+   
+   #output is twice the length of w and symmetric around 0
+   fw1=numpy.concatenate((numpy.flip(-fw[1:]),fw))
+   w=numpy.zeros(fw1.size)
+
+   if wx.size==0:
+      w[:len(wy)]=numpy.flip(wy)
+   if wy.size==0:
+      w[-len(wx):]=wx
+   if wx.size>0 and wy.size>0:
+      w=scipy.signal.correlate(wx,wy)
+      #w is a convolution of input LN distributions with a mean of A
+
+
+   #now shift it to the target meat muTarget
+   B=muTarget-(Ax-Ay)
+   print('B={}'.format(B))
+
+   #interpolate on w and fill y
+   f = scipy.interpolate.CubicSpline(fw1, w)
+     
+   y=numpy.zeros(fy.shape)
+   
+   for i in range(len(y)):
+      try:
+         y[i]=f(fy[i]-B)
+      except ValueError:
+         pass
+   
    return y/numpy.sum(y)
 
 def generateDistribution(fx,muTarget,mu,cv,dydp, dbg=0):
@@ -47,22 +123,23 @@ def generateDistribution(fx,muTarget,mu,cv,dydp, dbg=0):
    
    qa=mu*dydp
    A=numpy.sum(qa)
-   qa*=muTarget/A
-   cvPrime=cv*A/muTarget
+   #qa*=A
+   print('qa={}'.format(qa))
+   cvPrime=cv
    samples=[a*generateSample(1,c) for a,c in zip(qa,cvPrime)]
    mean=[numpy.mean(s) for s in samples]
    std=[numpy.std(s) for s in samples]
    convSample=numpy.zeros(len(samples[0]))
-   msg='[{}] mean={:.2f} std={:.2f} cv={:.2f}'
+   msg='[{}] mean={:.2g} std={:.2g} cv={:.2f}'
    for i in range(len(cv)):
       convSample+=samples[i]
       if dbg:
          print(msg.format(i,mean[i],std[i],std[i]/mean[i]))
-      
+   convSample+=muTarget-A   
    convM=numpy.mean(convSample)
    convS=numpy.std(convSample)
    if dbg:
-      code=x
+      code="x"
       print(msg.format(code,convM,convS,convS/convM))
    h,be=numpy.histogram(convSample,bins=getEdges(fx))
    #h should be len(fx)