fit_flow_density_with_FN1961_SEP3SigNS5pNuNS3pTauL2p.R

fit_flow_density_with_FN1961_SEP3SigNS5pNuNS3pTauL2p = function(traffic_data, ngrid, upper_density, output_files) {

# Description: This function fits a GAMLSS model to the flow-density values in "traffic_data", and it is designed to be called directly from the R
#              script "FitFun.R". The model component for the functional form of the flow-density relationship is the Franklin-Newell model (FN1961).
#              The model component for the noise in the flow-density relationship is defined as independent observations that follow a Skew
#              Exponential Power Type III distribution. The density dependence of the log of the scale parameter and the log of the skewness
#              parameter is modelled using natural cubic splines with five and three effective free parameters, respectively. The density dependence
#              of the log of the kurtosis parameter is modelled using a straight line function (i.e. an intercept and a gradient as the two free
#              parameters; SEP3SigNS5pNuNS3pTauL2p).
#                The input parameters "ngrid" and "upper_density" are used to define an equally spaced grid of "ngrid" density values ranging from
#              zero to "upper_density". The function employs this density grid to reconstruct the fitted model at the grid points for use in plots
#              and for estimating certain properties of the fitted model that are not directly accessible from the fitted parameter values.
#                The function creates various output files "output_files" including diagnostic plots (see "FitFun.R" for details). If the function
#              finishes successfully, then it returns the corresponding GAMLSS model fit object.
#
# Authors:
#
#   Dan Bramich (dan.bramich@hotmail.co.uk)
#   Lukas Ambuhl (lukas.ambuehl@ivt.baug.ethz.ch)
#
# Configuration Parameters:
#
par1_min = 0.0001      # Minimum acceptable value for the free parameter equivalent to exp(-v_bw*(k_jam/v_ff)) (must be positive and less than unity)
par1_max = 0.9999      # Maximum acceptable value for the free parameter equivalent to exp(-v_bw*(k_jam/v_ff)) (must be positive and less than unity)
par1_step = 0.0001     # Step size for the free parameter equivalent to exp(-v_bw*(k_jam/v_ff)) (must be positive)
par2_min = 0.0001      # Minimum acceptable value for the free parameter k_jam (must be positive)
par2_step = 0.0001     # Step size for the free parameter k_jam (must be positive)
inner_ccrit = 0.05     # Convergence criterion for the inner iteration of the GAMLSS fitting algorithm
inner_ncyc = 10        # Maximum number of cycles of the inner iteration of the GAMLSS fitting algorithm
outer_ccrit = 0.05     # Convergence criterion for the outer iteration of the GAMLSS fitting algorithm
outer_ncyc = 1000      # Maximum number of cycles of the outer iteration of the GAMLSS fitting algorithm


# Define some useful variables
functional_form_model = 'FN1961'
noise_model = 'SEP3SigNS5pNuNS3pTauL2p'

# Report on the GAMLSS model and the data
cat('\n')
cat('>-----------------------------------------------------------------------------<\n')
cat('\n')
cat('The following GAMLSS model will be fit to the flow-density data:\n')
cat('\n')
cat('Model component for the functional form:\n')
cat('  Franklin-Newell (FN1961)\n')
cat('\n')
cat('Model component for the noise:\n')
cat('  Independent observations\n')
cat('  Skew Exponential Power Type III distribution\n')
cat('  Scale is a smooth function of density\n')
cat('  Skewness is a smooth function of density\n')
cat('  Kurtosis is a smooth function of density (SEP3SigNS5pNuNS3pTauL2p)\n')
cat('\n')
cat('Data properties:\n')
tryCatch(
  { ntraffic_data = nrow(traffic_data)
    data_range_density = range(traffic_data$V2)
    data_min_density = data_range_density[1]
    data_max_density = data_range_density[2]
    data_range_flow = range(traffic_data$V3)
    data_min_flow = data_range_flow[1]
    data_max_flow = data_range_flow[2] },
  error = function(cond) { cat('ERROR - Failed to determine the data properties...\n')
                           q(save = 'no', status = 1) }
)
cat('  No. of flow-density measurement pairs (Ndat):', ntraffic_data, '\n')
cat('  Minimum density in the data:                 ', sprintf('%.8g', data_min_density), '\n')
cat('  Maximum density in the data:                 ', sprintf('%.8g', data_max_density), '\n')
cat('  Minimum flow in the data:                    ', sprintf('%.8g', data_min_flow), '\n')
cat('  Maximum flow in the data:                    ', sprintf('%.8g', data_max_flow), '\n')
cat('\n')
cat('Model reconstruction:\n')
tryCatch(
  { grid_density_step = upper_density/(ngrid - 1) },
  error = function(cond) { cat('ERROR - Failed to compute the grid density step...\n')
                           q(save = 'no', status = 1) }
)
cat('  No. of density grid points:', ngrid, '\n')
cat('  Grid lower density:         0\n')
cat('  Grid upper density:        ', sprintf('%.8g', upper_density), '\n')
cat('  Grid density step:         ', sprintf('%.8g', grid_density_step), '\n')

# Fit the GAMLSS model to the data
cat('\n')
cat('Fitting the GAMLSS model...\n')
tryCatch(

  # Fit a Greenshields model to estimate an initial value for v_ff
  { init_model_obj = gamlss(V3 ~ 0 + V2 + I(V2^2), sigma.formula = ~ 1, family = NO(), data = traffic_data)
    if (init_model_obj$converged != TRUE) {
      cat('ERROR - The initial fit of a Greenshields model did not converge...\n')
      q(save = 'no', status = 1)
    }
    v_ff_init = init_model_obj$mu.coefficients[1]

    # Fit a Greenberg model to estimate initial values for k_jam and exp(-v_bw*(k_jam/v_ff))
    init_model_obj = gamlss(V3 ~ 0 + V2 + I(V2*log(V2)), sigma.formula = ~ 1, family = NO(), data = traffic_data)
    if (init_model_obj$converged != TRUE) {
      cat('ERROR - The initial fit of a Greenberg model did not converge...\n')
      q(save = 'no', status = 1)
    }
    par2_init = max(exp(-init_model_obj$mu.coefficients[1]/init_model_obj$mu.coefficients[2]), par2_min + par2_step)
    par1_init = exp(init_model_obj$mu.coefficients[2]*(par2_init/v_ff_init))
    par1_init = max(par1_init, par1_min + par1_step)
    par1_init = min(par1_init, par1_max - par1_step)

    # Perform the intermediate fits
    inner_ccrit_use = data.frame(inner_ccrit_use = inner_ccrit)
    inner_ncyc_use = data.frame(inner_ncyc_use = inner_ncyc)
    outer_ccrit_use = data.frame(outer_ccrit_use = outer_ccrit)
    outer_ncyc_use = data.frame(outer_ncyc_use = outer_ncyc)
    model_formula = quote(gamlss(V3 ~ 0 + I(V2*(1.0 - (p[1]^((1.0/V2) - (1.0/p[2]))))), sigma.formula = ~ ns(V2, df = 4), nu.formula = ~ ns(V2, df = 2),
                                 tau.formula = ~ 1 + V2, family = SEP3(), c.crit = outer_ccrit_use, n.cyc = outer_ncyc_use,
                                 i.control = glim.control(cc = inner_ccrit_use, cyc = inner_ncyc_use)))
    attach(traffic_data)
    attach(inner_ccrit_use)
    attach(inner_ncyc_use)
    attach(outer_ccrit_use)
    attach(outer_ncyc_use)
    optim_obj = try(find.hyper(model = model_formula, parameters = c(par1_init, par2_init), k = 0.0, steps = c(par1_step, par2_step), lower = c(par1_min, par2_min),
                               upper = c(par1_max, Inf), maxit = 500))
    if (class(optim_obj) == 'try-error') { optim_obj = list(convergence = 1) }
    if (optim_obj$convergence != 0) {
      par1_min_use = data.frame(par1_min_use = par1_min)
      par1_max_use = data.frame(par1_max_use = par1_max)
      par2_min_use = data.frame(par2_min_use = par2_min)
      model_formula = quote(gamlss(V3 ~ 0 + I(V2*(1.0 - (max(min(p[1], par1_max_use), par1_min_use)^((1.0/V2) - (1.0/(par2_min_use + abs(p[2]))))))),
                                   sigma.formula = ~ ns(V2, df = 4), nu.formula = ~ ns(V2, df = 2), tau.formula = ~ 1 + V2, family = SEP3(),
                                   c.crit = outer_ccrit_use, n.cyc = outer_ncyc_use, i.control = glim.control(cc = inner_ccrit_use, cyc = inner_ncyc_use)))
      attach(par1_min_use)
      attach(par1_max_use)
      attach(par2_min_use)
      optim_obj = try(find.hyper(model = model_formula, parameters = c(par1_init, par2_init - par2_min), k = 0.0, steps = c(par1_step, par2_step),
                                 method = 'Nelder-Mead', maxit = 500))
      if (class(optim_obj) == 'try-error') { optim_obj = list(convergence = 1) }
      detach(par2_min_use)
      detach(par1_max_use)
      detach(par1_min_use)
      detach(outer_ncyc_use)
      detach(outer_ccrit_use)
      detach(inner_ncyc_use)
      detach(inner_ccrit_use)
      detach(traffic_data)
      if (optim_obj$convergence != 0) {
        cat('ERROR - The intermediate fits did not converge...\n')
        q(save = 'no', status = 1)
      }
      par1 = max(min(optim_obj$par[1], par1_max), par1_min)
      par2 = par2_min + abs(optim_obj$par[2])
    } else {
      detach(outer_ncyc_use)
      detach(outer_ccrit_use)
      detach(inner_ncyc_use)
      detach(inner_ccrit_use)
      detach(traffic_data)
      par1 = optim_obj$par[1]
      par2 = optim_obj$par[2]
    }

    # Perform the final fit
    model_obj = gamlss(V3 ~ 0 + I(V2*(1.0 - (par1^((1.0/V2) - (1.0/par2))))), sigma.formula = ~ ns(V2, df = 4), nu.formula = ~ ns(V2, df = 2),
                       tau.formula = ~ 1 + V2, family = SEP3(), data = traffic_data, c.crit = outer_ccrit, n.cyc = outer_ncyc,
                       i.control = glim.control(cc = inner_ccrit, cyc = inner_ncyc))
    if (model_obj$converged != TRUE) {
      cat('ERROR - The final fit did not converge...\n')
      q(save = 'no', status = 1)
    }
    model_obj$mu.df = model_obj$mu.df + 2
    model_obj$df.fit = model_obj$df.fit + 2
    model_obj$df.residual = model_obj$df.residual - 2
    model_obj$aic = model_obj$aic + 4.0
    model_obj$sbc = model_obj$sbc + 2.0*log(ntraffic_data) },
  error = function(cond) { cat('ERROR - Failed to fit the GAMLSS model...\n')
                           q(save = 'no', status = 1) }
)

# Store the predicted values for the fitted model at the density values in the data in the data table
cat('\n')
cat('Storing the predicted values for the fitted model in the data table...\n')
tryCatch(
  { if (!all(is.finite(model_obj$mu.fv))) {
      cat('ERROR - The predicted values for "mu" at the density values in the data include at least one value that is infinite...\n')
      q(save = 'no', status = 1)
    }
    if (!all(is.finite(model_obj$sigma.fv))) {
      cat('ERROR - The predicted values for "sigma" at the density values in the data include at least one value that is infinite...\n')
      q(save = 'no', status = 1)
    }
    if (any(model_obj$sigma.fv <= 0.0)) {
      cat('ERROR - The predicted values for "sigma" at the density values in the data include at least one value that is zero or negative...\n')
      q(save = 'no', status = 1)
    }
    if (!all(is.finite(model_obj$nu.fv))) {
      cat('ERROR - The predicted values for "nu" at the density values in the data include at least one value that is infinite...\n')
      q(save = 'no', status = 1)
    }
    if (any(model_obj$nu.fv <= 0.0)) {
      cat('ERROR - The predicted values for "nu" at the density values in the data include at least one value that is zero or negative...\n')
      q(save = 'no', status = 1)
    }
    if (!all(is.finite(model_obj$tau.fv))) {
      cat('ERROR - The predicted values for "tau" at the density values in the data include at least one value that is infinite...\n')
      q(save = 'no', status = 1)
    }
    if (any(model_obj$tau.fv <= 0.0)) {
      cat('ERROR - The predicted values for "tau" at the density values in the data include at least one value that is zero or negative...\n')
      q(save = 'no', status = 1)
    }
    traffic_data[, fitted_values_mu := model_obj$mu.fv]
    traffic_data[, fitted_values_sigma := model_obj$sigma.fv]
    traffic_data[, fitted_values_nu := model_obj$nu.fv]
    traffic_data[, fitted_values_tau := model_obj$tau.fv] },
  error = function(cond) { cat('ERROR - Failed to store the predicted values for the fitted model...\n')
                           q(save = 'no', status = 1) }
)

# Compute the normalised quantile residuals and store them in the data table. Note that the normalised quantile residuals may include some "-Inf"
# or "Inf" values for particularly bad outliers.
cat('Computing the normalised quantile residuals...\n')
tryCatch(
  { cumulative_probs_lower = pSEP3(traffic_data$V3, mu = model_obj$mu.fv, sigma = model_obj$sigma.fv, nu = model_obj$nu.fv, tau = model_obj$tau.fv)
    cumulative_probs_upper = pSEP3(traffic_data$V3, mu = model_obj$mu.fv, sigma = model_obj$sigma.fv, nu = model_obj$nu.fv, tau = model_obj$tau.fv, lower.tail = FALSE)
    traffic_data[, normalised_quantile_residuals := calculate_normalised_quantile_residuals(cumulative_probs_lower, cumulative_probs_upper)] },
  error = function(cond) { cat('ERROR - Failed to compute the normalised quantile residuals...\n')
                           q(save = 'no', status = 1) }
)

# Compute the percentiles for the fitted model at the density values in the data and store them in the data table
cat('Computing the percentiles for the fitted model at the density values in the data...\n')
tryCatch(
  { traffic_data[, percentile_m3sig := qSEP3(pNO(-3.0), mu = traffic_data$fitted_values_mu, sigma = traffic_data$fitted_values_sigma,
                                             nu = traffic_data$fitted_values_nu, tau = traffic_data$fitted_values_tau)]
    traffic_data[, percentile_m2sig := qSEP3(pNO(-2.0), mu = traffic_data$fitted_values_mu, sigma = traffic_data$fitted_values_sigma,
                                             nu = traffic_data$fitted_values_nu, tau = traffic_data$fitted_values_tau)]
    traffic_data[, percentile_m1sig := qSEP3(pNO(-1.0), mu = traffic_data$fitted_values_mu, sigma = traffic_data$fitted_values_sigma,
                                             nu = traffic_data$fitted_values_nu, tau = traffic_data$fitted_values_tau)]
    traffic_data[, percentile_0sig := qSEP3(0.5, mu = traffic_data$fitted_values_mu, sigma = traffic_data$fitted_values_sigma,
                                            nu = traffic_data$fitted_values_nu, tau = traffic_data$fitted_values_tau)]
    traffic_data[, percentile_p1sig := qSEP3(pNO(1.0), mu = traffic_data$fitted_values_mu, sigma = traffic_data$fitted_values_sigma,
                                             nu = traffic_data$fitted_values_nu, tau = traffic_data$fitted_values_tau)]
    traffic_data[, percentile_p2sig := qSEP3(pNO(2.0), mu = traffic_data$fitted_values_mu, sigma = traffic_data$fitted_values_sigma,
                                             nu = traffic_data$fitted_values_nu, tau = traffic_data$fitted_values_tau)]
    traffic_data[, percentile_p3sig := qSEP3(pNO(3.0), mu = traffic_data$fitted_values_mu, sigma = traffic_data$fitted_values_sigma,
                                             nu = traffic_data$fitted_values_nu, tau = traffic_data$fitted_values_tau)] },
  error = function(cond) { cat('ERROR - Failed to compute the percentiles for the fitted model at the density values in the data...\n')
                           q(save = 'no', status = 1) }
)

# Compute a set of distributional measures for the fitted model at the density values in the data and store them in the data table
cat('Computing a set of distributional measures for the fitted model at the density values in the data...\n')
tryCatch(
  { distributional_measures = calculate_distributional_measures_for_SEP3(traffic_data$fitted_values_mu, traffic_data$fitted_values_sigma,
                                                                         traffic_data$fitted_values_nu, traffic_data$fitted_values_tau)
    traffic_data[, mean := distributional_measures$mean]
    traffic_data[, median := double(length = ntraffic_data)]
    traffic_data[, mode := distributional_measures$mode]
    traffic_data[, standard_deviation := distributional_measures$standard_deviation]
    traffic_data[, moment_skewness := distributional_measures$moment_skewness]
    traffic_data[, moment_excess_kurtosis := distributional_measures$moment_excess_kurtosis] },
  error = function(cond) { cat('ERROR - Failed to compute a set of distributional measures for the fitted model at the density values in the data...\n')
                           q(save = 'no', status = 1) }
)

# Reconstruct the fitted model over the density range from zero to "upper_density"
cat('Reconstructing the fitted model over the density range from 0 to', sprintf('%.8g', upper_density), '...\n')
tryCatch(
  { reconstructed_model_fit = data.table(V2 = seq(from = 0.0, to = upper_density, length.out = ngrid))
    predicted_values_for_mu = double(length = ngrid)
    predicted_values_for_mu[2:ngrid] = predict(model_obj, what = 'mu', newdata = reconstructed_model_fit[2:ngrid], type = 'response', data = traffic_data)
    predicted_values_for_sigma = predict(model_obj, what = 'sigma', newdata = reconstructed_model_fit, type = 'response', data = traffic_data)
    predicted_values_for_nu = predict(model_obj, what = 'nu', newdata = reconstructed_model_fit, type = 'response', data = traffic_data)
    predicted_values_for_tau = predict(model_obj, what = 'tau', newdata = reconstructed_model_fit, type = 'response', data = traffic_data)
    predicted_values_for_sigma = fix_out_of_data_curves(reconstructed_model_fit$V2, predicted_values_for_sigma, data_min_density, data_max_density)
    predicted_values_for_nu = fix_out_of_data_curves(reconstructed_model_fit$V2, predicted_values_for_nu, data_min_density, data_max_density)
    predicted_values_for_tau = fix_out_of_data_curves(reconstructed_model_fit$V2, predicted_values_for_tau, data_min_density, data_max_density)
    if (!all(is.finite(predicted_values_for_mu))) {
      cat('ERROR - The reconstructed fitted model for "mu" includes at least one value that is infinite...\n')
      q(save = 'no', status = 1)
    }
    if (!all(is.finite(predicted_values_for_sigma))) {
      cat('ERROR - The reconstructed fitted model for "sigma" includes at least one value that is infinite...\n')
      q(save = 'no', status = 1)
    }
    if (any(predicted_values_for_sigma <= 0.0)) {
      cat('ERROR - The reconstructed fitted model for "sigma" includes at least one value that is zero or negative...\n')
      q(save = 'no', status = 1)
    }
    if (!all(is.finite(predicted_values_for_nu))) {
      cat('ERROR - The reconstructed fitted model for "nu" includes at least one value that is infinite...\n')
      q(save = 'no', status = 1)
    }
    if (any(predicted_values_for_nu <= 0.0)) {
      cat('ERROR - The reconstructed fitted model for "nu" includes at least one value that is zero or negative...\n')
      q(save = 'no', status = 1)
    }
    if (!all(is.finite(predicted_values_for_tau))) {
      cat('ERROR - The reconstructed fitted model for "tau" includes at least one value that is infinite...\n')
      q(save = 'no', status = 1)
    }
    if (any(predicted_values_for_tau <= 0.0)) {
      cat('ERROR - The reconstructed fitted model for "tau" includes at least one value that is zero or negative...\n')
      q(save = 'no', status = 1)
    }
    reconstructed_model_fit[, mu := predicted_values_for_mu]
    reconstructed_model_fit[, sigma := predicted_values_for_sigma]
    reconstructed_model_fit[, nu := predicted_values_for_nu]
    reconstructed_model_fit[, tau := predicted_values_for_tau] },
  error = function(cond) { cat('ERROR - Failed to reconstruct the fitted model over the required density range...\n')
                           q(save = 'no', status = 1) }
)

# Construct percentile curves for the fitted model over the density range from zero to "upper_density"
cat('Constructing percentile curves for the fitted model over the density range from 0 to', sprintf('%.8g', upper_density), '...\n')
tryCatch(
  { reconstructed_model_fit[, percentile_m3sig := qSEP3(pNO(-3.0), mu = reconstructed_model_fit$mu, sigma = reconstructed_model_fit$sigma,
                                                        nu = reconstructed_model_fit$nu, tau = reconstructed_model_fit$tau)]
    reconstructed_model_fit[, percentile_m2sig := qSEP3(pNO(-2.0), mu = reconstructed_model_fit$mu, sigma = reconstructed_model_fit$sigma,
                                                        nu = reconstructed_model_fit$nu, tau = reconstructed_model_fit$tau)]
    reconstructed_model_fit[, percentile_m1sig := qSEP3(pNO(-1.0), mu = reconstructed_model_fit$mu, sigma = reconstructed_model_fit$sigma,
                                                        nu = reconstructed_model_fit$nu, tau = reconstructed_model_fit$tau)]
    reconstructed_model_fit[, percentile_0sig := qSEP3(0.5, mu = reconstructed_model_fit$mu, sigma = reconstructed_model_fit$sigma,
                                                       nu = reconstructed_model_fit$nu, tau = reconstructed_model_fit$tau)]
    reconstructed_model_fit[, percentile_p1sig := qSEP3(pNO(1.0), mu = reconstructed_model_fit$mu, sigma = reconstructed_model_fit$sigma,
                                                        nu = reconstructed_model_fit$nu, tau = reconstructed_model_fit$tau)]
    reconstructed_model_fit[, percentile_p2sig := qSEP3(pNO(2.0), mu = reconstructed_model_fit$mu, sigma = reconstructed_model_fit$sigma,
                                                        nu = reconstructed_model_fit$nu, tau = reconstructed_model_fit$tau)]
    reconstructed_model_fit[, percentile_p3sig := qSEP3(pNO(3.0), mu = reconstructed_model_fit$mu, sigma = reconstructed_model_fit$sigma,
                                                        nu = reconstructed_model_fit$nu, tau = reconstructed_model_fit$tau)] },
  error = function(cond) { cat('ERROR - Failed to construct percentile curves for the fitted model over the required density range...\n')
                           q(save = 'no', status = 1) }
)

# Construct curves of a set of distributional measures for the fitted model over the density range from zero to "upper_density"
cat('Constructing curves of a set of distributional measures for the fitted model over the density range from 0 to', sprintf('%.8g', upper_density), '...\n')
tryCatch(
  { distributional_measures = calculate_distributional_measures_for_SEP3(reconstructed_model_fit$mu, reconstructed_model_fit$sigma, reconstructed_model_fit$nu, reconstructed_model_fit$tau)
    reconstructed_model_fit[, mean := distributional_measures$mean]
    reconstructed_model_fit[, median := double(length = ngrid)]
    reconstructed_model_fit[, mode := distributional_measures$mode]
    reconstructed_model_fit[, standard_deviation := distributional_measures$standard_deviation]
    reconstructed_model_fit[, moment_skewness := distributional_measures$moment_skewness]
    reconstructed_model_fit[, moment_excess_kurtosis := distributional_measures$moment_excess_kurtosis] },
  error = function(cond) { cat('ERROR - Failed to construct curves of a set of distributional measures for the fitted model over the required density range...\n')
                           q(save = 'no', status = 1) }
)

# Estimate useful properties of the fitted model over the density range from zero to the maximum observed density using the reconstruction
cat('Estimating useful properties of the fitted model using the reconstruction...\n')
tryCatch(
  { selection = reconstructed_model_fit$V2 < (data_max_density + grid_density_step)
    reconstructed_model_fit_selection = reconstructed_model_fit[selection]
    curve_properties_for_mu_over_data_range = get_curve_properties_for_mu(reconstructed_model_fit_selection, 'Flow.Density')
    curve_properties_for_sigma_over_data_range = get_curve_properties_for_sigma(reconstructed_model_fit_selection, curve_properties_for_mu_over_data_range)
    curve_properties_for_nu_over_data_range = get_curve_properties_for_nu(reconstructed_model_fit_selection, curve_properties_for_mu_over_data_range)
    curve_properties_for_tau_over_data_range = get_curve_properties_for_tau(reconstructed_model_fit_selection, curve_properties_for_mu_over_data_range) },
  error = function(cond) { cat('ERROR - Failed to estimate useful properties of the fitted model over the density range from zero to the maximum observed density...\n')
                           q(save = 'no', status = 1) }
)

# Estimate useful properties of the fitted model over the density range from zero to "upper_density" using the reconstruction
tryCatch(
  { curve_properties_for_mu_over_full_range = get_curve_properties_for_mu(reconstructed_model_fit, 'Flow.Density')
    curve_properties_for_sigma_over_full_range = get_curve_properties_for_sigma(reconstructed_model_fit, curve_properties_for_mu_over_full_range)
    curve_properties_for_nu_over_full_range = get_curve_properties_for_nu(reconstructed_model_fit, curve_properties_for_mu_over_full_range)
    curve_properties_for_tau_over_full_range = get_curve_properties_for_tau(reconstructed_model_fit, curve_properties_for_mu_over_full_range) },
  error = function(cond) { cat('ERROR - Failed to estimate useful properties of the fitted model over the density range from zero to "upper_density"...\n')
                           q(save = 'no', status = 1) }
)

# Extract information from the model fit object for the fit summary
cat('Extracting information from the model fit object for the fit summary...\n')
tryCatch(
  { npar_mu = model_obj$mu.df
    npar_sigma = model_obj$sigma.df
    npar_nu = model_obj$nu.df
    npar_tau = model_obj$tau.df
    npar_all = model_obj$df.fit
    gdev = model_obj$G.deviance
    aic = model_obj$aic
    bic = model_obj$sbc },
  error = function(cond) { cat('ERROR - Failed to extract information from the model fit object for the fit summary...\n')
                           q(save = 'no', status = 1) }
)

# Where possible, extract physical parameter values from the model fit object for the fit summary
tryCatch(
  { q_0 = 0.0
    v_ff = NA
    dvdk_0 = 0.0
    k_crit = NA
    k_vmax = NA
    q_cap = NA
    v_max = NA
    k_jam = NA
    v_bw = NA
    dvdk_kjam = NA
    if (model_obj$mu.coefficients[1] > 0.0) {
      v_ff = model_obj$mu.coefficients[1]
#      k_crit =     # Can be computed numerically - not yet implemented
      k_vmax = 0.0
#      q_cap =      # Can be computed once k_crit is available - not yet implemented
      v_max = v_ff
      k_jam = par2
      v_bw = (-v_ff/k_jam)*log(par1)
      dvdk_kjam = -v_bw/k_jam
    } },
  error = function(cond) { cat('ERROR - Failed to extract physical parameter values from the model fit object for the fit summary...\n')
                           q(save = 'no', status = 1) }
)

# Report a basic fit summary
cat('\n')
cat('Fit summary (abridged):\n')
cat('\n')
cat('Model parameter counts:\n')
cat('  No. of free parameters (mu):        ', sprintf('%.6f', npar_mu), '\n')
cat('  No. of free parameters (sigma):     ', sprintf('%.6f', npar_sigma), '\n')
cat('  No. of free parameters (nu):        ', sprintf('%.6f', npar_nu), '\n')
cat('  No. of free parameters (tau):       ', sprintf('%.6f', npar_tau), '\n')
cat('  Total no. of free parameters (Npar):', sprintf('%.6f', npar_all), '\n')
cat('\n')
cat('Fit quality:\n')
cat('  Global deviance [ -2*ln(Lmax) ]:    ', sprintf('%.4f', gdev), '\n')
cat('  AIC [ -2*ln(Lmax) + 2*Npar ]:       ', sprintf('%.4f', aic), '\n')
cat('  BIC [ -2*ln(Lmax) + Npar*ln(Ndat) ]:', sprintf('%.4f', bic), '\n')
cat('\n')
cat('Fitted physical parameters (where available):\n')
if (!is.na(q_0)) { cat('  Flow at zero density:                                  ', sprintf('%.8g', q_0), '\n') }
if (!is.na(v_ff)) { cat('  Free-flow speed:                                       ', sprintf('%.8g', v_ff), '\n') }
if (!is.na(dvdk_0)) { cat('  Gradient of the speed (w.r.t. density) at zero density:', sprintf('%.8g', dvdk_0), '\n') }
if (!is.na(k_crit)) { cat('  Critical density:                                      ', sprintf('%.8g', k_crit), '\n') }
if (!is.na(k_vmax)) { cat('  Density at maximum speed:                              ', sprintf('%.8g', k_vmax), '\n') }
if (!is.na(q_cap)) { cat('  Capacity:                                              ', sprintf('%.8g', q_cap), '\n') }
if (!is.na(v_max)) { cat('  Maximum speed:                                         ', sprintf('%.8g', v_max), '\n') }
if (!is.na(k_jam)) { cat('  Jam density:                                           ', sprintf('%.8g', k_jam), '\n') }
if (!is.na(v_bw)) { cat('  Back-propagating wave speed at jam density:            ', sprintf('%.8g', v_bw), '\n') }
if (!is.na(dvdk_kjam)) { cat('  Gradient of the speed (w.r.t. density) at jam density: ', sprintf('%.8g', dvdk_kjam), '\n') }
cat('\n')
cat('Fitted model parameters (see the accompanying papers by Bramich, Menendez & Ambuhl for details):\n')
cat('  v_ff:                   ', sprintf('%.8g', model_obj$mu.coefficients[1]), '\n')
cat('  exp(-v_bw*(k_jam/v_ff)):', sprintf('%.8g', par1), '\n')
cat('  k_jam:                  ', sprintf('%.8g', par2), '\n')

# Write out the fit summary file "Fit.Summary.<fd_type>.<functional_form_model>.<noise_model>.txt"
cat('\n')
cat('Writing out the fit summary file:    ', output_files[1], '\n')
tryCatch(
  { write_fit_summary(output_files[1], 'Flow.Density', ntraffic_data, data_min_density, data_max_density, data_min_flow, data_max_flow,
                      npar_mu, npar_sigma, npar_nu, npar_tau, npar_all, gdev, aic, bic,
                      q_0, v_ff, dvdk_0, k_crit, k_vmax, q_cap, v_max, k_jam, v_bw, dvdk_kjam,
                      curve_properties_for_mu_over_data_range, curve_properties_for_sigma_over_data_range,
                      curve_properties_for_nu_over_data_range, curve_properties_for_tau_over_data_range,
                      curve_properties_for_mu_over_full_range, curve_properties_for_sigma_over_full_range,
                      curve_properties_for_nu_over_full_range, curve_properties_for_tau_over_full_range)
    cat('######################################################################################################################\n',
        '# FITTED MODEL PARAMETERS (SEE THE ACCOMPANYING PAPERS BY BRAMICH, MENENDEZ & AMBUHL FOR DETAILS)\n',
        '# N.B: FITTED COEFFICIENTS FOR ANY NON-PARAMETRIC SMOOTHING FUNCTIONS IN THE MODEL ARE NOT REPORTED HERE\n',
        '######################################################################################################################\n',
        sprintf('%.8g', model_obj$mu.coefficients[1]), '           # v_ff\n',
        sprintf('%.8g', par1), '           # exp(-v_bw*(k_jam/v_ff))\n',
        sprintf('%.8g', par2), '           # k_jam\n',
        file = output_files[1], sep = '', append = TRUE) },
  error = function(cond) { cat('ERROR - Failed to write out the fit summary file...\n')
                           remove_file_list(output_files)
                           q(save = 'no', status = 1) }
)

# Write out the fit curves file "Fit.Curves.<fd_type>.<functional_form_model>.<noise_model>.txt"
cat('Writing out the fit curves file:     ', output_files[2], '\n')
tryCatch(
  { cat('# Density : Mu : Sigma : Nu : Tau : 0.135 Percentile (Corresponding To -3*Sigma In A Normal Distribution) : 2.28 Percentile (Corresponding To -2*Sigma In A',
        'Normal Distribution) : 15.87 Percentile (Corresponding To -1*Sigma In A Normal Distribution) : 50.00 Percentile (Median) : 84.13 Percentile (Corresponding',
        'To 1*Sigma In A Normal Distribution) : 97.72 Percentile (Corresponding To 2*Sigma In A Normal Distribution) : 99.865 Percentile (Corresponding To 3*Sigma',
        'In A Normal Distribution) : Mean : Median : Mode : Standard Deviation : Moment Skewness : Moment Excess Kurtosis\n', file = output_files[2])
    write.table(reconstructed_model_fit, file = output_files[2], append = TRUE, quote = FALSE, row.names = FALSE, col.names = FALSE) },
  error = function(cond) { cat('ERROR - Failed to write out the fit curves file...\n')
                           remove_file_list(output_files)
                           q(save = 'no', status = 1) }
)

# Write out the fit predictions file "Fit.Predictions.<fd_type>.<functional_form_model>.<noise_model>.txt"
cat('Writing out the fit predictions file:', output_files[3], '\n')
tryCatch(
  { cat('# Data Column 1 : Data Column 2 : Data Column 3 : Fitted Value For Mu : Fitted Value For Sigma : Fitted Value For Nu : Fitted Value For Tau :',
        'Normalised Quantile Residual ("-Inf" Or "Inf" Values May Be Present) : 0.135 Percentile (Corresponding To -3*Sigma In A Normal Distribution) :',
        '2.28 Percentile (Corresponding To -2*Sigma In A Normal Distribution) : 15.87 Percentile (Corresponding To -1*Sigma In A Normal Distribution) :',
        '50.00 Percentile (Median) : 84.13 Percentile (Corresponding To 1*Sigma In A Normal Distribution) : 97.72 Percentile (Corresponding To 2*Sigma',
        'In A Normal Distribution) : 99.865 Percentile (Corresponding To 3*Sigma In A Normal Distribution) : Mean : Median : Mode : Standard Deviation :',
        'Moment Skewness : Moment Excess Kurtosis\n', file = output_files[3])
    write.table(traffic_data, file = output_files[3], append = TRUE, quote = FALSE, row.names = FALSE, col.names = FALSE) },
  error = function(cond) { cat('ERROR - Failed to write out the fit predictions file...\n')
                           remove_file_list(output_files)
                           q(save = 'no', status = 1) }
)

# If required, then create the plots for the GAMLSS model fit
if (length(output_files) > 3) {
  cat('\n')
  cat('Creating the plots for the GAMLSS model fit...\n')
  tryCatch(
    { create_all_plots(traffic_data, ntraffic_data, data_max_density, upper_density, reconstructed_model_fit_selection, reconstructed_model_fit, ngrid,
                       'Flow.Density', functional_form_model, noise_model, output_files) },
    error = function(cond) { cat('ERROR - Failed to create the plot...\n')
                             remove_file_list(output_files)
                             q(save = 'no', status = 1) }
  )
}

# Return the GAMLSS model fit object
cat('\n')
cat('>-----------------------------------------------------------------------------<\n')
return(model_obj)
}