As you can see below, I am trying to use observed data in some of the inputs for my model but I am not sure my syntax is correct and while I am waiting for PYMC to reinstall on my work computer, I wanted some advice.
I am looking to model the probability curve of firetruck travel times. The model is based on the following function for Expected travel time:
ET = \alpha + \beta\left(\frac{A}{n-\lambda ES}\right)^\gamma
\text{where:}
\begin{align*}
\alpha &= \text{scaling parameter dependednt on city characterisitcs} \\
\beta &= \text{scaling parameter dependednt on city characterisitcs} \\
\gamma &= \text{scaling parameter dependednt on city characterisitcs} \\
\lambda &= \text{expected number of alarms in an hour} \\
A &= \text{physical area of region in square miles }\\
n &= \text{number of fire companies in region} \\
ED &= \text{Expected distance traveled by closest fire company} \\
ES &= \text{Expected total service time in hours} \\
ET &= \text{Expected travel time in minutes of closest fire company}
\end{align*}
With this in mind, my first model is envisioned as:
\begin{align*}
\\
\alpha &\sim \mathcal{normal}(0, 100) \\
\beta &\sim \mathcal{normal}(0, 100) \\
\gamma &\sim |\mathcal{N}(0, 1)| \\
\theta &\sim \mathcal{Gamma}(\tau X) \\
\lambda &= \mathcal{Poisson}(\theta) \\
n &= 33 \\
A &= 61.13 \text{ mi}^2 \\
ES &\sim \mathcal{Student}(\nu=3,\nu =1) \text{?}\\
ET &= \alpha + \beta\left(\frac{A}{n-\lambda ES}\right)^\gamma \\
y &\sim \mathcal{normal}(ET, \sigma)
\end{align*}
Which I have coded below with data which has observed values for Service Times and Travel times.
first_model = pm.Model()
with first_model:
# Priors for unknown model parameters
alpha = pm.Deterministic('alpha', pm.Normal(mu=0, sigma=100))
beta = pm.Deterministic('beta', pm.Normal(mu=0, sigma=100))
gamma = pm.Deterministic('gamma', pm.HalfNormal(sigma=1))
theta = pm.Gamma(mu = pm.Normal(mu=0, sigma=100), sigma=pm.Normal(mu=0, sigma=100), observed = data['travel_times'])
lambda_ = pm.Deterministic('lambda_', pm.Poisson(theta))
ES = pm.Student(3,0,1, observed=data['service_times'])
sigma = pm.HalfNormal(sigma=10)
# Expected value of outcome
ET = alpha + beta*(61.13/(33-lambda_*ES))**gamma
# Sampling distribution of observations
Y_obs = pm.Normal("Y_obs", mu=ET, sigma=pm.Normal(mu=0, sigma=100))
