diff --git a/examples/multiple.py b/examples/multiple.py
new file mode 100644
index 00000000..f2f76881
--- /dev/null
+++ b/examples/multiple.py
@@ -0,0 +1,52 @@
+import itertools
+
+import numpy as np
+import physbo
+
+# Make a set of candidates, test_X
+D = 2  # The number of params (the dimension of parameter space)
+Nx = 11  # The number of candidates
+N = Nx*Nx
+
+# score = "HVPI"
+score = "EHVI"
+
+a = np.linspace(-2, 2, Nx)
+test_X = np.array(list(itertools.product(a, a)))
+
+
+def vlmop2_minus(x):
+    n = x.shape[1]
+    y1 = 1 - np.exp(-1 * np.sum((x - 1 / np.sqrt(n)) ** 2, axis=1))
+    y2 = 1 - np.exp(-1 * np.sum((x + 1 / np.sqrt(n)) ** 2, axis=1))
+
+    return np.c_[-y1, -y2]
+
+
+class simulator(object):
+    def __init__(self, X):
+        self.t = vlmop2_minus(X)
+
+    def __call__(self, action):
+        return self.t[action]
+
+
+sim = simulator(test_X)
+
+policy = physbo.search.discrete_multi.policy(test_X, num_objectives=2)
+policy.set_seed(0)
+# Random search (10 times)
+policy.random_search(max_num_probes=10, simulator=sim)
+
+# Bayesian search (100 times)
+#   score function (acquition function): expectation of improvement (EI)
+policy.bayes_search(max_num_probes=0, simulator=sim, score=score, interval=0)
+
+print("Pareto fronts:")
+res = policy.history
+front, front_index = res.export_pareto_front()
+for fr, ifr in zip(front, front_index):
+    print("  action: ", ifr)
+    print("  X: ", test_X[ifr, :])
+    print("  f: ", fr)
+    print()
diff --git a/examples/multiple_score.py b/examples/multiple_score.py
new file mode 100644
index 00000000..ec9ecdfe
--- /dev/null
+++ b/examples/multiple_score.py
@@ -0,0 +1,64 @@
+import itertools
+
+import numpy as np
+import physbo
+
+# Make a set of candidates, test_X
+D = 2  # The number of params (the dimension of parameter space)
+Nx = 11  # The number of candidates
+N = Nx*Nx
+
+# score = "HVPI"
+score = "EHVI"
+
+a = np.linspace(-2, 2, Nx)
+test_X = np.array(list(itertools.product(a, a)))
+
+
+def vlmop2_minus(x):
+    n = x.shape[1]
+    y1 = 1 - np.exp(-1 * np.sum((x - 1 / np.sqrt(n)) ** 2, axis=1))
+    y2 = 1 - np.exp(-1 * np.sum((x + 1 / np.sqrt(n)) ** 2, axis=1))
+
+    return np.c_[-y1, -y2]
+
+
+class simulator(object):
+    def __init__(self, X):
+        self.t = vlmop2_minus(X)
+
+    def __call__(self, action):
+        return self.t[action]
+
+
+sim = simulator(test_X)
+
+# actions = np.arange(N)
+# np.random.shuffle(actions)
+# n = 10
+# actions = actions[0:n]
+# data = sim(actions)
+# policy = physbo.search.discrete_multi.policy(test_X, num_objectives=2, initial_data=[actions, data])
+
+policy = physbo.search.discrete_multi.policy(test_X, num_objectives=2)
+policy.set_seed(0)
+# Random search (10 times)
+policy.random_search(max_num_probes=10, simulator=sim)
+
+# Bayesian search (100 times)
+#   score function (acquisition function): expectation of improvement (EI)
+policy.bayes_search(max_num_probes=0, simulator=sim, score=score, interval=0)
+
+print("Mean values of prediction")
+scores = policy.get_post_fmean(xs=test_X)
+print(scores)
+print()
+
+print("Standard derivations of prediction")
+scores = policy.get_post_fcov(xs=test_X)
+print(np.sqrt(scores))
+print()
+
+print("Acquisition function")
+scores = policy.get_score(mode=score, xs=test_X)
+print(scores)
diff --git a/examples/simple_score.py b/examples/simple_score.py
new file mode 100644
index 00000000..713ea0cb
--- /dev/null
+++ b/examples/simple_score.py
@@ -0,0 +1,41 @@
+import numpy as np
+import physbo
+
+# Make a set of candidates, test_X
+D = 3  # The number of params (the dimension of parameter space)
+N = 100  # The number of candidates
+test_X = np.random.randn(N, D)  # Generated from Gaussian
+score = "EI"
+
+
+def simulator(actions: np.ndarray) -> np.ndarray:
+    """Objective function
+
+    Quadratic function, -Σ_i x_i^2
+    Recieves an array of actions (indices of candidates) and returns the corresponding results as an array
+    """
+    return -np.sum(test_X[actions, :] ** 2, axis=1)
+
+
+policy = physbo.search.discrete.policy(test_X)
+
+# Random search (10 times)
+policy.random_search(max_num_probes=10, simulator=simulator)
+
+# Bayesian search (100 times)
+#   score function (acquisition function): expectation of improvement (EI)
+policy.bayes_search(max_num_probes=0, simulator=simulator, score=score)
+
+print("Mean values of prediction")
+scores = policy.get_post_fmean(xs=test_X)
+print(scores)
+print()
+
+print("Standard derivations of prediction")
+scores = policy.get_post_fcov(xs=test_X)
+print(np.sqrt(scores))
+print()
+
+print("Acquisition function")
+scores = policy.get_score(mode=score, xs=test_X)
+print(scores)