Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SMC-ABC add distance, refactor and update notebook #3996

Merged
merged 10 commits into from
Jul 16, 2020
Prev Previous commit
Next Next commit
move dist functions out of simulator class
  • Loading branch information
aloctavodia committed Jul 6, 2020
commit fe3026074e0e3e4021a214c4150d65f4d65e9479
39 changes: 25 additions & 14 deletions pymc3/distributions/simulator.py
Original file line number Diff line number Diff line change
Expand Up @@ -59,12 +59,12 @@ def __init__(
self.epsilon = epsilon

if distance == "gaussian_kernel":
self.distance = self.gaussian_kernel
self.distance = gaussian_kernel
elif distance == "wasserstein":
self.distance = self.wasserstein
self.distance = wasserstein
sum_stat = "sort"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we do a assertion here or a warning that sum_stat need to be "sort" for wasserstein and energy?

elif distance == "energy":
self.distance = self.energy
self.distance = energy
sum_stat = "sort"
elif hasattr(distance, "__call__"):
self.distance = distance
Expand Down Expand Up @@ -107,17 +107,6 @@ def random(self, point=None, size=None):
else:
return np.array([self.function(*params) for _ in range(size)])

def gaussian_kernel(self, obs_data, sim_data):
return np.sum(-0.5 * ((obs_data - sim_data) / self.epsilon) ** 2)

# we are assuming obs_data and sim_data are already sorted!
def wasserstein(self, obs_data, sim_data):
return np.mean(np.abs((obs_data - sim_data) / self.epsilon))

# we are assuming obs_data and sim_data are already sorted!
def energy(self, obs_data, sim_data):
return 1.4142 * np.mean(((obs_data - sim_data) / self.epsilon) ** 2) ** 0.5

def _repr_latex_(self, name=None, dist=None):
if dist is None:
dist = self
Expand All @@ -130,4 +119,26 @@ def _repr_latex_(self, name=None, dist=None):


def identity(x):
"""identity function, used as a summary statistics"""
return x


def gaussian_kernel(epsilon, obs_data, sim_data):
"""gaussian distance function"""
return np.sum(-0.5 * ((obs_data - sim_data) / epsilon) ** 2)


def wasserstein(epsilon, obs_data, sim_data):
"""wasserstein distance function

We are assuming obs_data and sim_data are already sorted!
"""
return np.mean(np.abs((obs_data - sim_data) / epsilon))


def energy(epsilon, obs_data, sim_data):
"""Energy distance function

We are assuming obs_data and sim_data are already sorted!
"""
return 1.4142 * np.mean(((obs_data - sim_data) / epsilon) ** 2) ** 0.5
2 changes: 1 addition & 1 deletion pymc3/smc/smc.py
Original file line number Diff line number Diff line change
Expand Up @@ -366,4 +366,4 @@ def __call__(self, posterior):
sim_data = self.function(**func_parameters)
if self.save:
self.save_data(sim_data)
return self.distance(self.observations, self.sum_stat(sim_data))
return self.distance(self.epsilon, self.observations, self.sum_stat(sim_data))
4 changes: 2 additions & 2 deletions pymc3/tests/test_smc.py
Original file line number Diff line number Diff line change
Expand Up @@ -104,8 +104,8 @@ def normal_sim(a, b):
return np.random.normal(a, b, 1000)

with pm.Model() as self.SMABC_test:
a = pm.Normal("a", mu=0, sd=5)
b = pm.HalfNormal("b", sd=2)
a = pm.Normal("a", mu=0, sigma=5)
b = pm.HalfNormal("b", sigma=2)
s = pm.Simulator(
"s", normal_sim, params=(a, b), sum_stat="sort", epsilon=1, observed=self.data
)
Expand Down