Let's say I have data on 10 stores that sell widgets, each of which received
num_orders number of orders in a certain timeframe, and sold a total of
total_qty widgets in that timeframe. In other words, for store i, the average quantity per order is
My goal is to use this data produce a prediction interval for the
total_qty variable that accounts for the fact that there are two layers of uncertainty:
- We don't know how many orders a store will get - there is a distribution of
- Even once we know how many orders the store gets, there is uncertainty because we don't know what the quantity for each order will be.
What's the best way to do this?
Let's say the data look like this (in reality, of course, I don't know the distributions that gave rise to the data)
num_stores <- 10 df <- tibble(store_id = 1:num_stores, num_orders = rpois(num_stores, 50), total_qty = rpois(num_stores, 100))
My idea is to use a Bayesian approach:
- Write down a generative probability model
- Fit the model to the data and get samples of parameters from the posterior distribution
- Simulate a large number of data points with the samples from the posterior, and then take the 5th and 95th quantiles of that set of simulated datapoints.
Is this a sensible approach?
As an example, this is the probability model I might start with:
total_qty ~ dpois(lambda_qty*orders) # sum of Poissons is Poisson orders ~ dpois(lambda_orders) lambda_qty ~ dnorm(30, 2) # prior for avg qty ordered lambda_orders ~ dexp(1) # prior for avg num orders
I'm pretty new to Bayesian analysis, so would appreciate any tips.