Adsorption — Fit and Inspect

MODEL "Adsorption — Fit and Inspect"
DESCRIPTION "Fit NIST Misra1a observations and inspect coefficients, residuals, stopping reason, and local identifiability."
VERSION "1.0.0"
AUTHOR "Grid Team"
TAGS "chemistry", "adsorption", "fitting", "diagnostics"

# Observed NIST dental monomolecular adsorption data, Misra (1978).
# Source: https://www.itl.nist.gov/div898/strd/nls/data/LINKS/DATA/Misra1a.dat
# NIST names pressure and volume without specifying their measurement units.
PressureData = {77.6,114.9,141.1,190.8,239.9,289,332.8,378.4,434.8,477.3,536.8,593.1,689.1,760}
VolumeData = {10.07,14.73,17.94,23.93,29.61,35.18,40.02,44.82,50.76,55.05,61.01,66.40,75.47,81.78}

# These bounds guide initial guesses; they are not constraints on the fit.
input InitialCapacity = 250
  VALIDATE BETWEEN 1 AND 1000 MESSAGE "Choose an initial capacity from 1 to 1000 for this example."
input InitialScaledRate = 0.5
  VALIDATE BETWEEN 0.001 AND 5 MESSAGE "Choose an initial scaled rate from 0.001 to 5."
input IterationBudget = 1000
  VALIDATE BETWEEN 1 AND 1000 MESSAGE "Choose an integer iteration budget from 1 to 1000."

# Scaling the rate makes the parameter units explicit and improves scaling.
# The NIST coefficient b2 equals scaled_rate / 1000.
output Fit = CURVE_FIT.REPORT(
  (pressure, capacity, scaled_rate) => capacity*(1-EXP(-scaled_rate*pressure/1000)),
  PressureData, VolumeData, {InitialCapacity, InitialScaledRate}, 1e-12, IterationBudget
)
FitSucceeded = NOT ISERROR(Fit)
output Capacity = FitSucceeded THEN INDEX(Fit.parameters,1,1) ELSE #N/A
output RatePerPressure = FitSucceeded THEN INDEX(Fit.parameters,1,2)/1000 ELSE #N/A
output SumSquaredResiduals = FitSucceeded THEN Fit.ssr ELSE #N/A
output Interpretation = FitSucceeded THEN "Inspect Fit for the stopping reason, rank, conditioning, and residuals. A returned fit is not proof of model adequacy." ELSE "Fit failed. Inspect Fit for the native error and recovery advice."
output Assumptions = "Unweighted, unconstrained least squares. Residuals are fitted minus observed. No parameter confidence intervals are estimated. Scaling affects the reported Jacobian condition number."
output ScenarioGuide = "Save the initial state. Compare InitialCapacity 500 and InitialScaledRate 0.1 with the starting values 250 and 0.5. Set IterationBudget to 1 to inspect an explicit failure; restore 1000 to recover."

FitData = FitSucceeded THEN MAP(SEQUENCE(1,14), (observation) => {
  pressure: INDEX(PressureData,1,observation),
  measured: INDEX(VolumeData,1,observation),
  fitted: INDEX(Fit.fitted_values,1,observation),
  residual: INDEX(Fit.residuals,1,observation)
}) ELSE []

A1 = "Fitted capacity (NIST volume scale)"
B1 = Capacity
A2 = "Fitted b2 (per NIST pressure scale)"
B2 = RatePerPressure
A3 = "Sum of squared residuals"
B3 = SumSquaredResiduals
A5 = Interpretation
A6 = Assumptions
A7 = ScenarioGuide

Charts_1!type = "charts"
Charts_1!config = <toml>
  version = 1
  kind = "charts"
  title = "NIST adsorption fit"
  [layout]
  x_axis = "Pressure (NIST scale)"
  y_axis = "Volume (NIST scale)"
  [[charts]]
  type = "scatter"
  bind = "FitData"
  x = "pressure"
  y = "measured"
  label = "Observed"
  [[charts]]
  type = "line"
  smooth = false
  bind = "FitData"
  x = "pressure"
  y = "fitted"
  label = "Fitted at observation pressures"
</toml>

Charts_2!type = "charts"
Charts_2!config = <toml>
  version = 1
  kind = "charts"
  title = "Adsorption residuals"
  [layout]
  x_axis = "Pressure (NIST scale)"
  y_axis = "Fitted minus observed volume"
  [[charts]]
  type = "scatter"
  bind = "FitData"
  x = "pressure"
  y = "residual"
  label = "Residual"
</toml>

END MODEL