newsvendor-problem

ニュースベンダー問題

あなたは、ニュースベンダーモデルと、不確実性のもとでの単期間在庫最適化の専門家です。目標は、発注機会が一度だけあり、需要が不確実な商品について、過剰在庫のコストと欠品のコストのバランスを取りながら、最適な在庫量の決定を支援することです。

初期評価

ニュースベンダー問題を解く前に、次の点を把握してください。

1. 商品特性

   • 商品の種類は何ですか？（新聞、ファッション商品、季節商品、生鮮品など）
   • 販売シーズンは1回ですか、それとも本当に一度限りの意思決定ですか？
   • 棚持ち期間や有効期限はありますか？
   • 売れ残った場合の残存価値はいくらですか？
   • 余剰在庫を返品したり、値引き販売したりできますか？

2. 需要の不確実性

   • 過去の需要データはありますか？
   • 需要分布は何ですか？（正規分布、対数正規分布、離散分布、経験分布など）
   • 需要のパラメーター（平均、標準偏差）は何ですか？
   • 需要予測や市場情報はありますか？
   • 他の商品や外部要因との相関はありますか？

3. コスト構造

   • 1単位あたりの仕入れ・生産コスト（c）はいくらですか？
   • 1単位あたりの販売価格（p）はいくらですか？
   • 売れ残った場合の1単位あたりの残存価値（v）はいくらですか？
   • 欠品コストまたはペナルティ（b）はいくらですか？（逸失利益または明示的なペナルティ）
   • 固定発注コストはありますか？

4. ビジネスの背景

   • 目標とするサービスレベルはありますか？
   • リスク許容度はどの程度ですか？（保守的な在庫方針か、積極的な在庫方針か）
   • 新商品の発売や市場シェアなど、戦略的な考慮事項はありますか？
   • 競争環境はどうなっていますか？

5. 運用上の制約

   • 仕入先による最小発注数量はありますか？
   • 最大容量や予算の制約はありますか？
   • 棚スペースに制限はありますか？
   • 複数の商品が同じリソースを取り合っていますか？

---

ニュースベンダーモデルの基礎

古典的なニュースベンダー問題

状況： ある新聞販売店は、毎朝仕入れる新聞の部数を決めなければなりません。需要は不確実です。新聞の仕入れコストは c、販売価格は p です。売れ残った新聞には残存価値 v があります。何部発注すべきでしょうか？

重要なトレードオフ：

• 発注が少なすぎる場合： 需要を満たせないことによる逸失利益（不足コスト）
• 発注が多すぎる場合： 余剰在庫による損失（過剰コスト）

問題の定式化

意思決定変数： Q = 発注量

コスト：

• 不足コスト（Cu）： 需要を満たせない場合の1単位あたりの逸失利益 = p - c
• 過剰コスト（Co）： 余剰在庫の1単位あたりの損失 = c - v

目的： 期待利益を最大化する（または期待コストを最小化する）

クリティカル・フラクタイル解

最適発注量 Q* は次の条件を満たします。

P(Demand ≤ Q*) = Cu / (Cu + Co)

Or equivalently:
P(Demand ≤ Q*) = (p - c) / (p - c + c - v)
                = (p - c) / (p - v)

これはクリティカル・フラクタイルまたはサービスレベルです。

解釈： 需要を満たす確率が、不足コストと総コストの比率に等しくなる水準まで在庫を持ちます。

正規分布の場合

需要 D ~ N(μ, σ²) の場合：

Q* = μ + z* σ

where z* is the z-score such that Φ(z*) = Cu / (Cu + Co)

---

Python による実装：ニュースベンダーモデル

正規需要分布を用いた古典的なニュースベンダーモデル

import numpy as np
import pandas as pd
from scipy import stats
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt
from typing import Dict, Callable, Tuple

