为scipy.stats.rv_frozen添加分布专属损失函数的最优方案咨询

Original Problem #

I'm implementing inventory models in Python, using scipy.stats distributions (starting with Gamma) to represent lead-time demand. These models require specific loss functions for service level calculations, which I need to attach to the frozen distribution objects (rv_frozen).

For the Gamma distribution, I've already written the loss functions:

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)

I'll eventually use other distributions like Geometric and Negative Binomial, each with their own loss function formulas. Subclassing rv_frozen for each distribution works, but it seems I'd also need to modify the generator classes (e.g., gamma_gen, geom_gen). What's the optimal way to add these two loss functions to frozen SciPy distributions, both for Gamma now and other distributions later?

Recommended Solution: Wrapper Class (Composition Over Inheritance) #

The cleanest, most scalable approach is to create a wrapper class that wraps an rv_frozen object and adds your custom loss functions. This avoids modifying SciPy's internal code, works across all distributions, and keeps your loss logic organized.

Step 1: Base Wrapper + Distribution-Specific Subclasses #

First, create a base class that handles common delegation, then subclass it for each distribution to implement their unique loss functions:

from scipy.stats import gamma, geom, nbinom

class InventoryDistribution:
    def __init__(self, frozen_dist):
        self.dist = frozen_dist  # Wrap the rv_frozen object
        # Delegate all existing SciPy methods/properties to the wrapped distribution
        self.__dict__.update(frozen_dist.__dict__)
    
    def first_order_loss_fn(self, x):
        """Base method - override in subclasses for each distribution"""
        raise NotImplementedError("First-order loss function not implemented for this distribution")
    
    def second_order_loss_fn(self, x):
        """Base method - override in subclasses for each distribution"""
        raise NotImplementedError("Second-order loss function not implemented for this distribution")

# Gamma Distribution Subclass
class GammaInventoryDistribution(InventoryDistribution):
    def first_order_loss_fn(self, x):
        term1 = (self.dist.a - self.dist.scale * x) * self.dist.sf(x)
        term2 = x * self.dist.pdf(x)
        return (term1 + term2) / self.dist.scale
    
    def second_order_loss_fn(self, x):
        abx = self.dist.a - self.dist.scale * x
        term1 = self.dist.sf(x) * (abx ** 2 - self.dist.a)
        term2 = self.dist.pdf(x) * (abx - 1)
        return 0.5 * (term1 + term2) / (self.dist.scale ** 2)

# Geometric Distribution Subclass (example implementation)
class GeometricInventoryDistribution(InventoryDistribution):
    def first_order_loss_fn(self, x):
        p = self.dist.p
        if x < 0:
            return (1 - p)/p
        return ((1 - p)**(x+1)) / p
    
    def second_order_loss_fn(self, x):
        p = self.dist.p
        if x < 0:
            return (1 - p)*(2 - p)/(p**2)
        return ((1 - p)**(x+1) * (x + 1 - p)) / (p**2)

# Negative Binomial Subclass (example implementation)
class NegativeBinomialInventoryDistribution(InventoryDistribution):
    def first_order_loss_fn(self, x):
        n = self.dist.n
        p = self.dist.p
        return ((n * (1 - p))/p) * self.dist.sf(x) - x * self.dist.sf(x) + self.dist.expect(lambda y: y, lb=x+1)
    
    def second_order_loss_fn(self, x):
        # Implement your derived Zipkin formula here
        return 0.5 * self.dist.expect(lambda y: (y - x)*(y - x - 1), lb=x+1)

Step 2: Usage Example #

Create your frozen distribution as usual, then wrap it with your custom inventory class to access both standard SciPy methods and your loss functions:

# Create a frozen Gamma distribution
gamma_frozen = gamma(a=2.5, scale=3.0)
# Wrap it with your custom class
gamma_inv_dist = GammaInventoryDistribution(gamma_frozen)

# Use standard SciPy methods alongside your custom loss functions
print(gamma_inv_dist.mean())  # Standard SciPy mean calculation
print(gamma_inv_dist.first_order_loss_fn(10))  # Your custom first-order loss
print(gamma_inv_dist.second_order_loss_fn(10))  # Your custom second-order loss

Alternative: Monkey-Patching (Quick but Less Maintainable) #

If you want a faster, more lightweight solution without writing wrapper classes, you can monkey-patch the rv_frozen instance directly. Note this affects only the specific frozen distribution you patch, but it’s riskier long-term (may break with SciPy version updates).

Example for Gamma: #

from scipy.stats import gamma

def gamma_first_order_loss(self, x):
    term1 = (self.a - self.scale * x) * self.sf(x)
    term2 = x * self.pdf(x)
    return (term1 + term2) / self.scale

def gamma_second_order_loss(self, x):
    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)

# Monkey-patch the frozen Gamma distribution instance
gamma_frozen = gamma(a=2.5, scale=3.0)
gamma_frozen.first_order_loss_fn = gamma_first_order_loss.__get__(gamma_frozen, type(gamma_frozen))
gamma_frozen.second_order_loss_fn = gamma_second_order_loss.__get__(gamma_frozen, type(gamma_frozen))

# Usage
print(gamma_frozen.first_order_loss_fn(10))

Why Wrappers Are the Better Choice #

• Isolation: Your custom logic stays separate from SciPy’s code, so you don’t have to worry about breaking changes in future SciPy updates.
• Scalability: Adding a new distribution only requires writing a small subclass with its specific loss formulas.
• Clarity: The code is self-documenting—anyone reading it can immediately see which distributions have custom inventory-focused logic.
• Safety: You never risk modifying the behavior of other parts of SciPy that rely on rv_frozen.

内容的提问来源于stack exchange，提问作者Juliana Bright