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:

- Expected stock levels each day
- demand
- 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