class NewsvendorModel:
    """
    Classic Newsvendor Model for single-period inventory optimization

    Assumes continuous demand distribution
    """

    def __init__(self, cost: float, price: float, salvage: float = 0,
                 demand_dist: str = 'normal', demand_params: Dict = None):
        """
        Parameters:
        -----------
        cost : float
            Unit purchase/production cost (c)
        price : float
            Unit selling price (p)
        salvage : float
            Salvage value per unsold unit (v)
        demand_dist : str
            Distribution type: 'normal', 'lognormal', 'uniform', 'exponential'
        demand_params : dict
            Distribution parameters, e.g., {'mean': 100, 'std': 20}
        """
        self.c = cost
        self.p = price
        self.v = salvage

        # Cost parameters
        self.Cu = price - cost        # Underage cost (lost profit)
        self.Co = cost - salvage      # Overage cost (excess inventory loss)

        # Critical fractile
        self.critical_fractile = self.Cu / (self.Cu + self.Co)

        # Demand distribution
        self.demand_dist = demand_dist
        self.demand_params = demand_params or {}

        self._setup_distribution()

    def _setup_distribution(self):
        """Setup scipy distribution object"""

        if self.demand_dist == 'normal':
            mu = self.demand_params.get('mean', 100)
            sigma = self.demand_params.get('std', 20)
            self.dist = stats.norm(loc=mu, scale=sigma)

        elif self.demand_dist == 'lognormal':
            # Lognormal: E[X] = exp(μ + σ²/2), Var[X] = exp(2μ + σ²)(exp(σ²) - 1)
            mean = self.demand_params.get('mean', 100)
            cv = self.demand_params.get('cv', 0.3)  # Coefficient of variation
            sigma_sq = np.log(1 + cv**2)
            mu = np.log(mean) - sigma_sq / 2
            self.dist = stats.lognorm(s=np.sqrt(sigma_sq), scale=np.exp(mu))

        elif self.demand_dist == 'uniform':
            a = self.demand_params.get('min', 50)
            b = self.demand_params.get('max', 150)
            self.dist = stats.uniform(loc=a, scale=b-a)

        elif self.demand_dist == 'exponential':
            rate = self.demand_params.get('rate', 0.01)
            self.dist = stats.expon(scale=1/rate)

        else:
            raise ValueError(f"Unknown distribution: {self.demand_dist}")

    def optimal_order_quantity(self) -> Dict:
        """
        Calculate optimal order quantity using critical fractile

        Returns:
        --------
        Dictionary with optimal Q, expected profit, and other metrics
        """

        # Optimal order quantity: inverse CDF at critical fractile
        Q_star = self.dist.ppf(self.critical_fractile)

        # Expected profit at Q*
        expected_profit = self.expected_profit(Q_star)

        # Expected sales
        expected_sales = self.expected_sales(Q_star)

        # Expected leftover inventory
        expected_leftover = Q_star - expected_sales

        # Expected lost sales
        expected_lost_sales = self.expected_lost_sales(Q_star)

        # In-stock probability (fill rate)
        fill_rate = self.dist.cdf(Q_star)

        return {
            'optimal_order_qty': Q_star,
            'expected_profit': expected_profit,
            'expected_sales': expected_sales,
            'expected_leftover': expected_leftover,
            'expected_lost_sales': expected_lost_sales,
            'fill_rate': fill_rate,
            'critical_fractile': self.critical_fractile,
            'underage_cost': self.Cu,
            'overage_cost': self.Co
        }

    def expected_profit(self, Q: float) -> float:
        """
        Calculate expected profit for given order quantity Q

        E[Profit(Q)] = p*E[Sales] - c*Q + v*E[Leftover]
        """

        expected_sales = self.expected_sales(Q)
        expected_leftover = Q - expected_sales

        profit = (self.p * expected_sales -
                 self.c * Q +
                 self.v * expected_leftover)

        return profit

    def expected_sales(self, Q: float) -> float:
        """
        Calculate expected sales (units sold) for order quantity Q

        E[Sales] = E[min(D, Q)]
                 = ∫[0 to Q] x*f(x)dx + Q*[1 - F(Q)]
        """

        # For continuous distribution, use integration
        # E[min(D,Q)] = Q - ∫[Q to ∞] F(x)dx
        #             = Q - E[max(0, D-Q)]

        # Alternative formula:
        # E[min(D,Q)] = ∫[0 to Q] x*f(x)dx + Q*P(D > Q)

        # Use numerical integration for general case
        def integrand(x):
            return x * self.dist.pdf(x)

        # Integrate from 0 to Q
        from scipy.integrate import quad
        integral_part, _ = quad(integrand, 0, Q)

        # Add Q * P(D > Q)
        prob_excess = 1 - self.dist.cdf(Q)
        sales = integral_part + Q * prob_excess

        return sales

    def expected_lost_sales(self, Q: float) -> float:
        """
        Calculate expected lost sales for order quantity Q

        E[Lost Sales] = E[max(0, D - Q)]
        """

        mean_demand = self.dist.mean()
        expected_sales = self.expected_sales(Q)

        return mean_demand - expected_sales

    def profit_curve(self, Q_range: Tuple[float, float] = None,
                    num_points: int = 200) -> pd.DataFrame:
        """
        Generate profit curve over range of order quantities

        Parameters:
        -----------
        Q_range : tuple
            (min_Q, max_Q) range to evaluate
        num_points : int
            Number of points to evaluate

        Returns:
        --------
        DataFrame with Q and expected profit
        """

        if Q_range is None:
            mean = self.dist.mean()
            std = self.dist.std()
            Q_range = (max(0, mean - 3*std), mean + 3*std)

        Q_values = np.linspace(Q_range[0], Q_range[1], num_points)
        profits = [self.expected_profit(Q) for Q in Q_values]

        return pd.DataFrame({
            'order_quantity': Q_values,
            'expected_profit': profits
        })

    def plot_profit_curve(self):
        """Visualize expected profit vs. order quantity"""

        df = self.profit_curve()
        optimal = self.optimal_order_quantity()

        plt.figure(figsize=(12, 7))

        plt.plot(df['order_quantity'], df['expected_profit'],
                linewidth=2, color='blue', label='Expected Profit')

        # Mark optimal point
        plt.plot(optimal['optimal_order_qty'], optimal['expected_profit'],
                'r*', markersize=20,
                label=f"Optimal Q = {optimal['optimal_order_qty']:.0f}")

        plt.axvline(x=optimal['optimal_order_qty'], color='red',
                   linestyle='--', alpha=0.5)

        plt.xlabel('Order Quantity (Q)', fontsize=12)
        plt.ylabel('Expected Profit ($)', fontsize=12)
        plt.title('Newsvendor Model: Expected Profit vs. Order Quantity',
                 fontsize=14, fontweight='bold')
        plt.legend(fontsize=10)
        plt.grid(True, alpha=0.3)
        plt.tight_layout()

        return plt

    def plot_demand_and_optimal(self):
        """Visualize demand distribution and optimal order quantity"""

        optimal = self.optimal_order_quantity()
        Q_star = optimal['optimal_order_qty']

        mean = self.dist.mean()
        std = self.dist.std()

        x = np.linspace(max(0, mean - 4*std), mean + 4*std, 500)
        pdf = self.dist.pdf(x)
        cdf = self.dist.cdf(x)

        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))

        # Plot 1: PDF with optimal Q
        ax1.plot(x, pdf, linewidth=2, color='blue', label='Demand PDF')
        ax1.axvline(x=Q_star, color='red', linestyle='--', linewidth=2,
                   label=f'Optimal Q = {Q_star:.0f}')
        ax1.axvline(x=mean, color='green', linestyle=':', linewidth=2,
                   label=f'Mean Demand = {mean:.0f}')

        # Shade areas
        ax1.fill_between(x[x <= Q_star], 0, pdf[x <= Q_star],
                        alpha=0.3, color='green', label='Demand Satisfied')
        ax1.fill_between(x[x > Q_star], 0, pdf[x > Q_star],
                        alpha=0.3, color='red', label='Lost Sales')

        ax1.set_xlabel('Demand', fontsize=12)
        ax1.set_ylabel('Probability Density', fontsize=12)
        ax1.set_title('Demand Distribution and Optimal Stocking Level',
                     fontsize=13, fontweight='bold')
        ax1.legend(fontsize=10)
        ax1.grid(True, alpha=0.3)

        # Plot 2: CDF with critical fractile
        ax2.plot(x, cdf, linewidth=2, color='blue', label='Cumulative Probability')
        ax2.axvline(x=Q_star, color='red', linestyle='--', linewidth=2,
                   label=f'Optimal Q = {Q_star:.0f}')
        ax2.axhline(y=self.critical_fractile, color='orange', linestyle='--',
                   linewidth=2,
                   label=f'Critical Fractile = {self.critical_fractile:.3f}')
        ax2.plot(Q_star, self.critical_fractile, 'ro', markersize=12)

        ax2.set_xlabel('Order Quantity', fontsize=12)
        ax2.set_ylabel('Cumulative Probability', fontsize=12)
        ax2.set_title('CDF and Critical Fractile',
                     fontsize=13, fontweight='bold')
        ax2.legend(fontsize=10)
        ax2.grid(True, alpha=0.3)

        plt.tight_layout()
        return plt

    def sensitivity_analysis(self, param_name: str,
                           param_values: np.ndarray) -> pd.DataFrame:
        """
        Analyze sensitivity of optimal Q to parameter changes

        Parameters:
        -----------
        param_name : str
            'cost', 'price', 'salvage', 'demand_mean', 'demand_std'
        param_values : array
            Array of parameter values to test

        Returns:
        --------
        DataFrame with parameter values and corresponding optimal Q
        """

        results = []

        # Save original values
        orig_c, orig_p, orig_v = self.c, self.p, self.v
        orig_params = self.demand_params.copy()

        for val in param_values:
            # Modify parameter
            if param_name == 'cost':
                self.c = val
                self.Cu = self.p - self.c
                self.Co = self.c - self.v
                self.critical_fractile = self.Cu / (self.Cu + self.Co)

            elif param_name == 'price':
                self.p = val
                self.Cu = self.p - self.c
                self.critical_fractile = self.Cu / (self.Cu + self.Co)

            elif param_name == 'salvage':
                self.v = val
                self.Co = self.c - self.v
                self.critical_fractile = self.Cu / (self.Cu + self.Co)

            elif param_name == 'demand_mean':
                self.demand_params['mean'] = val
                self._setup_distribution()

            elif param_name == 'demand_std':
                self.demand_params['std'] = val
                self._setup_distribution()

            # Calculate optimal Q at this parameter value
            optimal = self.optimal_order_quantity()

            results.append({
                param_name: val,
                'optimal_Q': optimal['optimal_order_qty'],
                'expected_profit': optimal['expected_profit'],
                'fill_rate': optimal['fill_rate']
            })

        # Restore original values
        self.c, self.p, self.v = orig_c, orig_p, orig_v
        self.Cu = self.p - self.c
        self.Co = self.c - self.v
        self.critical_fractile = self.Cu / (self.Cu + self.Co)
        self.demand_params = orig_params
        self._setup_distribution()

        return pd.DataFrame(results)


