Multi-Regimes Fundamental Diagram – New Calibration Approach ¶
Highlights ¶
For almost a century, several models have been developed to calibrate the pairwise relationship among traffic flow variables, i.e., speed, density, and flow.
Multi-regime models are well known for being superior in fitting the speed-density relationship over single-regime models.
However, in modeling multi-regimes models, breakpoints separating the regimes are generally estimated by a model developer based on the subjective judgment of visual data characteristics.
A new calibration approach based on Bayesian inference is introduced for fitting the multi-regimes models.
The superiority of the new calibration approach is the ability to locate the breakpoints in the multi-regime considering the specific data characteristics.
More information: ¶
Python Code ¶
If you want to replicate tha analysis feel free to use the code.
In [1]:
import pandas as pd, matplotlib.pyplot as plt, seaborn as sns
%pylab inline
%config InlineBackend.figure_format = 'retina'
In [2]:
from rethinking import MAP, coef, extract_samples, link, precis, sim, vcov
In [3]:
# Import the "Default" dataset.
data = pd.read_csv("2262-3.csv")[['occupancy','speed']]
data.head()
Out[3]:
In [4]:
import time
import torch
import numpy as onp
import jax.numpy as np
from jax import jit, random
# NB: replace gpu by cpu to run this notebook in cpu
from jax.config import config; config.update("jax_platform_name", "cpu")
import numpyro.distributions as dist
from numpyro.diagnostics import summary
import arviz as az
from numpyro.handlers import sample
from numpyro.hmc_util import initialize_model
from numpyro.mcmc import hmc, mcmc
from numpyro.util import fori_collect
In [5]:
%config InlineBackend.figure_format = 'retina'
In [6]:
occupancy = data.occupancy
speed = data.speed
In [7]:
f, ax = plt.subplots(1, 1, figsize=(10, 5))
plt.scatter(occupancy,speed, color='g', s=50)
Out[7]:
In [8]:
def sigmoid_f(x):
return 1 / (1 + np.exp(-x))
def weight_f(w, mu1, mu2):
mu = (1.0 - w) * mu1 + w * mu2
return mu
def weight_three(w1, w2, mu1, mu2, mu3):
out_ = (1.0 - w1) * mu1 + w1 * mu2
mu = (1.0 - w2) * out_ + w2 * mu3
return mu
from jax.scipy.special import erf
def cumulative_normal(x):
'Cumulative distribution function for the standard normal distribution'
return 0.5 + 0.5 * erf(x / math.sqrt(2.0))
Two-regime model¶
In [9]:
def model_twoRegime(occupancy, speed):
α1 = sample('α1', dist.Normal(0, 10))
α2 = sample('α2', dist.Normal(0, 10))
β1 = sample('β1', dist.Normal(0, 10))
β2 = sample('β2', dist.Normal(0, 10))
σ1 = sample('σ1', dist.HalfCauchy(0.5)) #Noise
σ2 = sample('σ2', dist.HalfCauchy(0.5))
switch_point = sample("switch_point", dist.Uniform(occupancy.min(), occupancy.max()))
mu1 = α1 + β1 * occupancy
mu2 = α2 + β2 * occupancy
w = cumulative_normal(occupancy - switch_point)
mu = weight_f(w, mu1, mu2)
sigma = weight_f(w, σ1, σ2)
return sample("obs", dist.Normal(mu, sigma), obs = speed)
In [10]:
'''`hmc` is implemented (default).'''
init_params, potential_fn, constrain_fn = initialize_model(random.PRNGKey(1), model_twoRegime,
occupancy.values, speed.values)
num_warmup, num_samples = 1000, 1000
samples_twoRegime = mcmc(num_warmup, num_samples, init_params,
potential_fn=potential_fn,
constrain_fn=constrain_fn,
print_summary=False)
In [11]:
summary(samples_twoRegime, prob=0.95)
Three-regime model¶
In [12]:
def model_threeRegime(occupancy, speed):
α1 = sample('α1', dist.Normal(0, 10))
α2 = sample('α2', dist.Normal(0, 10))
α3 = sample('α3', dist.Normal(0, 10))
β1 = sample('β1', dist.Normal(0, 10))
β2 = sample('β2', dist.Normal(0, 10))
β3 = sample('β3', dist.Normal(0, 10))
σ1 = sample('σ1', dist.HalfCauchy(0.5)) #Noise
σ2 = sample('σ2', dist.HalfCauchy(0.5))
σ3 = sample('σ3', dist.HalfCauchy(0.5))
switch_point1 = sample("switch_point1", dist.Uniform(occupancy.min(), occupancy.max()))
switch_point2 = sample("switch_point2", dist.Uniform(17, occupancy.max()))
mu1 = α1 + β1 * occupancy
mu2 = α2 + β2 * occupancy
mu3 = α3 + β3 * occupancy
w1 =cumulative_normal(occupancy - switch_point1)
w2 = cumulative_normal(occupancy - switch_point2)
mu = weight_three(w1, w2, mu1, mu2, mu3)
sigma = weight_three(w1, w2, σ1, σ2, σ3)
return sample("obs", dist.Normal(mu, sigma), obs = speed)
In [13]:
'''`hmc` is implemented (default).'''
init_params, potential_fn, constrain_fn = initialize_model(random.PRNGKey(5), model_threeRegime,
occupancy.values, speed.values)
num_warmup, num_samples = 1000, 1000
samples_threeRegime = mcmc(num_warmup, num_samples, init_params,
potential_fn=potential_fn,
constrain_fn=constrain_fn,
print_summary=False)
In [14]:
summary(samples_threeRegime, prob=0.95)
Constant - Exponential model¶
In [15]:
def ConstantExponential(occupancy, speed):
α1 = sample('α1', dist.Normal(0, 10))
β1 = sample('β1', dist.Normal(0, 10))
theta = sample('theta', dist.Normal(0, 10))
σ1 = sample('σ1', dist.HalfCauchy(0.5)) #Noise
σ2 = sample('σ2', dist.HalfCauchy(0.5))
switch_point1 = sample("switch_point1", dist.Uniform(occupancy.min(), occupancy.max()))
mu1 = β1
mu2 = α1 * np.exp(-occupancy/theta)
w = cumulative_normal(occupancy - switch_point1)
mu = weight_f(w, mu1, mu2)
sigma = weight_f(w, σ1, σ2)
return sample("obs", dist.Normal(mu, sigma), obs = speed)
In [16]:
'''`hmc` is implemented (default).'''
init_params, potential_fn, constrain_fn = initialize_model(random.PRNGKey(10), ConstantExponential,
occupancy.values, speed.values)
num_warmup, num_samples = 1000, 1000
ConstantExponentials = mcmc(num_warmup, num_samples, init_params,
sampler='hmc',
potential_fn=potential_fn,
constrain_fn=constrain_fn,
#target_accept_prob=0.9, max_tree_depth=12,
print_summary=False)
In [17]:
summary(ConstantExponentials, prob=0.95)
Linear Exponential¶
In [18]:
def LinearExponential(occupancy, speed):
α1 = sample('α1', dist.Normal(0, 10))
β1 = sample('β1', dist.Normal(0, 10))
β2 = sample('β2', dist.Normal(0, 10))
theta = sample('theta', dist.Normal(0, 10))
σ1 = sample('σ1', dist.HalfCauchy(0.5)) #Noise
σ2 = sample('σ2', dist.HalfCauchy(0.5))
switch_point1 = sample("switch_point1", dist.Uniform(occupancy.min(), occupancy.max()))
mu1 = β1 + β2 * occupancy
mu2 = α1 * np.exp(-occupancy/theta)
w = cumulative_normal(occupancy - switch_point1)
mu = weight_f(w, mu1, mu2)
sigma = weight_f(w, σ1, σ2)
return sample("obs", dist.Normal(mu, sigma), obs = speed)
In [19]:
'''`hmc` is implemented (default).'''
init_params, potential_fn, constrain_fn = initialize_model(random.PRNGKey(10), LinearExponential,
occupancy.values, speed.values)
num_warmup, num_samples = 1000, 1000
LinearExponentialS = mcmc(num_warmup, num_samples, init_params,
sampler='hmc',
potential_fn=potential_fn,
constrain_fn=constrain_fn,
#target_accept_prob=0.9, max_tree_depth=12,
print_summary=False)
In [20]:
summary(LinearExponentialS, prob=0.95)
Modified Greenberg model¶
In [21]:
def ModifiedGreenbergmodel(occupancy, speed):
α1 = sample('α1', dist.Normal(0, 10))
β1 = sample('β1', dist.Normal(0, 10))
β2 = sample('β2', dist.Normal(0, 10))
σ1 = sample('σ1', dist.HalfCauchy(0.5)) #Noise
σ2 = sample('σ2', dist.HalfCauchy(0.5))
switch_point1 = sample("switch_point1", dist.Uniform(occupancy.min(), occupancy.max()))
mu1 = α1
mu2 = β1 + β2 * np.log(occupancy)
w = cumulative_normal(occupancy - switch_point1)
mu = weight_f(w, mu1, mu2)
sigma = weight_f(w, σ1, σ2)
return sample("obs", dist.Normal(mu, sigma), obs = speed)
In [22]:
'''`hmc` is implemented (default).'''
# Start from this source of randomness. We will split keys for subsequent operations.
rng = random.PRNGKey(5)
rng_, rng = random.split(rng)
init_params, potential_fn, constrain_fn = initialize_model(rng_, ModifiedGreenbergmodel,
occupancy.values, speed.values)
num_warmup, num_samples = 1000, 1000
ModifiedGreenbergsamples = mcmc(num_warmup, num_samples, init_params,
sampler='hmc',
potential_fn=potential_fn,
constrain_fn=constrain_fn,
#target_accept_prob=0.9, max_tree_depth=12,
print_summary=False)
In [23]:
summary(ModifiedGreenbergsamples, prob=0.95)
Eddies Model¶
In [24]:
def Eddiesmodel(occupancy, speed):
α1 = sample('α1', dist.Normal(0, 10))
theta = sample('theta', dist.Normal(0, 10))
β1 = sample('β1', dist.Normal(0, 10))
β2 = sample('β2', dist.Normal(0, 10))
σ1 = sample('σ1', dist.HalfCauchy(0.5)) #Noise
σ2 = sample('σ2', dist.HalfCauchy(0.5))
switch_point1 = sample("switch_point1", dist.Uniform(occupancy.min(), occupancy.max()))
mu1 = α1 * np.exp(-occupancy/theta)
mu2 = β1 + β2 * np.log(occupancy)
w = sigmoid_f(occupancy - switch_point1)
mu = weight_f(w, mu1, mu2)
sigma = weight_f(w, σ1, σ2)
return sample("obs", dist.Normal(mu, sigma), obs = speed)
In [25]:
'''`hmc` is implemented (default).'''
# Start from this source of randomness. We will split keys for subsequent operations.
rng = random.PRNGKey(2)
rng_, rng = random.split(rng)
init_params, potential_fn, constrain_fn = initialize_model(rng_, Eddiesmodel,
occupancy.values, speed.values)
num_warmup, num_samples = 1000, 1000
Eddies_samples = mcmc(num_warmup, num_samples, init_params,
sampler='hmc',
potential_fn=potential_fn,
constrain_fn=constrain_fn,
#target_accept_prob=0.9, max_tree_depth=12,
print_summary=False)
In [26]:
summary(Eddies_samples, prob=0.95)
Linear logarithm Model¶
In [27]:
def Linearlogarithm(occupancy, speed):
α1 = sample('α1', dist.Normal(0, 10))
α2 = sample('α2', dist.Normal(0, 10))
β1 = sample('β1', dist.Normal(0, 10))
β2 = sample('β2', dist.Normal(0, 10))
σ1 = sample('σ1', dist.HalfCauchy(0.5)) #Noise
σ2 = sample('σ2', dist.HalfCauchy(0.5))
switch_point1 = sample("switch_point1", dist.Uniform(occupancy.min(), occupancy.max()))
mu1 = α1 + α2 * occupancy
mu2 = β1 + β2 * np.log(occupancy)
w = cumulative_normal(occupancy - switch_point1)
mu = weight_f(w, mu1, mu2)
sigma = weight_f(w, σ1, σ2)
return sample("obs", dist.Normal(mu, sigma), obs = speed)
In [28]:
'''`hmc` is implemented (default).'''
# Start from this source of randomness. We will split keys for subsequent operations.
rng = random.PRNGKey(4)
rng_, rng = random.split(rng)
init_params, potential_fn, constrain_fn = initialize_model(rng_, Linearlogarithm,
occupancy.values, speed.values)
num_warmup, num_samples = 1000, 1000
Linearlogarithmsamples = mcmc(num_warmup, num_samples, init_params,
sampler='hmc',
potential_fn=potential_fn,
constrain_fn=constrain_fn,
#target_accept_prob=0.9, max_tree_depth=12,
print_summary=False)
In [29]:
summary(Linearlogarithmsamples, prob=0.95)
Posterior Predictive Distribution¶
In [31]:
def LinearExponential_predi(switch_point, α1,theta, β1, β2, occupancy):
mu1 = β1 + β2 * occupancy
mu2 = α1 * np.exp(-occupancy/theta)
w = cumulative_normal(occupancy - switch_point)
mu = weight_f(w, mu1, mu2)
return mu
def TwoRegime_predi(switch_point, α1, α2, β1, β2, occupancy):
mu1 = α1 + β1 * occupancy
mu2 = α2 + β2 * occupancy
w = cumulative_normal(occupancy - switch_point)
mu = weight_f(w, mu1, mu2)
return mu
def ThreeRegime_predi(switch_point1, switch_point2, α1, α2, α3, β1, β2, β3, occupancy):
mu1 = α1 + β1 * occupancy
mu2 = α2 + β2 * occupancy
mu3 = α3 + β3 * occupancy
w1 =cumulative_normal(occupancy - switch_point1)
w2 = cumulative_normal(occupancy - switch_point2)
mu = weight_three(w1, w2, mu1, mu2, mu3)
return mu
def ConstantExponential_predi(switch_point, α1,theta, β1, occupancy):
mu1 = β1
mu2 = α1 * np.exp(-occupancy/theta)
w = cumulative_normal(occupancy - switch_point)
mu = weight_f(w, mu1, mu2)
return mu
def ModifiedGreenberg_predi(switch_point, α, β1, β2, occupancy):
mu1 = α
mu2 = β1 + β2 * np.log(occupancy)
w = cumulative_normal(occupancy - switch_point)
mu = weight_f(w, mu1, mu2)
return mu
def Eddiesmodel_predi(switch_point, α1, theta, β1, β2, occupancy):
mu1 = α1 * np.exp(-occupancy/theta)
mu2 = β1 + β2 * np.log(occupancy)
w =sigmoid_f(occupancy - switch_point)
mu = weight_f(w, mu1, mu2)
return mu
def Linearlogarithm_predi(switch_point, α1, α2, β1, β2, occupancy):
mu1 = α1 + α2 * occupancy
mu2 = β1 + β2 * np.log(occupancy)
w = cumulative_normal(occupancy - switch_point)
mu = weight_f(w, mu1, mu2)
return mu
In [32]:
f, ax = plt.subplots(figsize=(14, 16))
plt.subplot(421)
plt.scatter(occupancy,speed, color='r', s=50, label = 'Data scatter')
#plt.xlabel("Occupancy (%)"),
plt.ylabel("Traffic speed (mph)")
x = np.linspace(occupancy.min(), occupancy.max(), 1000)
for i in range(100):
sns.lineplot(x, TwoRegime_predi(samples_twoRegime['switch_point'][i],
samples_twoRegime['α1'][i],
samples_twoRegime['α2'][i],
samples_twoRegime['β1'][i],
samples_twoRegime['β2'][i],
x),
color="k",
alpha=0.3)
occupancyBreakdown = samples_twoRegime['switch_point'].mean()
BreakpointSpeed = TwoRegime_predi(samples_twoRegime['switch_point'].mean(),
samples_twoRegime['α1'].mean(),
samples_twoRegime['α2'].mean(),
samples_twoRegime['β1'].mean(),
samples_twoRegime['β2'].mean(),
occupancyBreakdown)
plt.vlines(x=occupancyBreakdown, ymin=0, ymax=BreakpointSpeed, color='g',label='Occupancy threshold = {:.2f} %'.format(occupancyBreakdown),linewidth=2,linestyle='-.')
plt.hlines(BreakpointSpeed, xmin=0, xmax=occupancyBreakdown, color='b',label='Breakpoint speed = {:.2f} mph'.format(BreakpointSpeed),linewidth=2,linestyle='--')
plt.xlim(0, 65)
plt.ylim(0, 85)
plt.legend(loc=1,fontsize=12, frameon=False);
plt.title("Two-regime model", fontsize=14)
##########################@@@@ Three-regime model @@@@@@###########################################################
plt.subplot(422)
plt.scatter(occupancy,speed, color='b', s=50, label = 'Data scatter')
#plt.xlabel("Occupancy (%)"), plt.ylabel("Traffic speed")
# plot posterior predictive lines
x = np.linspace(occupancy.min(), occupancy.max(), 1000)
for i in range(100):
sns.lineplot(x, ThreeRegime_predi(samples_threeRegime['switch_point1'][i],
samples_threeRegime['switch_point2'][i],
samples_threeRegime['α1'][i],
samples_threeRegime['α2'][i],
samples_threeRegime['α3'][i],
samples_threeRegime['β1'][i],
samples_threeRegime['β2'][i],
samples_threeRegime['β3'][i],
x),
color="chocolate",
alpha=0.3)
occupancyBreakdown = samples_threeRegime['switch_point1'].mean()
BreakpointSpeed = ThreeRegime_predi(samples_threeRegime['switch_point1'][i],
samples_threeRegime['switch_point2'][i],
samples_threeRegime['α1'][i],
samples_threeRegime['α2'][i],
samples_threeRegime['α3'][i],
samples_threeRegime['β1'][i],
samples_threeRegime['β2'][i],
samples_threeRegime['β3'][i],
occupancyBreakdown)
occupancyBreakdown2 = samples_threeRegime['switch_point2'].mean()
BreakpointSpeed2 = ThreeRegime_predi(samples_threeRegime['switch_point1'].mean(),
samples_threeRegime['switch_point2'].mean(),
samples_threeRegime['α1'].mean(),
samples_threeRegime['α2'].mean(),
samples_threeRegime['α3'].mean(),
samples_threeRegime['β1'].mean(),
samples_threeRegime['β2'].mean(),
samples_threeRegime['β3'].mean(),
occupancyBreakdown2)
plt.vlines(x=occupancyBreakdown, ymin=0, ymax=BreakpointSpeed, color='r',label='Occupancy threshold 2 = {:.2f} %'.format(occupancyBreakdown),linewidth=2,linestyle='-.')
plt.hlines(BreakpointSpeed, xmin=0, xmax=occupancyBreakdown, color='g',label='Breakpoint speed 2 = {:.2f} mph'.format(BreakpointSpeed),linewidth=2,linestyle='--')
plt.vlines(x=occupancyBreakdown2, ymin=0, ymax=BreakpointSpeed2, color='#c71585',label='Occupancy threshold 1 = {:.2f} %'.format(occupancyBreakdown2),linewidth=2,linestyle='-.')
plt.hlines(BreakpointSpeed2, xmin=0, xmax=occupancyBreakdown2, color='orange',label='Breakpoint speed 1 = {:.2f} mph'.format(BreakpointSpeed2),linewidth=2,linestyle='--')
plt.xlim(0, 65)
plt.ylim(0, 85)
plt.legend(loc=1,fontsize=12, frameon=False);
plt.title("Three-regime model", fontsize=14)
##########################@@@@ Modified Greenberg model@@@@@@###########################################################
plt.subplot(423)
plt.scatter(occupancy,speed, color='m', s=50, label = 'Data scatter')
#plt.xlabel("Occupancy (%)")
plt.ylabel("Traffic speed (mph)")
# plot posterior predictive lines
x = np.linspace(occupancy.min(), occupancy.max(), 1000)
for i in range(100):
sns.lineplot(x, ModifiedGreenberg_predi(ModifiedGreenbergsamples['switch_point1'][i],
ModifiedGreenbergsamples['α1'][i],
ModifiedGreenbergsamples['β1'][i],
ModifiedGreenbergsamples['β2'][i],
x),
color="k",
alpha=0.3)
occupancyBreakdown = ModifiedGreenbergsamples['switch_point1'].mean()
BreakpointSpeed = ModifiedGreenberg_predi(ModifiedGreenbergsamples['switch_point1'].mean(),
ModifiedGreenbergsamples['α1'].mean(),
ModifiedGreenbergsamples['β1'].mean(),
ModifiedGreenbergsamples['β2'].mean(),
occupancyBreakdown)
plt.vlines(x=occupancyBreakdown, ymin=0, ymax=BreakpointSpeed, color='r',label='Occupancy threshold = {:.2f} %'.format(occupancyBreakdown),linewidth=2,linestyle='-.')
plt.hlines(BreakpointSpeed, xmin=0, xmax=occupancyBreakdown, color='b',label='Breakpoint speed = {:.2f} mph'.format(BreakpointSpeed),linewidth=2,linestyle='--')
plt.xlim(0, 65)
plt.ylim(0, 85)
plt.legend(loc=1,fontsize=12, frameon=False);
plt.title("Modified Greenberg model", fontsize=14)
##########################@@@@ Eddies model@@@@@@###########################################################
plt.subplot(424)
plt.scatter(occupancy,speed, color='tan', s=50, label = 'Data scatter')
#plt.xlabel("Occupancy (%)"), plt.ylabel("Traffic speed")
# plot posterior predictive lines
x = np.linspace(occupancy.min(), occupancy.max(), 1000)
for i in range(100):
sns.lineplot(x, Eddiesmodel_predi(Eddies_samples['switch_point1'][i],
Eddies_samples['α1'][i],
Eddies_samples['theta'][i],
Eddies_samples['β1'][i],
Eddies_samples['β2'][i],
x),
color="k",
alpha=0.3)
occupancyBreakdown = Eddies_samples['switch_point1'].mean()
BreakpointSpeed = Eddiesmodel_predi(Eddies_samples['switch_point1'].mean(),
Eddies_samples['α1'].mean(),
Eddies_samples['theta'].mean(),
Eddies_samples['β1'].mean(),
Eddies_samples['β2'].mean(),
occupancyBreakdown)
plt.vlines(x=occupancyBreakdown, ymin=0, ymax=BreakpointSpeed, color='r',label='Occupancy threshold = {:.2f} %'.format(occupancyBreakdown),linewidth=2,linestyle='-.')
plt.hlines(BreakpointSpeed, xmin=0, xmax=occupancyBreakdown, color='b',label='Breakpoint speed = {:.2f} mph'.format(BreakpointSpeed),linewidth=2,linestyle='--')
plt.xlim(0, 65)
plt.ylim(0, 85)
plt.legend(loc=1,fontsize=12, frameon=False);
plt.title("Edie model", fontsize=14)
##########################@@@@ Constant-Exponential model@@@@@@###########################################################
plt.subplot(425)
plt.scatter(occupancy,speed, color='k', s=50, label = 'Data scatter')
#plt.xlabel("Occupancy (%)")
plt.ylabel("Traffic speed (mph)")
# plot posterior predictive lines
x = np.linspace(occupancy.min(), occupancy.max(), 1000)
for i in range(100):
sns.lineplot(x, ConstantExponential_predi(ConstantExponentials['switch_point1'][i],
ConstantExponentials['α1'][i],
ConstantExponentials['theta'][i],
ConstantExponentials['β1'][i],
x),
color="chocolate",
alpha=0.3)
occupancyBreakdown = ConstantExponentials['switch_point1'].mean()
BreakpointSpeed = ConstantExponential_predi(ConstantExponentials['switch_point1'].mean(),
ConstantExponentials['α1'].mean(),
ConstantExponentials['theta'].mean(),
ConstantExponentials['β1'].mean(),
occupancyBreakdown)
plt.vlines(x=occupancyBreakdown, ymin=0, ymax=BreakpointSpeed, color='r',label='Occupancy threshold = {:.2f} %'.format(occupancyBreakdown),linewidth=2,linestyle='-.')
plt.hlines(BreakpointSpeed, xmin=0, xmax=occupancyBreakdown, color='b',label='Breakpoint speed = {:.2f} mph'.format(BreakpointSpeed),linewidth=2,linestyle='--')
plt.xlim(0, 65)
plt.ylim(0, 85)
plt.legend(loc=1,fontsize=12, frameon=False);
plt.title("Constant-Exponential model", fontsize=14)
##########################@@@@ Linear-Exponential model@@@@@@###########################################################
plt.subplot(426)
plt.scatter(occupancy,speed, color='g', s=50, label = 'Data scatter')
plt.xlabel("Occupancy (%)"), #plt.ylabel("Traffic speed")
# plot posterior predictive lines
x = np.linspace(occupancy.min(), occupancy.max(), 1000)
for i in range(100):
sns.lineplot(x, LinearExponential_predi(LinearExponentialS['switch_point1'][i],
LinearExponentialS['α1'][i],
LinearExponentialS['theta'][i],
LinearExponentialS['β1'][i],
LinearExponentialS['β2'][i],
x),
color="r",
alpha=0.3)
occupancyBreakdown = LinearExponentialS['switch_point1'].mean()
BreakpointSpeed = LinearExponential_predi(LinearExponentialS['switch_point1'].mean(),
LinearExponentialS['α1'].mean(),
LinearExponentialS['theta'].mean(),
LinearExponentialS['β1'].mean(),
LinearExponentialS['β2'].mean(),
occupancyBreakdown)
plt.vlines(x=occupancyBreakdown, ymin=0, ymax=BreakpointSpeed, color='r',label='Occupancy threshold = {:.2f} %'.format(occupancyBreakdown),linewidth=2,linestyle='-.')
plt.hlines(BreakpointSpeed, xmin=0, xmax=occupancyBreakdown, color='b',label='Breakpoint speed = {:.2f} mph'.format(BreakpointSpeed),linewidth=2,linestyle='--')
plt.xlim(0, 65)
plt.ylim(0, 85)
plt.legend(loc=1,fontsize=12, frameon=False);
plt.title("Linear-Exponential model", fontsize=14)
##########################@@@@ Linear-Logarithmic model@@@@@@###########################################################
plt.subplot(427)
plt.scatter(occupancy,speed, color='grey', s=50, label = 'Data scatter')
plt.xlabel("Occupancy (%)"), plt.ylabel("Traffic speed (mph)")
# plot posterior predictive lines
x = np.linspace(occupancy.min(), occupancy.max(), 1000)
for i in range(100):
sns.lineplot(x, Linearlogarithm_predi(Linearlogarithmsamples['switch_point1'][i],
Linearlogarithmsamples['α1'][i],
Linearlogarithmsamples['α2'][i],
Linearlogarithmsamples['β1'][i],
Linearlogarithmsamples['β2'][i],
x),
color="chocolate",
alpha=0.3)
occupancyBreakdown = Linearlogarithmsamples['switch_point1'].mean()
BreakpointSpeed = Linearlogarithm_predi(Linearlogarithmsamples['switch_point1'][i],
Linearlogarithmsamples['α1'][i],
Linearlogarithmsamples['α2'][i],
Linearlogarithmsamples['β1'][i],
Linearlogarithmsamples['β2'][i],
occupancyBreakdown)
plt.vlines(x=occupancyBreakdown, ymin=0, ymax=BreakpointSpeed, color='r',label='Occupancy threshold = {:.2f} %'.format(occupancyBreakdown),linewidth=2,linestyle='-.')
plt.hlines(BreakpointSpeed, xmin=0, xmax=occupancyBreakdown, color='b',label='Breakpoint speed = {:.2f} mph'.format(BreakpointSpeed),linewidth=2,linestyle='--')
plt.xlim(0, 65)
plt.ylim(0, 85)
plt.legend(loc=1,fontsize=12, frameon=False);
plt.title("Linear-Logarithmic model", fontsize=14)
plt.tight_layout()
Posterior distribution of Breakpoints¶
In [33]:
from numpyro.diagnostics import hpdi
plt.style.use("seaborn-dark")
In [34]:
f, ax = plt.subplots(figsize=(14, 16))
plt.subplot(421)
cred_mass = 0.95
breakpoints = TwoRegime_predi(samples_twoRegime['switch_point'],
samples_twoRegime['α1'],
samples_twoRegime['α2'],
samples_twoRegime['β1'],
samples_twoRegime['β2'],
samples_twoRegime['switch_point'])
MEAN = breakpoints.mean()
# Computes “highest posterior density interval” (HPDI)
HDI = hpdi(breakpoints, cred_mass)
sns.distplot(breakpoints, hist = False, kde = True,
kde_kws={"lw": 2})
#plt.ylabel("Data density")
plt.title("Two-regime speed threshold")
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
plot_height = frame.get_ylim()[1]
plt.plot(HDI, [0, 0], linewidth=6, color='r')
plt.text(HDI[1],
plot_height * 0.07,
HDI[1].round(2),
horizontalalignment="center",
)
plt.text(HDI[0],
plot_height * 0.07,
HDI[0].round(2),
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.3,
'95%' + " HPD",
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.8,
'Mean = %.2f' % MEAN,
horizontalalignment="center",
)
# plt.xlim(0, 65)
# plt.ylim(-0.01, 1.30)
# ##########################@@@@ Three-regime model @@@@@@###########################################################
plt.subplot(422)
cred_mass = 0.95
breakpoints1 = ThreeRegime_predi(samples_threeRegime['switch_point1'],
samples_threeRegime['switch_point2'],
samples_threeRegime['α1'],
samples_threeRegime['α2'],
samples_threeRegime['α3'],
samples_threeRegime['β1'],
samples_threeRegime['β2'],
samples_threeRegime['β3'],
samples_threeRegime['switch_point1'])
MEAN = breakpoints1.mean()
plot_height = ax.get_ylim()[1]
# Computes “highest posterior density interval” (HPDI)
HDI = hpdi(breakpoints1, cred_mass)
sns.distplot(breakpoints1, hist = False, kde = True,
kde_kws={"lw": 2})
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
plot_height = frame.get_ylim()[1]
plt.title("Three-regime speed threshold 1")
plt.plot(HDI, [0, 0], linewidth=6, color='r')
plt.text(HDI[1],
plot_height * 0.07,
HDI[1].round(2),
horizontalalignment="center",
)
plt.text(HDI[0],
plot_height * 0.07,
HDI[0].round(2),
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.3,
'95%' + " HPD",
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.8,
'Mean = %.2f' % MEAN,
horizontalalignment="center",
)
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
# plt.xlim(0, 65)
# plt.ylim(-0.01, 1.75)
# ##########################@@@@ Three-regime model @@@@@@###########################################################
plt.subplot(423)
cred_mass = 0.95
breakpoints2 = ThreeRegime_predi(samples_threeRegime['switch_point1'],
samples_threeRegime['switch_point2'],
samples_threeRegime['α1'],
samples_threeRegime['α2'],
samples_threeRegime['α3'],
samples_threeRegime['β1'],
samples_threeRegime['β2'],
samples_threeRegime['β3'],
samples_threeRegime['switch_point2'])
MEAN = breakpoints2.mean()
# plot_height = ax.get_ylim()[1]
# Computes “highest posterior density interval” (HPDI)
HDI = hpdi(breakpoints2, cred_mass)
sns.distplot(breakpoints2, hist = False, kde = True,
kde_kws={"lw": 2})
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
plot_height = frame.get_ylim()[1]
plt.title("Three-regime speed threshold 2")
plt.plot(HDI, [0, 0], linewidth=6, color='r')
plt.text(HDI[1],
plot_height * 0.07,
HDI[1].round(2),
horizontalalignment="center",
)
plt.text(HDI[0],
plot_height * 0.07,
HDI[0].round(2),
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.3,
'95%' + " HPD",
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.8,
'Mean = %.2f' % MEAN,
horizontalalignment="center",
)
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
# plt.xlim(0, 65)
# plt.ylim(-0.01, 1.75)
# ##########################@@@@ Modified Greenberg model@@@@@@###########################################################
plt.subplot(424)
cred_mass = 0.95
breakpoints = ModifiedGreenberg_predi(ModifiedGreenbergsamples['switch_point1'],
ModifiedGreenbergsamples['α1'],
ModifiedGreenbergsamples['β1'],
ModifiedGreenbergsamples['β2'],
ModifiedGreenbergsamples['switch_point1'])
MEAN = breakpoints.mean()
# plot_height = ax.get_ylim()[1]
# Computes “highest posterior density interval” (HPDI)
HDI = hpdi(breakpoints, cred_mass)
sns.distplot(breakpoints, hist = False, kde = True,
kde_kws={"lw": 2})
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
plot_height = frame.get_ylim()[1]
plt.title("Modified Greenberg model")
plt.plot(HDI, [0, 0], linewidth=6, color='r')
plt.text(HDI[1],
plot_height * 0.07,
HDI[1].round(2),
horizontalalignment="center",
)
plt.text(HDI[0],
plot_height * 0.07,
HDI[0].round(2),
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.3,
'95%' + " HPD",
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.8,
'Mean = %.2f' % MEAN,
horizontalalignment="center",
)
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
# plt.xlim(0, 65)
# plt.ylim(-0.01, 1.85)
# ##########################@@@@ Eddies model@@@@@@###########################################################
plt.subplot(425)
breakpoints = Eddiesmodel_predi(Eddies_samples['switch_point1'],
Eddies_samples['α1'],
Eddies_samples['theta'],
Eddies_samples['β1'],
Eddies_samples['β2'],
Eddies_samples['switch_point1'])
MEAN = breakpoints.mean()
# plot_height = ax.get_ylim()[1]
# Computes “highest posterior density interval” (HPDI)
HDI = hpdi(breakpoints, cred_mass)
sns.distplot(breakpoints, hist = False, kde = True,
kde_kws={"lw": 2})
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
plot_height = frame.get_ylim()[1]
plt.plot(HDI, [0, 0], linewidth=6, color='r')
plt.text(HDI[1],
plot_height * 0.07,
HDI[1].round(2),
horizontalalignment="center",
)
plt.text(HDI[0],
plot_height * 0.07,
HDI[0].round(2),
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.3,
'95%' + " HPD",
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.8,
'Mean = %.2f' % MEAN,
horizontalalignment="center",
)
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
# plt.xlim(0, 65)
#plt.ylim(-0.01, 1.75)
plt.title("Edie model")
# ##########################@@@@ Constant-Exponential model@@@@@@###########################################################
plt.subplot(426)
breakpoints = ConstantExponential_predi(ConstantExponentials['switch_point1'],
ConstantExponentials['α1'],
ConstantExponentials['theta'],
ConstantExponentials['β1'],
ConstantExponentials['switch_point1'])
MEAN = breakpoints.mean()
plot_height = ax.get_ylim()[1]
# Computes “highest posterior density interval” (HPDI)
HDI = hpdi(breakpoints, cred_mass)
sns.distplot(breakpoints, hist = False, kde = True,
kde_kws={"lw": 2})
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
plot_height = frame.get_ylim()[1]
plt.plot(HDI, [0, 0], linewidth=6, color='r')
plt.text(HDI[1],
plot_height * 0.07,
HDI[1].round(2),
horizontalalignment="center",
)
plt.text(HDI[0],
plot_height * 0.07,
HDI[0].round(2),
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.3,
'95%' + " HPDI",
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.8,
'Mean = %.2f' % MEAN,
horizontalalignment="center",
)
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
# plt.xlim(0, 65)
#plt.ylim(-0.01, 1.75)
plt.title("Constant-Exponential model")
# ##########################@@@@ Linear-Exponential model@@@@@@###########################################################
plt.subplot(427)
breakpoints = LinearExponential_predi(LinearExponentialS['switch_point1'],
LinearExponentialS['α1'],
LinearExponentialS['theta'],
LinearExponentialS['β1'],
LinearExponentialS['β2'],
LinearExponentialS['switch_point1'])
MEAN = breakpoints.mean()
# plot_height = ax.get_ylim()[1]
# Computes “highest posterior density interval” (HPDI)
HDI = hpdi(breakpoints, cred_mass)
sns.distplot(breakpoints, hist = False, kde = True,
kde_kws={"lw": 2})
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
plot_height = frame.get_ylim()[1]
plt.xlabel("Traffic speed (mph)")
plt.plot(HDI, [0, 0], linewidth=6, color='r')
plt.text(HDI[1],
plot_height * 0.07,
HDI[1].round(2),
horizontalalignment="center",
)
plt.text(HDI[0],
plot_height * 0.07,
HDI[0].round(2),
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.3,
'95%' + " HPD",
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.8,
'Mean = %.2f' % MEAN,
horizontalalignment="center",
)
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
# plt.xlim(0, 65)
#plt.ylim(-0.01, 1.75)
plt.title("Linear-Exponential model")
# ##########################@@@@ Linear-Logarithmic model@@@@@@###########################################################
plt.subplot(428)
breakpoints = Linearlogarithm_predi(Linearlogarithmsamples['switch_point1'],
Linearlogarithmsamples['α1'],
Linearlogarithmsamples['α2'],
Linearlogarithmsamples['β1'],
Linearlogarithmsamples['β2'],
Linearlogarithmsamples['switch_point1'])
MEAN = breakpoints.mean()
# Computes “highest posterior density interval” (HPDI)
HDI = hpdi(breakpoints, cred_mass)
sns.distplot(breakpoints, hist = False, kde = True,
kde_kws={"lw": 2})
# plot_height = ax.get_ylim()[1]
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
plot_height = frame.get_ylim()[1]
plt.xlabel("Traffic speed (mph)")
plt.plot(HDI, [0, 0], linewidth=8, color='r')
plt.text(HDI[1],
plot_height * 0.07,
HDI[1].round(2),
horizontalalignment="center",
)
plt.text(HDI[0],
plot_height * 0.07,
HDI[0].round(2),
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.3,
'95%' + " HPD",
horizontalalignment="center",
)
plt.text(
(HDI[0] + HDI[1]) / 2,
plot_height * 0.8,
'Mean = %.2f' % MEAN,
horizontalalignment="center",
)
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
# plt.xlim(0, 65)
plt.title("Linear-Logarithmic model")
plt.tight_layout()
Prediction Accuracy¶
In [35]:
import scipy
def rsquared(x, y):
""" Return R^2 where x and y are array-like."""
slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(x, y)
return r_value**2
In [36]:
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error
In [39]:
sns.set_context("notebook", font_scale= 1.7)
f, ax = plt.subplots(figsize=(16, 16))
plt.subplot(331)
posterior_mu = TwoRegime_predi(samples_twoRegime['switch_point'].mean(),
samples_twoRegime['α1'].mean(),
samples_twoRegime['α2'].mean(),
samples_twoRegime['β1'].mean(),
samples_twoRegime['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'r'},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2,fontsize=13);
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("Observed traffic speed $(mph)$", fontsize=18)
plt.title("Two-regime model", fontsize=18)
plt.xlabel('')
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy = pd.DataFrame({'Metrics':['Mean squared error',
'Mean absolute error',
'Coefficient of determination',
'Pearson correlation']})
Accuracy['TwoRegime']= [mean_squarederror, meanabsolut, R2, r2]
#################################ThreeRegime###########################################
plt.subplot(332)
posterior_mu = ThreeRegime_predi(samples_threeRegime['switch_point1'].mean(),
samples_threeRegime['switch_point2'].mean(),
samples_threeRegime['α1'].mean(),
samples_threeRegime['α2'].mean(),
samples_threeRegime['α3'].mean(),
samples_threeRegime['β1'].mean(),
samples_threeRegime['β2'].mean(),
samples_threeRegime['β3'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
sns.set_context("notebook", font_scale= 1.7)
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'b'},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2,fontsize=13);
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("", fontsize=18)
plt.xlabel("", fontsize=18)
plt.title("Three-regime model", fontsize=18)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['ThreeRegime'] = [mean_squarederror, meanabsolut, R2, r2]
#################################ModifiedGreenberg###########################################
plt.subplot(333)
posterior_mu = ModifiedGreenberg_predi(ModifiedGreenbergsamples['switch_point1'].mean(),
ModifiedGreenbergsamples['α1'].mean(),
ModifiedGreenbergsamples['β1'].mean(),
ModifiedGreenbergsamples['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
sns.set_context("notebook", font_scale= 1.7)
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'m'},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2,fontsize=13);
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("", fontsize=18)
plt.xlabel("", fontsize=18)
plt.title("Modified Greenberg model", fontsize=18)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['ModifiedGreenberg'] = [mean_squarederror, meanabsolut, R2, r2]
#################################Eddiesmodel###########################################
plt.subplot(334)
posterior_mu = Eddiesmodel_predi(Eddies_samples['switch_point1'].mean(),
Eddies_samples['α1'].mean(),
Eddies_samples['theta'].mean(),
Eddies_samples['β1'].mean(),
Eddies_samples['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
sns.set_context("notebook", font_scale= 1.7)
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'m'},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2,fontsize=13);
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("Observed traffic speed $(mph)$", fontsize=18)
plt.xlabel("", fontsize=18)
plt.title("Eddies model", fontsize=18)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['Eddiesmodel'] = [mean_squarederror, meanabsolut, R2, r2]
#################################ConstantExponential###########################################
plt.subplot(335)
posterior_mu = ConstantExponential_predi(ConstantExponentials['switch_point1'].mean(),
ConstantExponentials['α1'].mean(),
ConstantExponentials['theta'].mean(),
ConstantExponentials['β1'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
sns.set_context("notebook", font_scale= 1.7)
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'k', 'alpha':1},
line_kws={"lw":4,
'ls':'-',
'color':'chocolate', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2,fontsize=13);
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("", fontsize=18)
plt.xlabel("Predicted traffic speed $(mph)$", fontsize=18)
plt.title("Constant-Exponential model", fontsize=18)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['ConstantExponential'] = [mean_squarederror, meanabsolut, R2, r2]
#################################LinearExponential###########################################
plt.subplot(336)
posterior_mu = LinearExponential_predi(LinearExponentialS['switch_point1'].mean(),
LinearExponentialS['α1'].mean(),
LinearExponentialS['theta'].mean(),
LinearExponentialS['β1'].mean(),
LinearExponentialS['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
sns.set_context("notebook", font_scale= 1.7)
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'g', 'alpha':1},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2,fontsize=13);
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("", fontsize=18)
plt.xlabel("Predicted traffic speed $(mph)$", fontsize=18)
plt.title("Linear-Exponential model", fontsize=18)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['LinearExponential'] = [mean_squarederror, meanabsolut, R2, r2]
#################################Linearlogarithm###########################################
plt.subplot(337)
posterior_mu = Linearlogarithm_predi(Linearlogarithmsamples['switch_point1'].mean(),
Linearlogarithmsamples['α1'].mean(),
Linearlogarithmsamples['α2'].mean(),
Linearlogarithmsamples['β1'].mean(),
Linearlogarithmsamples['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
sns.set_context("notebook", font_scale= 1.7)
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'grey', 'alpha':1},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=4,fontsize=13);
plt.title("Linear-Logarithm model", fontsize=18)
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("", fontsize=18)
plt.xlabel("Predicted traffic speed $(mph)$", fontsize=18)
plt.ylabel("Observed traffic speed $(mph)$", fontsize=18)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['Linearlogarithm'] = [mean_squarederror, meanabsolut, R2, r2]
In [38]:
Accuracy
Out[38]:
In [105]:
sns.set_context("notebook", font_scale= 2.2)
f, ax = plt.subplots(figsize=(30, 30))
plt.subplot(441)
posterior_mu = TwoRegime_predi(samples_twoRegime['switch_point'].mean(),
samples_twoRegime['α1'].mean(),
samples_twoRegime['α2'].mean(),
samples_twoRegime['β1'].mean(),
samples_twoRegime['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'r'},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2#,fontsize=13
);
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("Observed traffic speed $(mph)$"#, fontsize=18
)
plt.title("Two-regime model"#, fontsize=18
)
plt.xlabel('')
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy = pd.DataFrame({'Metrics':['Mean squared error',
'Mean absolute error',
'Coefficient of determination',
'Pearson correlation']})
Accuracy['TwoRegime']= [mean_squarederror, meanabsolut, R2, r2]
#################################ThreeRegime###########################################
plt.subplot(442)
posterior_mu = ThreeRegime_predi(samples_threeRegime['switch_point1'].mean(),
samples_threeRegime['switch_point2'].mean(),
samples_threeRegime['α1'].mean(),
samples_threeRegime['α2'].mean(),
samples_threeRegime['α3'].mean(),
samples_threeRegime['β1'].mean(),
samples_threeRegime['β2'].mean(),
samples_threeRegime['β3'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'b'},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2#,fontsize=13
);
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("", fontsize=18)
plt.xlabel("", fontsize=18)
plt.title("Three-regime model"#, fontsize=18
)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['ThreeRegime'] = [mean_squarederror, meanabsolut, R2, r2]
#################################ModifiedGreenberg###########################################
plt.subplot(443)
posterior_mu = ModifiedGreenberg_predi(ModifiedGreenbergsamples['switch_point1'].mean(),
ModifiedGreenbergsamples['α1'].mean(),
ModifiedGreenbergsamples['β1'].mean(),
ModifiedGreenbergsamples['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'m'},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2#,fontsize=13
);
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("", fontsize=18)
plt.xlabel("", fontsize=18)
plt.title("Modified Greenberg model"#, fontsize=18
)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['ModifiedGreenberg'] = [mean_squarederror, meanabsolut, R2, r2]
#################################Eddiesmodel###########################################
plt.subplot(444)
posterior_mu = Eddiesmodel_predi(Eddies_samples['switch_point1'].mean(),
Eddies_samples['α1'].mean(),
Eddies_samples['theta'].mean(),
Eddies_samples['β1'].mean(),
Eddies_samples['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'tan'},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2#,fontsize=13
);
plt.xlim(0, 80)
plt.ylim(0, 80)
#plt.ylabel("Observed traffic speed $(mph)$", fontsize=18)
plt.xlabel("Predicted traffic speed $(mph)$"#, fontsize=18
)
plt.title("Eddies model"#, fontsize=18
)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['Eddiesmodel'] = [mean_squarederror, meanabsolut, R2, r2]
#################################ConstantExponential###########################################
plt.subplot(445)
posterior_mu = ConstantExponential_predi(ConstantExponentials['switch_point1'].mean(),
ConstantExponentials['α1'].mean(),
ConstantExponentials['theta'].mean(),
ConstantExponentials['β1'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'k', 'alpha':1},
line_kws={"lw":4,
'ls':'-',
'color':'chocolate', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2#,fontsize=13
);
plt.xlim(0, 80)
plt.ylim(0, 80)
#plt.ylabel("", fontsize=18)
plt.xlabel("Predicted traffic speed $(mph)$"#, fontsize=18
)
plt.ylabel("Observed traffic speed $(mph)$"#, fontsize=18
)
plt.title("Constant-Exponential model"#, fontsize=18
)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['ConstantExponential'] = [mean_squarederror, meanabsolut, R2, r2]
#################################LinearExponential###########################################
plt.subplot(446)
posterior_mu = LinearExponential_predi(LinearExponentialS['switch_point1'].mean(),
LinearExponentialS['α1'].mean(),
LinearExponentialS['theta'].mean(),
LinearExponentialS['β1'].mean(),
LinearExponentialS['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'g', 'alpha':1},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=2#,fontsize=13
);
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("", fontsize=18)
plt.xlabel("Predicted traffic speed $(mph)$"#, fontsize=18
)
plt.title("Linear-Exponential model"#, fontsize=18
)
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['LinearExponential'] = [mean_squarederror, meanabsolut, R2, r2]
#################################Linearlogarithm###########################################
plt.subplot(447)
posterior_mu = Linearlogarithm_predi(Linearlogarithmsamples['switch_point1'].mean(),
Linearlogarithmsamples['α1'].mean(),
Linearlogarithmsamples['α2'].mean(),
Linearlogarithmsamples['β1'].mean(),
Linearlogarithmsamples['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
R2 = r2_score(Pred.pred,
Pred.Observed).round(2)
r2 = rsquared(Pred.pred,
Pred.Observed).round(2)
sns.regplot(x="pred", y="Observed", data=Pred,
scatter_kws={'color':'grey', 'alpha':1},
line_kws={"lw":4,
'ls':'-',
'color':'k', 'label':'Pearson correlation: {} '.format(r2)});
plt.legend(loc=4#,fontsize=13
);
plt.title("Linear-Logarithm model"#, fontsize=18
)
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.ylabel("", fontsize=18)
plt.xlabel("Predicted traffic speed $(mph)$"#, fontsize=18
)
# plt.ylabel("Observed traffic speed $(mph)$"#, fontsize=18
# )
plt.tight_layout()
mean_squarederror = mean_squared_error(posterior_mu, speed).round(2)
meanabsolut = mean_absolute_error(posterior_mu, speed).round(2)
Accuracy['Linearlogarithm'] = [mean_squarederror, meanabsolut, R2, r2]
In [115]:
dashList = [(5,2),(2,5),(4,10),(3,3,2,2),(5,2,20,2)]
Model Accuracy¶
In [148]:
sns.set_context("notebook", font_scale= 1.7)
f, ax = plt.subplots(figsize=(10, 7))
# plt.subplot(441)
posterior_mu = TwoRegime_predi(samples_twoRegime['switch_point'].mean(),
samples_twoRegime['α1'].mean(),
samples_twoRegime['α2'].mean(),
samples_twoRegime['β1'].mean(),
samples_twoRegime['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
Pred['oupany_ran'] = [0 if 0 <= i <= 20 else
1 if 20 < i <= 30 else
2 if 30 < i <= 40 else
3 if 40 < i <= 50 else
4 if 50 < i <= 60 else
5 if 60 < i <= 70 else
6 for i in Pred.Observed
]
Pred['Label'] = ['0-20' if 0 <= i <= 20 else
'20-30' if 20 < i <= 30 else
'30-40' if 30 < i <= 40 else
'40-50' if 40 < i <= 50 else
'50-60' if 50 < i <= 60 else
'60-70' if 60 < i <= 70 else
'> 70' for i in Pred.Observed
]
cols = ['Mean_squared_error',
'Mean_absolute_error',
'Coefficient_of_determination',
'Pearson_correlation',
'Number_of_observations (n) ']
Mod_a = []
for i in arange(7):
dxxf = Pred[Pred.oupany_ran==i]
observ = len(dxxf)
#print(i)
R2 = r2_score(dxxf.pred,
dxxf.Observed).round(2)
r2 = rsquared(dxxf.pred,
dxxf.Observed).round(2)
mean_squarederror = mean_squared_error(dxxf.pred, dxxf.Observed).round(2)
meanabsolut = mean_absolute_error(dxxf.pred, dxxf.Observed).round(2)
#print(mean_squarederror, meanabsolut, R2, r2, observ)
Mod_a.append([mean_squarederror, meanabsolut, R2, r2, observ])
Accuracy = pd.DataFrame({'Metrics':['0-20',
'20-30',
'30-40',
'40-50',
'50-60',
'60-70',
'> 70']})
TwoRegime_Accuracy = pd.concat([Accuracy,pd.DataFrame(Mod_a, columns=cols)], axis=1)
#################################ThreeRegime###########################################
posterior_mu = ThreeRegime_predi(samples_threeRegime['switch_point1'].mean(),
samples_threeRegime['switch_point2'].mean(),
samples_threeRegime['α1'].mean(),
samples_threeRegime['α2'].mean(),
samples_threeRegime['α3'].mean(),
samples_threeRegime['β1'].mean(),
samples_threeRegime['β2'].mean(),
samples_threeRegime['β3'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
Pred['oupany_ran'] = [0 if 0 <= i <= 20 else
1 if 20 < i <= 30 else
2 if 30 < i <= 40 else
3 if 40 < i <= 50 else
4 if 50 < i <= 60 else
5 if 60 < i <= 70 else
6 for i in Pred.Observed
]
Pred['Label'] = ['0-20' if 0 <= i <= 20 else
'20-30' if 20 < i <= 30 else
'30-40' if 30 < i <= 40 else
'40-50' if 40 < i <= 50 else
'50-60' if 50 < i <= 60 else
'60-70' if 60 < i <= 70 else
'> 70' for i in Pred.Observed
]
cols = ['Mean_squared_error',
'Mean_absolute_error',
'Coefficient_of_determination',
'Pearson_correlation',
'Number_of_observations (n) ']
Mod_a = []
for i in arange(7):
dxxf = Pred[Pred.oupany_ran==i]
observ = len(dxxf)
#print(i)
R2 = r2_score(dxxf.pred,
dxxf.Observed).round(2)
r2 = rsquared(dxxf.pred,
dxxf.Observed).round(2)
mean_squarederror = mean_squared_error(dxxf.pred, dxxf.Observed).round(2)
meanabsolut = mean_absolute_error(dxxf.pred, dxxf.Observed).round(2)
#print(mean_squarederror, meanabsolut, R2, r2, observ)
Mod_a.append([mean_squarederror, meanabsolut, R2, r2, observ])
Accuracy = pd.DataFrame({'Metrics':['0-20',
'20-30',
'30-40',
'40-50',
'50-60',
'60-70',
'> 70']})
ThreeRegime_Accuracy = pd.concat([Accuracy,pd.DataFrame(Mod_a, columns=cols)], axis=1)
#################################ModifiedGreenberg###########################################
posterior_mu = ModifiedGreenberg_predi(ModifiedGreenbergsamples['switch_point1'].mean(),
ModifiedGreenbergsamples['α1'].mean(),
ModifiedGreenbergsamples['β1'].mean(),
ModifiedGreenbergsamples['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
Pred['oupany_ran'] = [0 if 0 <= i <= 20 else
1 if 20 < i <= 30 else
2 if 30 < i <= 40 else
3 if 40 < i <= 50 else
4 if 50 < i <= 60 else
5 if 60 < i <= 70 else
6 for i in Pred.Observed
]
Pred['Label'] = ['0-20' if 0 <= i <= 20 else
'20-30' if 20 < i <= 30 else
'30-40' if 30 < i <= 40 else
'40-50' if 40 < i <= 50 else
'50-60' if 50 < i <= 60 else
'60-70' if 60 < i <= 70 else
'> 70' for i in Pred.Observed
]
cols = ['Mean_squared_error',
'Mean_absolute_error',
'Coefficient_of_determination',
'Pearson_correlation',
'Number_of_observations (n) ']
Mod_a = []
for i in arange(7):
dxxf = Pred[Pred.oupany_ran==i]
observ = len(dxxf)
#print(i)
R2 = r2_score(dxxf.pred,
dxxf.Observed).round(2)
r2 = rsquared(dxxf.pred,
dxxf.Observed).round(2)
mean_squarederror = mean_squared_error(dxxf.pred, dxxf.Observed).round(2)
meanabsolut = mean_absolute_error(dxxf.pred, dxxf.Observed).round(2)
#print(mean_squarederror, meanabsolut, R2, r2, observ)
Mod_a.append([mean_squarederror, meanabsolut, R2, r2, observ])
Accuracy = pd.DataFrame({'Metrics':['0-20',
'20-30',
'30-40',
'40-50',
'50-60',
'60-70',
'> 70']})
ModifiedGreenberg_Accuracy = pd.concat([Accuracy,pd.DataFrame(Mod_a, columns=cols)], axis=1)
#################################Eddiesmodel###########################################
posterior_mu = Eddiesmodel_predi(Eddies_samples['switch_point1'].mean(),
Eddies_samples['α1'].mean(),
Eddies_samples['theta'].mean(),
Eddies_samples['β1'].mean(),
Eddies_samples['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
Pred['oupany_ran'] = [0 if 0 <= i <= 20 else
1 if 20 < i <= 30 else
2 if 30 < i <= 40 else
3 if 40 < i <= 50 else
4 if 50 < i <= 60 else
5 if 60 < i <= 70 else
6 for i in Pred.Observed
]
Pred['Label'] = ['0-20' if 0 <= i <= 20 else
'20-30' if 20 < i <= 30 else
'30-40' if 30 < i <= 40 else
'40-50' if 40 < i <= 50 else
'50-60' if 50 < i <= 60 else
'60-70' if 60 < i <= 70 else
'> 70' for i in Pred.Observed
]
cols = ['Mean_squared_error',
'Mean_absolute_error',
'Coefficient_of_determination',
'Pearson_correlation',
'Number_of_observations (n) ']
Mod_a = []
for i in arange(7):
dxxf = Pred[Pred.oupany_ran==i]
observ = len(dxxf)
#print(i)
R2 = r2_score(dxxf.pred,
dxxf.Observed).round(2)
r2 = rsquared(dxxf.pred,
dxxf.Observed).round(2)
mean_squarederror = mean_squared_error(dxxf.pred, dxxf.Observed).round(2)
meanabsolut = mean_absolute_error(dxxf.pred, dxxf.Observed).round(2)
#print(mean_squarederror, meanabsolut, R2, r2, observ)
Mod_a.append([mean_squarederror, meanabsolut, R2, r2, observ])
Accuracy = pd.DataFrame({'Metrics':['0-20',
'20-30',
'30-40',
'40-50',
'50-60',
'60-70',
'> 70']})
Eddiesmodel_Accuracy = pd.concat([Accuracy,pd.DataFrame(Mod_a, columns=cols)], axis=1)
#################################ConstantExponential###########################################
posterior_mu = ConstantExponential_predi(ConstantExponentials['switch_point1'].mean(),
ConstantExponentials['α1'].mean(),
ConstantExponentials['theta'].mean(),
ConstantExponentials['β1'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
Pred['oupany_ran'] = [0 if 0 <= i <= 20 else
1 if 20 < i <= 30 else
2 if 30 < i <= 40 else
3 if 40 < i <= 50 else
4 if 50 < i <= 60 else
5 if 60 < i <= 70 else
6 for i in Pred.Observed
]
Pred['Label'] = ['0-20' if 0 <= i <= 20 else
'20-30' if 20 < i <= 30 else
'30-40' if 30 < i <= 40 else
'40-50' if 40 < i <= 50 else
'50-60' if 50 < i <= 60 else
'60-70' if 60 < i <= 70 else
'> 70' for i in Pred.Observed
]
cols = ['Mean_squared_error',
'Mean_absolute_error',
'Coefficient_of_determination',
'Pearson_correlation',
'Number_of_observations (n) ']
Mod_a = []
for i in arange(7):
dxxf = Pred[Pred.oupany_ran==i]
observ = len(dxxf)
#print(i)
R2 = r2_score(dxxf.pred,
dxxf.Observed).round(2)
r2 = rsquared(dxxf.pred,
dxxf.Observed).round(2)
mean_squarederror = mean_squared_error(dxxf.pred, dxxf.Observed).round(2)
meanabsolut = mean_absolute_error(dxxf.pred, dxxf.Observed).round(2)
#print(mean_squarederror, meanabsolut, R2, r2, observ)
Mod_a.append([mean_squarederror, meanabsolut, R2, r2, observ])
Accuracy = pd.DataFrame({'Metrics':['0-20',
'20-30',
'30-40',
'40-50',
'50-60',
'60-70',
'> 70']})
ConstantExponential_Accuracy = pd.concat([Accuracy,pd.DataFrame(Mod_a, columns=cols)], axis=1)
#################################LinearExponential###########################################
posterior_mu = LinearExponential_predi(LinearExponentialS['switch_point1'].mean(),
LinearExponentialS['α1'].mean(),
LinearExponentialS['theta'].mean(),
LinearExponentialS['β1'].mean(),
LinearExponentialS['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
Pred['oupany_ran'] = [0 if 0 <= i <= 20 else
1 if 20 < i <= 30 else
2 if 30 < i <= 40 else
3 if 40 < i <= 50 else
4 if 50 < i <= 60 else
5 if 60 < i <= 70 else
6 for i in Pred.Observed
]
Pred['Label'] = ['0-20' if 0 <= i <= 20 else
'20-30' if 20 < i <= 30 else
'30-40' if 30 < i <= 40 else
'40-50' if 40 < i <= 50 else
'50-60' if 50 < i <= 60 else
'60-70' if 60 < i <= 70 else
'> 70' for i in Pred.Observed
]
cols = ['Mean_squared_error',
'Mean_absolute_error',
'Coefficient_of_determination',
'Pearson_correlation',
'Number_of_observations (n) ']
Mod_a = []
for i in arange(7):
dxxf = Pred[Pred.oupany_ran==i]
observ = len(dxxf)
#print(i)
R2 = r2_score(dxxf.pred,
dxxf.Observed).round(2)
r2 = rsquared(dxxf.pred,
dxxf.Observed).round(2)
mean_squarederror = mean_squared_error(dxxf.pred, dxxf.Observed).round(2)
meanabsolut = mean_absolute_error(dxxf.pred, dxxf.Observed).round(2)
#print(mean_squarederror, meanabsolut, R2, r2, observ)
Mod_a.append([mean_squarederror, meanabsolut, R2, r2, observ])
Accuracy = pd.DataFrame({'Metrics':['0-20',
'20-30',
'30-40',
'40-50',
'50-60',
'60-70',
'> 70']})
LinearExponential_Accuracy = pd.concat([Accuracy,pd.DataFrame(Mod_a, columns=cols)], axis=1)
#################################Linearlogarithm###########################################
posterior_mu = Linearlogarithm_predi(Linearlogarithmsamples['switch_point1'].mean(),
Linearlogarithmsamples['α1'].mean(),
Linearlogarithmsamples['α2'].mean(),
Linearlogarithmsamples['β1'].mean(),
Linearlogarithmsamples['β2'].mean(),
occupancy.values)
Pred = pd.DataFrame({'pred':posterior_mu, 'Observed':speed})
Pred['oupany_ran'] = [0 if 0 <= i <= 20 else
1 if 20 < i <= 30 else
2 if 30 < i <= 40 else
3 if 40 < i <= 50 else
4 if 50 < i <= 60 else
5 if 60 < i <= 70 else
6 for i in Pred.Observed
]
Pred['Label'] = ['0-20' if 0 <= i <= 20 else
'20-30' if 20 < i <= 30 else
'30-40' if 30 < i <= 40 else
'40-50' if 40 < i <= 50 else
'50-60' if 50 < i <= 60 else
'60-70' if 60 < i <= 70 else
'> 70' for i in Pred.Observed
]
cols = ['Mean_squared_error',
'Mean_absolute_error',
'Coefficient_of_determination',
'Pearson_correlation',
'Number_of_observations (n) ']
Mod_a = []
for i in arange(7):
dxxf = Pred[Pred.oupany_ran==i]
observ = len(dxxf)
#print(i)
R2 = r2_score(dxxf.pred,
dxxf.Observed).round(2)
r2 = rsquared(dxxf.pred,
dxxf.Observed).round(2)
mean_squarederror = mean_squared_error(dxxf.pred, dxxf.Observed).round(2)
meanabsolut = mean_absolute_error(dxxf.pred, dxxf.Observed).round(2)
#print(mean_squarederror, meanabsolut, R2, r2, observ)
Mod_a.append([mean_squarederror, meanabsolut, R2, r2, observ])
Accuracy = pd.DataFrame({'Metrics':['0-20',
'20-30',
'30-40',
'40-50',
'50-60',
'60-70',
'> 70']})
Linearlogarithm_Accuracy = pd.concat([Accuracy,pd.DataFrame(Mod_a, columns=cols)], axis=1)
########################################################################################################
plt.plot(TwoRegime_Accuracy.Metrics, TwoRegime_Accuracy.Mean_absolute_error,
label='Two-regime', linewidth=2, linestyle =':',marker='*')
plt.plot(ThreeRegime_Accuracy.Metrics, ThreeRegime_Accuracy.Mean_absolute_error,
color='r',label='Three Regime', linewidth=2, dashes=dashList[1],marker='s')
plt.plot(ModifiedGreenberg_Accuracy.Metrics, ModifiedGreenberg_Accuracy.Mean_absolute_error,
color='k',label='Modified Greenberg', linewidth=2, linestyle='-.',marker='+')
plt.plot(Eddiesmodel_Accuracy.Metrics, Eddiesmodel_Accuracy.Mean_absolute_error,
label='Eddies model', linewidth=2, dashes=dashList[3],marker='o',)
# plt.plot(ConstantExponential_Accuracy.Metrics, ConstantExponential_Accuracy.Mean_squared_error,
# label='Constant-Exponential')
plt.plot(LinearExponential_Accuracy.Metrics, LinearExponential_Accuracy.Mean_absolute_error,
label='Linear-Exponential',marker='v',)
# plt.plot(Linearlogarithm_Accuracy.Metrics, Linearlogarithm_Accuracy.Mean_squared_error,
# color='k', label='Linear-logarithm', linewidth=2, dashes=dashList[4],)
plt.legend(loc=2, fontsize=14)
plt.xlabel("Traffic speed $(mph)$"#, fontsize=18
)
plt.ylabel("Mean squared error"#, fontsize=18
)
plt.legend(loc=1, fontsize=12)
Out[148]:
In [149]:
#plt.style.use("seaborn-dark")
plt.style.use("seaborn-darkgrid")
In [150]:
f, ax = plt.subplots(figsize=(10, 7))
plt.plot(TwoRegime_Accuracy.Metrics, TwoRegime_Accuracy.Mean_squared_error,
label='Two-regime', linewidth=2, linestyle =':',marker='*')
plt.plot(ThreeRegime_Accuracy.Metrics, ThreeRegime_Accuracy.Mean_squared_error,
color='r',label='Three Regime', linewidth=2, dashes=dashList[1],marker='s')
plt.plot(ModifiedGreenberg_Accuracy.Metrics, ModifiedGreenberg_Accuracy.Mean_squared_error,
color='k',label='Modified Greenberg', linewidth=2, linestyle='-.',marker='+')
plt.plot(Eddiesmodel_Accuracy.Metrics, Eddiesmodel_Accuracy.Mean_squared_error,
label='Eddies model', linewidth=2, dashes=dashList[3],marker='o',)
# plt.plot(ConstantExponential_Accuracy.Metrics, ConstantExponential_Accuracy.Mean_squared_error,
# label='Constant-Exponential')
plt.plot(LinearExponential_Accuracy.Metrics, LinearExponential_Accuracy.Mean_squared_error,
label='Linear-Exponential',marker='v',)
# plt.plot(Linearlogarithm_Accuracy.Metrics, Linearlogarithm_Accuracy.Mean_squared_error,
# color='k', label='Linear-logarithm', linewidth=2, dashes=dashList[4],)
plt.legend(loc=2, fontsize=14)
plt.xlabel("Traffic speed $(mph)$"#, fontsize=18
)
plt.ylabel("Mean squared error"#, fontsize=18
)
plt.legend(loc=1, fontsize=12)
Out[150]:
Information Criterion¶
In [39]:
def log_lk(rng, params, model, *args, **kwargs):
model = substitute(seed(model, rng), params)
model_trace = trace(model).get_trace(*args, **kwargs)
obs_node = model_trace['obs']
#print(obs_node)
return obs_node['fn'].log_prob(obs_node['value']).sum()
def expected_log_likelihood(rng, params, model, *args, **kwargs):
p_waic = list(params.values())[0].shape[0]
log_lk_fn = vmap(lambda rng, params: log_lk(rng, params, model, *args, **kwargs))
log_lk_vals = log_lk_fn(random.split(rng, p_waic), params)
return log_lk_vals
def waic(log_like, return_lppd_p_waic=False):
"""Compute WAIC from a set of log-likelihoods."""
lppd = np.sum(scipy.special.logsumexp(log_like, axis=0, b=1/log_like.shape[0]))
p_waic = np.sum(np.var(log_like, axis=0))
if return_lppd_p_waic:
return onp.array(-2*(lppd - p_waic)).tolist(), onp.array(lppd).tolist(), onp.array(p_waic).tolist()
else:
return onp.array(-2*(lppd - p_waic)).tolist()
def log_likelihood(rng, params, model, *args, **kwargs):
n = list(params.values())[0].shape[0]
log_lk_fn = vmap(lambda rng, params: log_lk(rng, params, model, *args, **kwargs))
log_lk_vals = log_lk_fn(random.split(rng, n), params)
return logsumexp(log_lk_vals) - np.log(n)
In [40]:
from jax import random, vmap
from numpyro.handlers import sample, seed, substitute, trace
from jax.scipy.special import logsumexp
import scipy.special
In [42]:
rng, rng_ = random.split(random.PRNGKey(1))
TwoRegime = expected_log_likelihood(rng_,
samples_twoRegime,
model_twoRegime,
occupancy.values, speed.values)
In [43]:
rng, rng_ = random.split(random.PRNGKey(5))
ThreeRegime = expected_log_likelihood(rng_,
samples_threeRegime,
model_threeRegime,
occupancy.values, speed.values)
In [44]:
rng, rng_ = random.split(random.PRNGKey(10))
ConstantExponentiall = expected_log_likelihood(rng_,
ConstantExponentials,
ConstantExponential,
occupancy.values, speed.values)
In [45]:
rng, rng_ = random.split(random.PRNGKey(10))
LinearExponentiall = expected_log_likelihood(rng_,
LinearExponentialS,
LinearExponential,
occupancy.values, speed.values)
In [46]:
rng, rng_ = random.split(random.PRNGKey(5))
ModifiedGreenberg = expected_log_likelihood(rng_,
ModifiedGreenbergsamples,
ModifiedGreenbergmodel,
occupancy.values, speed.values)
In [47]:
rng, rng_ = random.split(random.PRNGKey(2))
Eddiesmodell = expected_log_likelihood(rng_,
Eddies_samples,
Eddiesmodel,
occupancy.values, speed.values)
In [48]:
rng, rng_ = random.split(random.PRNGKey(4))
Linearlogarithmm = expected_log_likelihood(rng_,
Linearlogarithmsamples,
Linearlogarithm,
occupancy.values, speed.values)
In [49]:
log_likes = {'Two - Regime': TwoRegime,
'Three - Regime':ThreeRegime,
'Constant - Exponential': ConstantExponentiall,
'Linear - Exponential': LinearExponentiall,
'Modified Greenberg': ModifiedGreenberg,
'Eddies model': Eddiesmodell,
'Linear logarithm': Linearlogarithmm,
}
In [50]:
waics = {key: onp.round(waic(val, return_lppd_p_waic=True), 1) for key, val in log_likes.items()}
In [51]:
model_compare = pd.DataFrame(waics.items(), columns=['Model','WAIC'])
outcompare = model_compare.merge(model_compare.WAIC.apply(pd.Series),
left_index = True, right_index = True).drop(["WAIC"], axis = 1)
outcompare.columns = ['Model', 'waic', 'lppd', 'p_waic']
outcompare = outcompare.sort_values(by="waic")
outcompare["dWAIC"] = outcompare["waic"] - outcompare.iloc[0, 1]
outcompare["weight"] = onp.round(torch.nn.functional.softmax(-1/2 * torch.tensor(outcompare["dWAIC"]),
dim=0), 1)
outcompare
Out[51]:
In [ ]:
# https://tonysyu.github.io/raw_content/matplotlib-style-gallery/gallery.html