I'm trying to implement some inventory models in python. We're using the Gamma distribution from scipy.stats for the lead time demand distribution for each part. The service level metrics for these models involve certain loss functions, which need to be added to the rv_frozen object representing the lead time demand distribution. The loss functions for the Gamma distribution are:
def first_order_loss_fn(self, x):
"""F^1(x)=E{[X-x]^+} Used in many inventory calculations. Implemented from Zipkin Appendix C (p. 457)"""
term1 = (self.a - self.scale * x) * self.sf(x)
term2 = x * self.pdf(x)
return (term1 + term2) / self.scale
def second_order_loss_fn(self, x):
"""F^2(x)=1/2*E{[X-x]^+[X-x-1]^+. Using the form given by Zipkin (2000) Appendix C (p. 457)"""
abx = self.a - self.scale * x
term1 = self.sf(x) * (abx ** 2 - self.a)
term2 = self.pdf(x) * (abx - 1)
return 0.5 * (term1 + term2) / (self.scale ** 2)
We may eventually use other distributions from scipy.stats, each of which has their own equations for the two loss functions. Subclassing rv_frozen for each distribution wouldn't be an issue on its own, but it seems like you'd have to mess with the dist_gen classes (eg gamma_gen, geom_gen, etc.) as well.
What is the best way to add these two loss functions to a frozen scipy.stats distribution -- initially Gamma, but geometric and negative binomial are both on the horizon?