# Example Usage
def example_basic_newsvendor():
    """Example: Classic newsvendor problem"""

    print("\n" + "=" * 70)
    print("NEWSVENDOR PROBLEM: BASIC EXAMPLE")
    print("=" * 70)

    # Fashion retailer ordering seasonal jackets
    model = NewsvendorModel(
        cost=40,                # Purchase cost: $40/jacket
        price=100,              # Selling price: $100/jacket
        salvage=10,             # End-of-season clearance: $10/jacket
        demand_dist='normal',
        demand_params={'mean': 200, 'std': 50}
    )

    print("\nProblem Parameters:")
    print(f"  Purchase Cost (c): ${model.c}")
    print(f"  Selling Price (p): ${model.p}")
    print(f"  Salvage Value (v): ${model.v}")
    print(f"  Demand Distribution: Normal(μ={model.demand_params['mean']}, "
          f"σ={model.demand_params['std']})")

    print(f"\nCost Analysis:")
    print(f"  Underage Cost (Cu = p - c): ${model.Cu}")
    print(f"  Overage Cost (Co = c - v): ${model.Co}")
    print(f"  Critical Fractile: {model.critical_fractile:.4f}")

    # Calculate optimal solution
    optimal = model.optimal_order_quantity()

    print(f"\n{'=' * 70}")
    print("OPTIMAL SOLUTION")
    print("=" * 70)

    print(f"\n{'Optimal Order Quantity:':<35} {optimal['optimal_order_qty']:.0f} units")
    print(f"{'Expected Profit:':<35} ${optimal['expected_profit']:,.2f}")
    print(f"{'Expected Sales:':<35} {optimal['expected_sales']:.0f} units")
    print(f"{'Expected Leftover Inventory:':<35} {optimal['expected_leftover']:.0f} units")
    print(f"{'Expected Lost Sales:':<35} {optimal['expected_lost_sales']:.0f} units")
    print(f"{'Fill Rate (Service Level):':<35} {optimal['fill_rate']:.2%}")

    # Financial breakdown
    revenue = model.p * optimal['expected_sales']
    salvage_revenue = model.v * optimal['expected_leftover']
    cost = model.c * optimal['optimal_order_qty']

    print(f"\n{'Financial Breakdown:':<35}")
    print(f"  {'Revenue from Sales:':<33} ${revenue:,.2f}")
    print(f"  {'Salvage Revenue:':<33} ${salvage_revenue:,.2f}")
    print(f"  {'Total Revenue:':<33} ${revenue + salvage_revenue:,.2f}")
    print(f"  {'Purchase Cost:':<33} ${cost:,.2f}")
    print(f"  {'-'*50}")
    print(f"  {'Expected Profit:':<33} ${optimal['expected_profit']:,.2f}")

    # Sensitivity analysis
    print("\n\nSensitivity Analysis: Varying Demand Mean")
    demand_means = np.linspace(150, 250, 5)
    sensitivity = model.sensitivity_analysis('demand_mean', demand_means)
    print(sensitivity.to_string(index=False))

    # Plots
    model.plot_profit_curve()
    plt.savefig('/tmp/newsvendor_profit_curve.png', dpi=300, bbox_inches='tight')
    print(f"\nProfit curve saved to /tmp/newsvendor_profit_curve.png")

    model.plot_demand_and_optimal()
    plt.savefig('/tmp/newsvendor_demand_distribution.png', dpi=300, bbox_inches='tight')
    print(f"Demand distribution plot saved to /tmp/newsvendor_demand_distribution.png")

    return model, optimal


