Mesquite

Predicting the yields of mesquite bushes. See Chapter 12 in Regression and Other Stories.

Source: Mesquite/mesquite.Rmd.

The Python version keeps data handling explicit and uses lapylace for Stan-backed generalized linear models, so the statistical model can be read from a formula rather than from handwritten Stan.

from pathlib import Path
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lapylace as lp

sys.path.append(str(Path.cwd().parent / 'python'))
from data import ros_path

Data

mesquite = pd.read_csv(ros_path('Mesquite/data', 'mesquite.dat'), sep=r'\s+')
mesquite.head()
obs group diam1 diam2 total_height canopy_height density weight
0 1 MCD 1.8 1.15 1.30 1.00 1 401.3
1 2 MCD 1.7 1.35 1.35 1.33 1 513.7
2 3 MCD 2.8 2.55 2.16 0.60 1 1179.2
3 4 MCD 1.3 0.85 1.80 1.20 1 308.0
4 5 MCD 3.3 1.90 1.55 1.05 1 855.2

Models

The formulas below are the Python counterparts of the model formulas in the source example. Use lapylace for the Stan-backed Bayesian fit with the same formula interface.

# fit = lp.stan_glm('weight ~ diam1 + diam2 + canopy_height + total_height +', data=mesquite, family=lp.gaussian(), chains=4, iter_sampling=1000)

# fit = lp.stan_glm('log(weight) ~ log(diam1) + log(diam2) + log(canopy_height) +', data=mesquite, family=lp.gaussian(), chains=4, iter_sampling=1000)

# fit = lp.stan_glm('log(weight) ~ log(canopy_volume)', data=mesquite, family=lp.gaussian(), chains=4, iter_sampling=1000)

# fit = lp.stan_glm('log(weight) ~ log(canopy_volume) +', data=mesquite, family=lp.gaussian(), chains=4, iter_sampling=1000)

# fit = lp.stan_glm('log(weight) ~ log(canopy_volume) + log(canopy_shape) +', data=mesquite, family=lp.gaussian(), chains=4, iter_sampling=1000)

Notes

  • Source computation blocks represented: 27.
  • Data paths are expressed through the shared ros_path() helper.
  • Formula-based Bayesian regressions are routed through lapylace.stan_glm().
  • Plotting and simulation work uses NumPy, pandas, matplotlib idioms.