PBPK_public/pythonScripts/propagateErrorLN.py

import numpy
import convolveLN

def generateSample(mu,cv):
   return numpy.array([convolveLN.ln(mu,cv) for i in range(10000)])

def getEdges(a):
   #convert array of bin centers a to array of bin edges with one more element than a
   a1=numpy.zeros(len(a)+1)
   a2=numpy.zeros(a1.shape)
   a1[1:]=a
   a2[:-1]=a
   e=0.5*(a1+a2)
   e[0]=2*e[1]-e[2]
   e[-1]=2*e[-2]-e[-3]
   return e

def calculateDistribution(fy,muTarget,mu,cv,dydp):
   #calculate distribution of random variable y with a functional dependence 
   #on random parameters p, each with a log-normal 
   #distribution with a mean mu(p) and coefficients of variation cv(p) 
   #and dydp giving partial derivatives of y on  p
   #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
    
   qa=mu*dydp
   A=numpy.sum(qa)
   qa*=muTarget/A
   cvPrime=cv*A/muTarget
   sigmaS=numpy.sqrt(numpy.log(1+cvPrime*cvPrime))
   sign=qa>=0
   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)
   return y/numpy.sum(y)

def generateDistribution(fx,muTarget,mu,cv,dydp, dbg=0):
   #sum randomly generated variables to approximate distribution of 
   #variable y with a functional dependence y(p) on random parameters p, 
   #each with a log-normal (LN) distribution 
   #with a mean mu(p) and coefficient of variation cv(p)
   #target mean of y is muTarget
   #distribution is evaluated at points fy, 
   #dydp are derivatives of y with respect to parameters p
   
   qa=mu*dydp
   A=numpy.sum(qa)
   qa*=muTarget/A
   cvPrime=cv*A/muTarget
   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}'
   for i in range(len(cv)):
      convSample+=samples[i]
      if dbg:
         print(msg.format(i,mean[i],std[i],std[i]/mean[i]))
      
   convM=numpy.mean(convSample)
   convS=numpy.std(convSample)
   if dbg:
      code=x
      print(msg.format(code,convM,convS,convS/convM))
   h,be=numpy.histogram(convSample,bins=getEdges(fx))
   #h should be len(fx)
   return h/numpy.sum(h)