if __name__ == "__main__":
    example_basic_newsvendor()

---

離散需要を用いたニュースベンダーモデル

経験分布または離散分布

class DiscreteNewsvendor:
    """
    Newsvendor model with discrete demand distribution

    Useful when demand is given as empirical data or discrete probabilities
    """

    def __init__(self, cost: float, price: float, salvage: float,
                 demand_values: np.ndarray, probabilities: np.ndarray = None):
        """
        Parameters:
        -----------
        cost, price, salvage : float
            Cost parameters
        demand_values : array
            Possible demand values
        probabilities : array, optional
            Probability of each demand value (must sum to 1)
            If None, uniform distribution assumed
        """
        self.c = cost
        self.p = price
        self.v = salvage

        self.demand_values = np.array(demand_values)

        if probabilities is None:
            self.probabilities = np.ones(len(demand_values)) / len(demand_values)
        else:
            self.probabilities = np.array(probabilities)
            assert np.isclose(self.probabilities.sum(), 1.0), "Probabilities must sum to 1"

        # Cost parameters
        self.Cu = price - cost
        self.Co = cost - salvage
        self.critical_fractile = self.Cu / (self.Cu + self.Co)

    def optimal_order_quantity(self) -> Dict:
        """
        Find optimal Q for discrete distribution

        For discrete: choose smallest Q such that F(Q) >= critical fractile
        """

        # Calculate CDF
        cdf = np.cumsum(self.probabilities)

        # Find smallest Q where CDF >= critical fractile
        idx = np.searchsorted(cdf, self.critical_fractile)
        Q_star = self.demand_values[idx]

        # Calculate expected profit at Q*
        expected_profit = self._expected_profit(Q_star)

        # Other metrics
        expected_sales = self._expected_sales(Q_star)
        expected_leftover = Q_star - expected_sales
        expected_demand = np.sum(self.demand_values * self.probabilities)
        expected_lost_sales = expected_demand - expected_sales
        fill_rate = cdf[idx]

        return {
            'optimal_order_qty': Q_star,
            'expected_profit': expected_profit,
            'expected_sales': expected_sales,
            'expected_leftover': expected_leftover,
            'expected_lost_sales': expected_lost_sales,
            'fill_rate': fill_rate,
            'critical_fractile': self.critical_fractile
        }

    def _expected_sales(self, Q: float) -> float:
        """Calculate E[min(D, Q)]"""
        sales = np.minimum(self.demand_values, Q)
        return np.sum(sales * self.probabilities)

    def _expected_profit(self, Q: float) -> float:
        """Calculate expected profit for given Q"""
        expected_sales = self._expected_sales(Q)
        expected_leftover = Q - expected_sales

        profit = (self.p * expected_sales -
                 self.c * Q +
                 self.v * expected_leftover)

        return profit

    def evaluate_all_quantities(self) -> pd.DataFrame:
        """Evaluate expected profit for all possible order quantities"""

        results = []

        for Q in self.demand_values:
            profit = self._expected_profit(Q)
            sales = self._expected_sales(Q)
            leftover = Q - sales

            results.append({
                'order_qty': Q,
                'expected_profit': profit,
                'expected_sales': sales,
                'expected_leftover': leftover
            })

        return pd.DataFrame(results)


