Scan with sampling

Questions v5
1

I’m hoping to simulate inventory stock levels for a rental product with scan, but I’m very new to using scan.

Basically the pseudocode for the problem is:

initial_stock = 450
stock = initial_stock
returns = np.zeros(days)
for day in days:
    stock += returns[day]
    demand ~ Poisson(lambda)
    rentals = where(stock < demand, stock, demand)
    new_returns ~ Lognormal(params, size=rentals)
    stock -= rentals
    # new return times would have to get transformed to be an array of size ndays
    # where each index is a date, and the value is the number of returns that date
    returns += new_returns

I’m hoping to simulate data from that process, and I’m hoping to fit a model on the data to learn:

  1. Expected stock levels each day
  2. demand
  3. typical time an item is rented out before its returned

I have a basic model where returns come back at some constant rate below, but I want to change the model to reflect the pseudocode above. I’m hoping

SEED = 99
np.random.seed(SEED)

n_days = 10
stock0 = 450

coords={
    "date":np.arange(n_days),
}

with pm.Model(coords=coords) as sim:
    
    # assume constant demand rate
    demand = pm.Poisson("demand", 200, dims=("date"))   
    # assume constant return rate
    returns = pm.Poisson("returns", 120, dims=("date"))
    
    s0 = pm.MutableData(name="s0", value=np.zeros(n_days), dims=("date"))
    s0 = pt.set_subtensor(s0[0], stock0)
    
    def stock_position(returns_t, demand_t, stock_tm1):
        rentals = pt.where(stock_tm1 < demand_t, stock_tm1, demand_t)
        stock = pt.clip( stock_tm1 + returns_t - demand_t, 0, np.inf)
        # what I was hoping might work, but Im not sure how i'd fit parameters to that
        # return_times = pt.eq(
        #    pm.LogNormal.dist(1.25, 0.5, size=rentals).ceil(),
        #    (pt.arange(ndays))[:,None]
        # ).sum(1)
        return stock, rentals
    
    [outputs, rental_counts], _ = pytensor.scan(
        fn=stock_position, 
        sequences=[
            dict(input=returns, taps=[0]),
            dict(input=demand, taps=[0]),
        ],
        outputs_info=[dict(initial=s0, taps=[-1]), None],

        n_steps=n_days-1,
        strict=True
    )
    
    rental_counts = pm.Deterministic("rental_counts", rental_counts[:,0])
    stock = pm.Deterministic("stock", var=pt.set_subtensor(s0[1:], outputs[:,0]), dims="date")
    
    # sample prior predictive
    pp = pm.sample_prior_predictive()

I’m wondering if there’s a way to do this all cleanly in one model with scan, or if I’m better off just fitting separate models and simulating stock levels from those

2

Spent alot more time on this, and I think I’m approaching the problem wrong in general.

I think this could be a good fit for a hidden markov model with time varying transitions, where there are 2 states, rented and not rented, something like below, but the transition probabilities vary over time

\begin{equation*} \begin{bmatrix} 1- P_{rent} & P_{rent} \\ P_{return} & 1-P_{return} \end{bmatrix} \end{equation*}

Probability of being rented would change over time depending on demand and number of other units in stock (seems difficult, I’d have to keep track of all of the other states), and probability of being returned would be based on a lognormal distribution (where time to return from rental start is lognormally distributed).

This example seems similar.

I’m probably going to take the simple approach here and model demand and rental time separately, but if anyone has any pointers for the time varying markov model approach above I’d love to hear it