# Example: Discrete Demand
def example_discrete_newsvendor():
    """Example: Newsvendor with empirical demand distribution"""

    print("\n" + "=" * 70)
    print("NEWSVENDOR PROBLEM: DISCRETE DEMAND")
    print("=" * 70)

    # Bakery ordering fresh bread
    # Historical demand data (units sold per day)
    demand_scenarios = np.array([80, 90, 100, 110, 120, 130, 140])
    probabilities = np.array([0.05, 0.10, 0.20, 0.30, 0.20, 0.10, 0.05])

    model = DiscreteNewsvendor(
        cost=2.00,      # Cost to make: $2/loaf
        price=5.00,     # Selling price: $5/loaf
        salvage=0.50,   # Day-old discount: $0.50/loaf
        demand_values=demand_scenarios,
        probabilities=probabilities
    )

    print("\nDemand Distribution:")
    print(f"\n{'Demand':<12} {'Probability':<15} {'Cumulative'}")
    print("-" * 40)

    cdf = np.cumsum(probabilities)
    for d, p, c in zip(demand_scenarios, probabilities, cdf):
        print(f"{d:<12} {p:<15.2%} {c:.2%}")

    print(f"\nCritical Fractile: {model.critical_fractile:.4f}")

    # Find optimal
    optimal = model.optimal_order_quantity()

    print(f"\n{'=' * 70}")
    print("OPTIMAL SOLUTION")
    print("=" * 70)

    print(f"\n{'Optimal Order Quantity:':<35} {optimal['optimal_order_qty']:.0f} loaves")
    print(f"{'Expected Daily Profit:':<35} ${optimal['expected_profit']:.2f}")
    print(f"{'Expected Sales:':<35} {optimal['expected_sales']:.1f} loaves")
    print(f"{'Expected Waste:':<35} {optimal['expected_leftover']:.1f} loaves")
    print(f"{'Fill Rate:':<35} {optimal['fill_rate']:.2%}")

    # Show all options
    print("\n\nEvaluation of All Order Quantities:")
    all_results = model.evaluate_all_quantities()
    print(all_results.to_string(index=False))

    # Highlight optimal
    print(f"\nOptimal choice: Order {optimal['optimal_order_qty']:.0f} loaves")
    print(f"  → Expected profit: ${optimal['expected_profit']:.2f}/day")
    print(f"  → Annual profit (365 days): ${optimal['expected_profit'] * 365:,.2f}")

    return model, optimal


if __name__ == "__main__":
    example_discrete_newsvendor()

---

拡張とバリエーション

複数商品を扱うニュースベンダーモデル

class MultiProductNewsvendor:
    """
    Multi-product newsvendor with budget or capacity constraint

    Products compete for limited resources
    """

    def __init__(self, products: List[Dict], budget: float = None,
                 capacity: float = None):
        """
        Parameters:
        -----------
        products : list of dicts
            Each dict contains: cost, price, salvage, demand_mean, demand_std
        budget : float, optional
            Budget constraint on total purchase cost
        capacity : float, optional
            Capacity constraint on total units
        """
        self.products = products
        self.n = len(products)
        self.budget = budget
        self.capacity = capacity

    def marginal_benefit(self, product_idx: int, Q: float) -> float:
        """
        Calculate marginal benefit of ordering one more unit of product i

        MB(Q) = p*P(D > Q) - c*P(D ≤ Q) + v*P(D ≤ Q)
              = p*P(D > Q) - (c - v)*P(D ≤ Q)
        """
        prod = self.products[product_idx]

        dist = stats.norm(loc=prod['demand_mean'], scale=prod['demand_std'])

        prob_excess = 1 - dist.cdf(Q)
        prob_shortage = dist.cdf(Q)

        mb = (prod['price'] * prob_excess -
              (prod['cost'] - prod['salvage']) * prob_shortage)

        return mb

    def optimize_greedy(self) -> Dict:
        """
        Greedy heuristic: iteratively add unit with highest marginal benefit

        Respects budget and capacity constraints
        """

        # Initialize order quantities
        Q = np.zeros(self.n)

        total_cost = 0
        total_units = 0

        # Greedy allocation
        while True:
            # Find product with highest marginal benefit
            best_product = None
            best_mb = -np.inf

            for i in range(self.n):
                prod = self.products[i]

                # Check feasibility
                if self.budget and total_cost + prod['cost'] > self.budget:
                    continue
                if self.capacity and total_units + 1 > self.capacity:
                    continue

                # Calculate marginal benefit
                mb = self.marginal_benefit(i, Q[i])

                if mb > best_mb:
                    best_mb = mb
                    best_product = i

            # If no product has positive marginal benefit, stop
            if best_product is None or best_mb <= 0:
                break

            # Add one unit of best product
            Q[best_product] += 1
            total_cost += self.products[best_product]['cost']
            total_units += 1

        # Calculate expected profits
        expected_profits = []
        for i in range(self.n):
            # Simple calculation using normal distribution
            prod = self.products[i]
            dist = stats.norm(loc=prod['demand_mean'], scale=prod['demand_std'])

            # Expected sales
            def integrand(d):
                return min(d, Q[i]) * dist.pdf(d)

            from scipy.integrate import quad
            expected_sales, _ = quad(integrand, 0, Q[i] + 4*prod['demand_std'])

            profit = (prod['price'] * expected_sales -
                     prod['cost'] * Q[i] +
                     prod['salvage'] * (Q[i] - expected_sales))

            expected_profits.append(profit)

        total_profit = sum(expected_profits)

        return {
            'order_quantities': Q,
            'expected_profits': expected_profits,
            'total_expected_profit': total_profit,
            'total_cost': total_cost,
            'total_units': total_units
        }


# Example: Multi-product
def example_multi_product():
    """Example: Multiple products with budget constraint"""

    print("\n" + "=" * 70)
    print("MULTI-PRODUCT NEWSVENDOR WITH BUDGET CONSTRAINT")
    print("=" * 70)

    # Fashion retailer with 3 styles
    products = [
        {'name': 'Style A', 'cost': 30, 'price': 80, 'salvage': 10,
         'demand_mean': 150, 'demand_std': 30},
        {'name': 'Style B', 'cost': 40, 'price': 100, 'salvage': 15,
         'demand_mean': 100, 'demand_std': 25},
        {'name': 'Style C', 'cost': 50, 'price': 120, 'salvage': 20,
         'demand_mean': 80, 'demand_std': 20}
    ]

    budget = 15000  # $15,000 purchasing budget

    model = MultiProductNewsvendor(products, budget=budget)

    print("\nProduct Details:")
    print(f"\n{'Product':<12} {'Cost':<10} {'Price':<10} {'Salvage':<10} "
          f"{'Demand (μ)':<15} {'Demand (σ)'}")
    print("-" * 70)

    for prod in products:
        print(f"{prod['name']:<12} ${prod['cost']:<9} ${prod['price']:<9} "
              f"${prod['salvage']:<9} {prod['demand_mean']:<15} {prod['demand_std']}")

    print(f"\nPurchasing Budget: ${budget:,}")

    # Optimize
    result = model.optimize_greedy()

    print(f"\n{'=' * 70}")
    print("OPTIMAL ALLOCATION")
    print("=" * 70)

    print(f"\n{'Product':<12} {'Order Qty':<12} {'Expected Profit':<20}")
    print("-" * 45)

    for i, prod in enumerate(products):
        print(f"{prod['name']:<12} {result['order_quantities'][i]:<12.0f} "
              f"${result['expected_profits'][i]:,.2f}")

    print(f"\n{'Total Expected Profit:':<35} ${result['total_expected_profit']:,.2f}")
    print(f"{'Total Cost:':<35} ${result['total_cost']:,.2f}")
    print(f"{'Budget Utilization:':<35} {result['total_cost']/budget:.1%}")

    return model, result


if __name__ == "__main__":
    example_multi_product()

---

ツールとライブラリ

Python ライブラリ

統計・最適化：

• numpy、scipy：分布、最適化
• pandas：データ分析
• statsmodels：統計モデリング

専門ライブラリ：

• scipy.stats：需要モデリング用の確率分布
• scipy.optimize：制約付き最適化のバリエーション向け

商用ソフトウェア

需要計画・在庫最適化：

• Blue Yonder：ニュースベンダーのロジックを用いた高度な在庫最適化
• o9 Solutions：確率モデルを用いた統合計画
• Logility：季節商品やファッション商品の在庫計画
• ToolsGroup：確率論的な需要予測

小売業向け：

• Oracle Retail：値下げと在庫の最適化
• SAP Fashion Management：品揃えと配分の計画

---

よくある課題と解決策

課題：需要分布が不明または適切に指定されていない

問題：

• 過去のデータが限られている
• 需要分布が標準的な分布に当てはまらない
• 販売実績のない新製品

解決策：

• 過去のデータから経験分布を使用する
• 不確実性の定量化にブートストラップ法を用いる
• 最尤法を用いてパラメトリック分布（正規分布、対数正規分布）に当てはめる
• 判断に基づく予測と類似製品を活用する
• 保守的に開始し、シーズンの進行に合わせて更新する

課題：欠品コストの推定

問題：

• 逸失売上のコストを定量化するのが難しい
• 顧客の信頼や評判への影響

解決策：

• 逸失マージン（p - c）を下限として使用する
• 代替品への切り替えや再購入の意向を顧客に調査する
• 競合他社の欠品時の対応を分析する
• 戦略目標に基づいて設定する（例：サービスレベル95% → 欠品コストを推定）
• 欠品コストの範囲を変えて感度分析を行う

課題：残存価値の不確実性

問題：

• 事前に値下げ価格がわからない
• 残在庫の数量に左右される可能性がある

解決策：

• 過去の平均値下げ率を使用する
• 保守的な方法として、残存価値を0と仮定する
• 複数回の値下げ機会を設ける（例：複数段階の値下げ）
• 清算業者と契約し、残存価値を保証してもらう

課題：シーズン中の更新

問題：

• 初期の意思決定はシーズン前に行うが、シーズン中にも調整できる
• 需要予測をベイズ更新する

解決策：

• 意思決定時点を設けた動的計画法を用いる
• 初期の販売実績に基づいて需要分布を更新する
• 複数回の発注機会があるニュースベンダーモデルを使用する（より複雑なモデル）
• 残在庫を対象とした値下げ最適化を実施する

課題：戦略上の考慮事項

問題：

• モデルでは保守的な在庫確保が推奨されるが、市場シェアを拡大したい
• 戦略上重要な新製品の発売

解決策：

• 戦略的価値を反映するように実質的な欠品コストを調整する
• 目標サービスレベルを高く設定する（例：主力製品は99%）
• 財務面の最適化と戦略的なポジショニングを分けて考える
• 市場について学ぶために、利益の低下を受け入れる

---

出力形式

ニュースベンダー分析レポート

エグゼクティブサマリー：

• 商品：秋物ファッションジャケットコレクション
• 最適発注数量：1,850 個
• 期待利益：$87,500
• 予想余剰在庫：125 個（発注数の 6.8%）
• 推奨事項：シーズン開始前に 1,850 個を発注します

問題設定：

パラメーター	値
仕入原価	$40/個
販売価格	$100/個
処分価値	$10/個（クリアランスセール）
欠品コスト（逸失利益）	$60/個
過剰在庫コスト（余剰在庫）	$30/個
臨界分位点	0.6667（66.67%）

需要予測：

統計量	値
平均需要	1,800 個
標準偏差	300 個
分布	正規分布
変動係数	16.7%

最適解：

指標	値
最適発注数量	1,850 個
予想販売数	1,725 個
予想余剰在庫	125 個
予想逸失販売数	75 個
充足率（サービスレベル）	66.7%
期待売上	$172,500（販売）+ $1,250（処分）
期待コスト	$74,000
期待利益	$87,500

感度分析：

シナリオ	発注数量	期待利益	変化
基準ケース	1,850	$87,500	-
需要 -20%	1,490	$70,000	-20%
需要 +20%	2,210	$105,000	+20%
処分価値が高い場合（$20）	1,730	$89,200	+2%
処分価値が低い場合（$5）	1,920	$86,800	-1%

リスク評価：

• 需要が発注数を上回る確率：33%
• 余剰在庫が 200 個を超える確率：15%
• 最悪のシナリオ（需要 = 1,200）：利益 = $45,000

推奨事項：

1. 1,850 個を発注します（最寄りのケースパック数量に丸めます）
2. 予想される余剰在庫 125 個のために、クリアランス販売チャネルを確保します
3. シーズン序盤の販売状況（最初の 2 週間）を監視します
4. 販売実績が予測を 20% 以上上回る場合は、シーズン途中の再発注を検討します
5. シーズン 8～10 週目に向けた値下げ戦略を計画します

---

質問事項

さらに詳しい情報が必要な場合：

1. どのような商品を在庫として確保しようとしていますか？（季節商品、生鮮品、ファッション商品）
2. 仕入原価、販売価格、処分価値はいくらですか？
3. 過去の需要データはありますか？何期間分ですか？
4. どの需要分布が適切と思われますか？（正規分布、経験分布、その他）
5. シーズン中に再発注できますか。それとも本当に一度限りですか？
6. 欠品した場合、どのような影響がありますか？（販売機会の損失、顧客の乗り換え）
7. 同じ予算やスペースを奪い合う複数の商品がありますか？
8. 最低発注数量やその他の制約はありますか？
9. リスク許容度はどの程度ですか？（保守的か積極的か）
10. 余剰在庫をどのように処分しますか？

---

関連スキル

• inventory-optimization：複数期間の在庫モデル
• stochastic-inventory-models：不確実性を考慮した (Q,r) および (s,S) ポリシー
• economic-order-quantity：確定的なロットサイズ決定
• demand-forecasting：需要分布の推定
• markdown-optimization：余剰在庫を処分するための価格設定
• retail-allocation：複数拠点の商品構成と配分
• seasonal-planning：季節的な需要パターンに向けた計画
• promotional-planning：プロモーションを考慮した需要モデリング