Optimize a Mystery Function with BO

Use Bayesian optimization to find function optima faster than random search.

Introduction

30 Second Summary

Imagine tuning the settings on a machine learning model, but each experiment takes minutes to run. You can't afford to try hundreds of random combinations hoping to stumble on the best one.

In this project, you will build a Python tool that uses Bayesian optimization to find optimal solutions with far fewer attempts than random guessing. You'll go from optimizing a tricky math function to tuning a real machine learning model's hyperparameters.

What You'll Build

You'll run a single command, watch the optimizer intelligently zero in on the best answer to a mystery function, and see a side-by-side convergence plot proving it beat random search with the same budget of attempts.

By the end of this project, you'll have:

  • A convergence comparison plot proving Bayesian optimization finds a better optimum than random search using the same 25 evaluation budget.
  • Publication-quality visualizations of the Gaussian Process surrogate model, uncertainty bands, and acquisition function revealing exactly how the optimizer decides where to sample next.
  • A real-world tuning script that optimizes a Random Forest classifier on the Wine dataset and prints the accuracy improvement over random search.
  • Secret Mission: Experiment with different exploration strategies by tuning the acquisition function's kappa parameter and visualizing how the exploration/exploitation tradeoff changes the optimizer's behavior.

Are there any prerequisites?

You'll need basic Python knowledge (writing functions, running scripts from the terminal). Familiarity with what a machine learning model is and what hyperparameters are will help, but you don't need to be an ML expert.

Before We Start

Before we dive in, let's lock in what you're about to build and why it matters.

Set Up Python and Your Project

Every experiment in this project runs in Python and relies on specialized optimization libraries. Before you can watch random search struggle or Bayesian optimization succeed, you need the right version of Python and the right packages ready to go.

In this step, you'll install Python 3.13, create an isolated project environment, and pull in the libraries that power Bayesian optimization.

In this step, get ready to:
  • Install Python 3.13 and verify it's on your PATH.
  • Create a virtual environment and install the optimization libraries.
  • Set up version control with git.
Check your Python version

Python 3.13 is the version this project targets. It includes the scientific libraries we need and has stable support across all major platforms.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type Terminal and press Enter to open it.
  • Check if Python 3.13 is already installed by running the version check for your platform:

macOS

python3 --version

Windows

python --version

Check your output and pick the tab that matches what you see:

✔️ I see Python 3.13.x

Your Python version is correct. Skip ahead to the next substep.

ⓧ I see an older version

You have Python installed, but it's an older version. The optimization libraries in this project require 3.13 or higher.

  • Go to https://www.python.org/downloads/ in your browser.
  • Download the Python 3.13 installer for your platform (macOS 64-bit universal2 installer for Mac, Windows installer 64-bit for Windows).
  • Run the installer.

Important installer settings

On Windows, check the box that says Add python.exe to PATH on the very first screen of the installer. Without this, your terminal will not find the python command.

On macOS, the installer adds Python to your PATH automatically.

  • After installation completes, close and reopen your terminal.
  • Run the version check again to confirm you see Python 3.13.x in the output.

ⓧ Command not found

Python is not installed on your system yet.

  • Go to https://www.python.org/downloads/ in your browser.
  • Download the Python 3.13 installer for your platform (macOS 64-bit universal2 installer for Mac, Windows installer 64-bit for Windows).
  • Run the installer.

Important installer settings

On Windows, check the box that says Add python.exe to PATH on the very first screen of the installer. Without this, your terminal will not find the python command.

On macOS, the installer adds Python to your PATH automatically.

  • After installation completes, close and reopen your terminal.
  • Run the version check again to confirm you see Python 3.13.x in the output.
Create your project and install libraries

With Python confirmed, you'll create a dedicated folder for this project and install the libraries inside an isolated virtual environment.

  • Move to your Desktop by running this command:
cd ~/Desktop
  • Create your project folder and move into it by running these commands:
mkdir bayesian-opt
cd bayesian-opt
  • Create a virtual environment by running the command for your platform:

macOS

python3 -m venv venv

Windows

python -m venv venv
  • Confirm the virtual environment was created by running:
ls

You should see a venv folder listed in the output.

What is a virtual environment?

A virtual environment is an isolated Python installation inside your project folder. Packages you install here stay separate from the rest of your system, preventing version conflicts between projects.

  • Activate the virtual environment by running the command for your platform:

macOS

source venv/bin/activate

Windows

venv\Scripts\activate

Your terminal prompt now shows (venv) at the beginning, confirming the environment is active.

  • Install the optimization and plotting libraries by running:
pip install bayesian-optimization==3.3.0 matplotlib

You'll see pip downloading and installing several packages including numpy and scikit-learn (they're automatic dependencies).

  • Before you verify, what packages do you expect to see when you run pip list?
  • Verify the installation by running:
pip list

You should see bayesian-optimization and matplotlib in the list, along with their dependencies like numpy and scikit-learn.

Don't see the packages?

Make sure your virtual environment is active. Your terminal prompt should start with (venv). If it doesn't, run the activate command again for your platform.

If pip install failed with a network error, check your internet connection and try again. Help me troubleshoot my pip install

Set up version control

A requirements.txt file records your exact dependencies so anyone can recreate your environment. A .gitignore keeps generated files out of your repository.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type VS Code and press Enter to open it.
  • Click File in the top menu bar.
  • Click Open Folder.
  • Select the bayesian-opt folder on your Desktop, then click Open.
  • Create a new file by clicking the new file icon in the Explorer sidebar (or press Cmd+N on macOS / Ctrl+N on Windows).
  • Name the file requirements.txt and paste the following content:
bayesian-optimization==3.3.0
matplotlib

What does this file do?

This file pins the exact versions of your project's direct dependencies. Anyone who clones your repository can run pip install -r requirements.txt to get the same setup.

  • Save requirements.txt (Cmd+S on macOS / Ctrl+S on Windows).
  • Create another new file named .gitignore and paste the following content:
venv/
__pycache__/
*.pyc
*.png
.DS_Store

What does this file do?

This tells Git to ignore your virtual environment folder, compiled Python bytecode, generated plot images, and macOS system files. These are all files that can be regenerated and don't belong in a repository.

  • Save .gitignore.
  • Open the integrated terminal in VS Code by pressing Ctrl+` (backtick key, same shortcut on both platforms).
  • Initialize a git repository by running:
git init

What do you expect git status to show right now?

  • Verify your repository by running:
git status

You should see requirements.txt and .gitignore listed as untracked files. The venv/ folder should NOT appear because the .gitignore is already filtering it out.

Seeing venv/ in the untracked files?

Double-check that your .gitignore file contains venv/ on its first line and that you saved the file before running git status. Help me fix my .gitignore

✔️ Awesome, I've got everything!

Great. Double-check you've saved both files.

ⓧ I'd like to double check the full code

Here are both files in full.

bayesian-optimization==3.3.0
matplotlib
venv/
__pycache__/
*.pyc
*.png
.DS_Store

Your environment is ready. Next up, you'll write a random search script and see firsthand why brute-force optimization wastes evaluations.

Try Random Search and Watch It Struggle

Your environment is set up and ready to go. Now it's time to confront the problem that motivates this entire project.

Imagine you have a function that's expensive to evaluate and you want to find its maximum. The naive approach is random search: pick points at random and hope one lands near the optimum. In this step, you'll watch that approach struggle on a tricky function with multiple peaks.

In this step, get ready to:
  • Create a mystery function with multiple peaks that's hard to optimize.
  • Run random search and observe how it misses the true optimum.
  • Visualize random search's convergence with a matplotlib plot.
Define the mystery function and random search logic

The function you're optimizing has multiple peaks and valleys. Random search has no way to learn from previous samples, so it just guesses blindly across the entire range.

  • Create a new file called random_search.py in your project folder.
  • Add the imports and the two core functions by pasting this code:
"""Random search baseline: optimize a tricky function by sampling blindly."""

import random
import math
import matplotlib.pyplot as plt


def mystery_function(x):
    """A tricky function with multiple peaks. Pretend we don't know this formula."""
    return math.sin(x) * math.exp(-0.1 * x) + math.cos(3 * x) * 0.5


def run_random_search(n_evaluations, seed=None):
    """Run random search and return history of (x, target) pairs."""
    if seed is not None:
        random.seed(seed)

    history = []
    for _ in range(n_evaluations):
        x = random.uniform(0, 10)
        target = mystery_function(x)
        history.append((x, target))

    return history

What does this code do?

You've defined two functions that form the foundation of the experiment.

  • mystery_function(x) combines a decaying sine wave with a cosine wave. This creates multiple peaks across the range [0, 10], with the global maximum near x = 1.0.
  • run_random_search() picks random x values between 0 and 10, evaluates the function at each one, and returns the full history. It has no strategy and no memory of past results.
Add the main function to run the search

The main function calls the random search and tracks the best result found at each evaluation. This lets you see how (or whether) random search improves over time.

  • Add this code below the run_random_search function in random_search.py:
def main():
    n_evaluations = 25

    print("=" * 50)
    print(f"RANDOM SEARCH ({n_evaluations} evaluations)")
    print("=" * 50)

    history = run_random_search(n_evaluations, seed=42)

    best_target = float("-inf")
    best_x = None
    best_so_far = []

    for i, (x, target) in enumerate(history):
        if target > best_target:
            best_target = target
            best_x = x
            print(f"  New best at eval {i + 1}: x = {x:.4f}, target = {target:.4f}")
        best_so_far.append(best_target)

    print(f"\nFinal result after {n_evaluations} evaluations:")
    print(f"  Best x: {best_x:.4f}")
    print(f"  Best target: {best_target:.4f}")
    print(f"  (True optimum is approximately 1.13 at x ~ 1.0)")


if __name__ == "__main__":
    main()

What does this code do?

  • The loop iterates through all 25 random evaluations and prints whenever a new best value is found.
  • best_so_far records the running maximum at each step. You'll use this list to plot convergence shortly.
  • The final print compares what random search found against the true optimum (approximately 1.13 at x ~ 1.0) so you can see the gap.
  • Save random_search.py.
Run the script and watch random search struggle

The function has its true optimum near x = 1.0, with a target value of approximately 1.13. Before you run the script, think about this: do you expect 25 random guesses spread across a range of 0 to 10 to land close to that peak?

  • Run the script by entering this command in your terminal:
python random_search.py

You should see output like this, showing each time random search stumbles onto a better value:

==================================================
RANDOM SEARCH (25 evaluations)
==================================================
  New best at eval 1: x = 6.3946, target = -0.2863
  New best at eval 2: x = 0.2861, target = 0.7578
  New best at eval 5: x = 1.5602, target = 0.8650

Final result after 25 evaluations:
  Best x: 1.5602
  Best target: 0.8650
  (True optimum is approximately 1.13 at x ~ 1.0)

Notice the gap. The best target random search found (around 0.87) is noticeably below the true optimum of 1.13. It also found only 3 improvements across 25 evaluations. Most guesses were wasted.

Why does random search struggle here?

The mystery function has multiple peaks spread across the range [0, 10]. The global maximum occupies a narrow region near x = 1.0.

Random search picks points uniformly across the full range. It has no memory of where it has already looked and no model of which regions look promising. Each guess is independent of all previous ones.

If you change the seed value from 42 to any other number and re-run, you'll get completely different results. Some seeds get lucky and land closer to the optimum. Others do worse. This inconsistency is the core weakness of random search.

Seeing a ModuleNotFoundError?

  • Make sure your virtual environment is active. You should see (venv) at the start of your terminal prompt.
  • If it's not active, run venv\Scripts\activate (Windows) or source venv/bin/activate (macOS/Linux) and try again.

Why is my random_search.py throwing a ModuleNotFoundError?

Add a convergence plot to visualize the waste

Numbers in the terminal tell part of the story. A convergence plot makes the waste obvious: you'll see the best-found value flatline early because random search runs out of luck.

  • In random_search.py, find this line at the end of the main() function:
    print(f"  (True optimum is approximately 1.13 at x ~ 1.0)")
  • Add this plotting code directly below that line (still inside the main() function, above the if __name__ block):
    # Plot convergence
    plt.figure(figsize=(10, 5))
    plt.plot(range(1, n_evaluations + 1), best_so_far, "ro-", markersize=4)
    plt.xlabel("Number of Evaluations")
    plt.ylabel("Best Target Found")
    plt.title("Random Search Convergence")
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig("random_search_convergence.png", dpi=150)
    plt.show()
    print("\nSaved: random_search_convergence.png")

What does this code do?

  • plt.plot() draws the best-found-so-far value at each evaluation number. This is a convergence curve.
  • plt.savefig() saves the plot as a PNG file so you can reference it later without re-running.
  • A flat convergence curve means random search stopped improving. Every evaluation after that flatline was wasted compute.
  • Save random_search.py.
  • Run the script again:
python random_search.py

A matplotlib window pops up showing the convergence plot. You should see the curve jump up in the first few evaluations, then flatline for the remaining ones. Close the plot window to let the script finish.

  • Confirm that random_search_convergence.png was saved by running:
ls

You should see random_search_convergence.png listed alongside random_search.py and your other project files.

Plot window not appearing?

  • If you're running in a headless environment or over SSH, plt.show() may not display. The PNG file is still saved. Open random_search_convergence.png directly from your file explorer to see the plot.
  • If you see a TclError or display-related error, add import matplotlib; matplotlib.use('Agg') at the very top of the file before import matplotlib.pyplot as plt. This disables the interactive window but still saves the PNG.

Help me fix my matplotlib plot display issue in random_search.py.

✔️ Awesome, I've got everything!

Great. Double check that you've saved random_search.py and that random_search_convergence.png exists in your project folder.

ⓧ I'd like to double check the full code

Here's the complete random_search.py file. Compare it against yours:

"""Random search baseline: optimize a tricky function by sampling blindly."""

import random
import math
import matplotlib.pyplot as plt


def mystery_function(x):
    """A tricky function with multiple peaks. Pretend we don't know this formula."""
    return math.sin(x) * math.exp(-0.1 * x) + math.cos(3 * x) * 0.5


def run_random_search(n_evaluations, seed=None):
    """Run random search and return history of (x, target) pairs."""
    if seed is not None:
        random.seed(seed)

    history = []
    for _ in range(n_evaluations):
        x = random.uniform(0, 10)
        target = mystery_function(x)
        history.append((x, target))

    return history
def main():
    n_evaluations = 25

    print("=" * 50)
    print(f"RANDOM SEARCH ({n_evaluations} evaluations)")
    print("=" * 50)

    history = run_random_search(n_evaluations, seed=42)

    best_target = float("-inf")
    best_x = None
    best_so_far = []

    for i, (x, target) in enumerate(history):
        if target > best_target:
            best_target = target
            best_x = x
            print(f"  New best at eval {i + 1}: x = {x:.4f}, target = {target:.4f}")
        best_so_far.append(best_target)

    print(f"\nFinal result after {n_evaluations} evaluations:")
    print(f"  Best x: {best_x:.4f}")
    print(f"  Best target: {best_target:.4f}")
    print(f"  (True optimum is approximately 1.13 at x ~ 1.0)")

    # Plot convergence
    plt.figure(figsize=(10, 5))
    plt.plot(range(1, n_evaluations + 1), best_so_far, "ro-", markersize=4)
    plt.xlabel("Number of Evaluations")
    plt.ylabel("Best Target Found")
    plt.title("Random Search Convergence")
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig("random_search_convergence.png", dpi=150)
    plt.show()
    print("\nSaved: random_search_convergence.png")


if __name__ == "__main__":
    main()

You've now experienced the fundamental problem with random search: it wastes evaluations because it can't learn from past results. In the next step, you'll replace it with Bayesian optimization, which builds a model of the function and intelligently decides where to look next.

Replace Random Search with Bayesian Optimization

You saw random search struggle in the last step. It sampled 25 points blindly and often landed far from the true optimum because it has no memory and no strategy.

Bayesian optimization takes a fundamentally different approach. After each evaluation, it builds a model of the function and uses that model to decide where to look next. In this step, you'll replace random search with BO and see it find a better result with the same 25-evaluation budget.

In this step, get ready to:
  • Create a Bayesian optimization script targeting the same mystery function.
  • Run a side-by-side comparison with random search using the same evaluation budget.
  • Generate a convergence plot showing how quickly each method finds the optimum.
Set up the optimizer script

The new script needs the same mystery function from before, plus a streamlined random search helper for comparison. It also imports the BayesianOptimization class from the library you installed in Step 1.

  • Create a new file called bayesian_opt.py in your project folder and add these imports and helper functions:
"""Bayesian optimization vs random search on the same mystery function."""

import random
import math
import matplotlib.pyplot as plt
from bayes_opt import BayesianOptimization


def mystery_function(x):
    """A tricky function with multiple peaks. Pretend we don't know this formula."""
    return math.sin(x) * math.exp(-0.1 * x) + math.cos(3 * x) * 0.5


def run_random_search(n_evaluations, seed=42):
    """Run random search and return best-so-far history."""
    random.seed(seed)
    best_target = float("-inf")
    best_so_far = []

    for _ in range(n_evaluations):
        x = random.uniform(0, 10)
        target = mystery_function(x)
        if target > best_target:
            best_target = target
        best_so_far.append(best_target)

    return best_so_far

What does this code do?

  • from bayes_opt import BayesianOptimization imports the optimizer class from the bayesian-optimization library.
  • mystery_function is the same tricky multi-peak function from your random search script.
  • run_random_search is a simplified version that returns only the best-so-far list (instead of all x/target pairs). This makes the comparison plot easier to build later.

Now add the main() function that runs the Bayesian optimizer. This is the core of the script.

  • Below the run_random_search function, add this code:
def main():
    n_evaluations = 25

    # --- Bayesian Optimization ---
    print("=" * 50)
    print(f"BAYESIAN OPTIMIZATION ({n_evaluations} evaluations)")
    print("=" * 50)

    optimizer = BayesianOptimization(
        f=mystery_function,
        pbounds={"x": (0, 10)},
        random_state=42,
        verbose=2,
    )

    optimizer.maximize(
        init_points=5,
        n_iter=20,
    )

    bo_best = optimizer.max
    bo_best_so_far = []
    running_best = float("-inf")
    for res in optimizer.res:
        if res["target"] > running_best:
            running_best = res["target"]
        bo_best_so_far.append(running_best)


if __name__ == "__main__":
    main()

What does this code do?

  • BayesianOptimization(f=mystery_function, pbounds={"x": (0, 10)}) tells the optimizer what function to maximize and over what range.
  • init_points=5 means 5 random samples to start (building an initial picture of the function).
  • n_iter=20 means 20 guided samples after that. Total: 25 evaluations, the same budget as random search.
  • verbose=2 prints a table showing each iteration's target and parameters as the optimizer runs.
  • The loop over optimizer.res reconstructs the best-so-far curve for plotting later.

Why this library instead of Optuna?

Optuna is a popular hyperparameter tuning framework, but it hides the Gaussian Process internals behind a high-level API. The bayesian-optimization library exposes the surrogate model, acquisition function, and sampling decisions directly. That transparency is essential for understanding how BO actually works.

  • Save bayesian_opt.py.
  • Run the script to see the optimizer in action:
python bayesian_opt.py

You should see a table printed with columns for iteration number, target value, and the x parameter. Each row shows where the optimizer chose to sample and what value it found there.

Seeing ModuleNotFoundError: No module named 'bayes_opt'?

  • Make sure your virtual environment is active. You should see (venv) at the start of your terminal prompt.
  • If not, activate it again with venv\Scripts\activate (Windows) or source venv/bin/activate (macOS/Linux), then run the script again.

Still stuck? Help me fix this ModuleNotFoundError for bayes_opt

Add a side-by-side comparison

The iteration table shows BO working, but you need a direct comparison with random search to prove it actually does better with the same budget.

  • In bayesian_opt.py, find the line that says bo_best_so_far.append(running_best) at the end of the for loop.
  • After the for loop (but still inside main()), add this comparison code:
    # --- Random Search ---
    print("\n" + "=" * 50)
    print(f"RANDOM SEARCH ({n_evaluations} evaluations)")
    print("=" * 50)

    random_best_so_far = run_random_search(n_evaluations)

    # --- Comparison ---
    print("\n" + "=" * 50)
    print("COMPARISON")
    print("=" * 50)
    print(f"  Random Search best: {random_best_so_far[-1]:.4f}")
    print(f"  Bayesian Opt best:  {bo_best['target']:.4f}")
    print(f"  BO found x = {bo_best['params']['x']:.4f}")
    print(f"  Improvement: {bo_best['target'] - random_best_so_far[-1]:.4f}")

What does this code do?

  • run_random_search(n_evaluations) runs the same 25-evaluation random search from Step 2, returning its best-so-far curve.
  • The COMPARISON section prints both results side by side. optimizer.max gives you the best target and the x value where it was found.
  • The improvement line shows exactly how much better BO did with the same budget.
  • Save bayesian_opt.py and run it again:
python bayesian_opt.py

After the iteration table, you should now see a COMPARISON section showing BO's best target is higher than random search's best target. The improvement value should be positive.

Seeing a negative improvement?

  • Check that the comparison code is placed after the for loop that builds bo_best_so_far, not inside it.
  • Make sure random_state=42 is set in the BayesianOptimization constructor for reproducible results.

Need help? Help me debug my BO comparison output

Generate the convergence comparison plot

Numbers prove BO wins, but a convergence plot shows how it wins. You can see BO pulling ahead early because its guided samples home in on the optimum while random search plateaus.

  • In bayesian_opt.py, add this plotting code at the end of main() (after the Improvement print line, still inside the function):
    # --- Convergence comparison plot ---
    plt.figure(figsize=(10, 5))
    evals = range(1, n_evaluations + 1)
    plt.plot(evals, random_best_so_far, "ro-", markersize=4, label="Random Search")
    plt.plot(evals, bo_best_so_far, "b^-", markersize=4, label="Bayesian Optimization")
    plt.xlabel("Number of Evaluations")
    plt.ylabel("Best Target Found")
    plt.title("Convergence: Random Search vs Bayesian Optimization")
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig("convergence_comparison.png", dpi=150)
    plt.show()
    print("\nSaved: convergence_comparison.png")

What does this code do?

  • Both convergence curves are plotted on the same axes so you can visually compare how quickly each method improves.
  • Random search uses red circles ("ro-") and BO uses blue triangles ("b^-") so you can tell them apart at a glance.
  • The plot is saved as convergence_comparison.png in your project folder.
  • Save bayesian_opt.py.

Before you run the final version, think about what the convergence plot will look like. Do you expect the BO curve to pull ahead of random search early, or only at the very end?

  • Run the complete script:
python bayesian_opt.py

You should see the iteration table, the COMPARISON section showing a positive improvement, and a matplotlib window displaying both convergence curves. The BO curve (blue triangles) should climb higher and faster than the random search curve (red circles). After closing the plot window, you should see the message confirming convergence_comparison.png was saved.

Plot window not appearing?

  • If you're running via SSH or a headless environment, plt.show() may not work. The PNG file is still saved to disk. Open convergence_comparison.png directly from your file explorer.
  • Make sure the plotting code is inside the main() function (indented by 4 spaces), not at the top level of the file.

Still stuck? Help me fix my matplotlib plot in bayesian_opt.py

✔️ Awesome, I've got everything!

Great. Double check you've saved bayesian_opt.py and confirmed the plot file exists.

ⓧ I'd like to double check the full code

Scroll up and confirm your bayesian_opt.py contains each of these sections in order:

  • The imports (including from bayes_opt import BayesianOptimization) and both helper functions (mystery_function and run_random_search).
  • The main() function starting with the Bayesian Optimization section (optimizer creation and maximize call, then the bo_best_so_far loop).
  • The Random Search and COMPARISON print sections inside main().
  • The convergence comparison plot code at the end of main().
  • The if __name__ == "__main__": block at the very bottom.

Your convergence plot proves that BO consistently outperforms random search with the same budget. Next up, you'll look under the hood and visualize the Gaussian Process model that makes BO so much smarter about where to sample.

Visualize the Gaussian Process and Acquisition Function

Your convergence plot proves that Bayesian optimization finds better optima. But it does not explain why.

In this step, you will crack open the black box. You will visualize the Gaussian Process surrogate model that BO builds internally and the acquisition function it uses to decide where to sample next. These two plots reveal the mechanism: the model predicts where good values might hide, and uncertainty tells it where to explore.

In this step, get ready to:
  • Fit a Gaussian Process regressor to BO's observed points.
  • Plot the GP surrogate model with confidence bands.
  • Plot the Upper Confidence Bound acquisition function and its next suggested sample.
Run BO and fit a Gaussian Process to the observations

The bayesian-optimization library builds a GP internally but does not expose it for plotting. You will reconstruct the surrogate yourself using scikit-learn's GaussianProcessRegressor with the same Matern 5/2 kernel the library uses.

  • Create a new file called visualize_bo.py in your project folder and paste the following code:
"""Visualize the Gaussian Process surrogate and acquisition function."""

import math
import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern
from bayes_opt import BayesianOptimization


def mystery_function(x):
    """A tricky function with multiple peaks."""
    return math.sin(x) * math.exp(-0.1 * x) + math.cos(3 * x) * 0.5


def main():
    # Run BO to collect observations
    optimizer = BayesianOptimization(
        f=mystery_function,
        pbounds={"x": (0, 10)},
        random_state=42,
        verbose=0,
    )
    optimizer.maximize(init_points=3, n_iter=7)

    # Extract observed points
    x_observed = np.array([[res["params"]["x"]] for res in optimizer.res])
    y_observed = np.array([res["target"] for res in optimizer.res])

What does this code do?

This runs a short BO session (3 initial random points + 7 guided iterations = 10 total observations) with verbose=0 so no iteration table is printed.

After the run, you extract the x-values and targets from optimizer.res into NumPy arrays. These are the observations you will use to reconstruct the surrogate model.

  • Directly below the code you just pasted, add the GP fitting and prediction code:
    # Fit a GP to the observed data (same kernel family the library uses)
    gp = GaussianProcessRegressor(
        kernel=Matern(nu=2.5),
        alpha=1e-6,
        normalize_y=True,
        n_restarts_optimizer=5,
    )
    gp.fit(x_observed, y_observed)

    # Create a dense grid for plotting
    x_grid = np.linspace(0, 10, 300).reshape(-1, 1)
    y_true = np.array([mystery_function(x[0]) for x in x_grid])
    y_pred, y_std = gp.predict(x_grid, return_std=True)

What does this code do?

  • GaussianProcessRegressor with Matern(nu=2.5) fits a smooth curve through the 10 observed points, producing both a mean prediction and an uncertainty estimate at every location.
  • x_grid creates 300 evenly spaced test points across [0, 10] so the plot looks smooth.
  • gp.predict(x_grid, return_std=True) returns both the GP mean (y_pred) and standard deviation (y_std) at each grid point. You will use these for the confidence bands.
Plot the GP surrogate model

The GP mean shows what the model thinks the function looks like. The confidence interval shows where it is uncertain. Where the bands are wide, the model has not explored yet.

  • Below the prediction code, add the first plot and save it as gp_surrogate.png:
    # --- Plot 1: GP Surrogate Model ---
    fig, ax = plt.subplots(figsize=(12, 6))
    ax.plot(x_grid, y_true, "k--", linewidth=1.5, label="True Function (unknown)")
    ax.plot(x_grid, y_pred, "b-", linewidth=2, label="GP Mean Prediction")
    ax.fill_between(
        x_grid.flatten(),
        y_pred - 1.96 * y_std,
        y_pred + 1.96 * y_std,
        alpha=0.2,
        color="blue",
        label="95% Confidence Interval",
    )
    ax.scatter(
        x_observed, y_observed, c="red", s=80, zorder=5, label="Observed Points"
    )
    ax.set_xlabel("x", fontsize=12)
    ax.set_ylabel("f(x)", fontsize=12)
    ax.set_title("Gaussian Process Surrogate Model", fontsize=14)
    ax.legend(fontsize=11)
    ax.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig("gp_surrogate.png", dpi=150)
    plt.show()
    print("Saved: gp_surrogate.png")

What does this plot show?

  • The dashed black line is the true mystery function (which in real problems you would not know).
  • The solid blue line is the GP's best guess based on only 10 observations.
  • The shaded blue region is the 95% confidence interval. Where the band is narrow, the model is confident. Where it is wide, the model is uncertain.
  • Red dots are the actual points BO sampled.
  • Save visualize_bo.py.
  • Run the script to confirm the first plot saves correctly:
python visualize_bo.py

A matplotlib window pops up showing the GP surrogate. After you close it, you should see Saved: gp_surrogate.png printed in your terminal.

Not seeing the plot window?

  • If the script hangs without showing a window, your matplotlib backend may not support interactive display. Close the process with Ctrl+C and check that gp_surrogate.png was still saved in your project folder. The file saves before plt.show() is called.
  • If you see an import error for sklearn, confirm your virtual environment is active. Scikit-learn is installed automatically as a dependency of bayesian-optimization.

Why isn't my GP surrogate plot rendering?

Plot the Upper Confidence Bound acquisition function

The acquisition function answers the question: given what the GP predicts, where should we sample next? The Upper Confidence Bound (UCB) formula is simply the GP mean plus 1.96 times the standard deviation. It peaks where the model expects high values OR where uncertainty is large.

  • Below the first plot code in visualize_bo.py, add the two-panel acquisition function figure:
    # --- Plot 2: GP + Acquisition Function ---
    ucb = y_pred + 1.96 * y_std  # Upper Confidence Bound

    fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)

    # Top: GP model
    axes[0].plot(x_grid, y_true, "k--", linewidth=1.5, label="True Function (unknown)")
    axes[0].plot(x_grid, y_pred, "b-", linewidth=2, label="GP Mean")
    axes[0].fill_between(
        x_grid.flatten(),
        y_pred - 1.96 * y_std,
        y_pred + 1.96 * y_std,
        alpha=0.2,
        color="blue",
        label="95% Confidence",
    )
    axes[0].scatter(
        x_observed, y_observed, c="red", s=80, zorder=5, label="Observed Points"
    )
    axes[0].set_ylabel("f(x)", fontsize=12)
    axes[0].set_title("GP Surrogate Model", fontsize=14)
    axes[0].legend(fontsize=10)
    axes[0].grid(True, alpha=0.3)

What is the top panel?

The top panel is the same GP surrogate from your first plot, repeated here so you can see how the acquisition function below relates to the model's predictions and uncertainty bands.

  • Directly below the top panel code, add the bottom panel with the UCB curve and the entry point for the script:
    # Bottom: Acquisition function
    axes[1].plot(x_grid, ucb, "g-", linewidth=2, label="UCB Acquisition Function")
    next_x_idx = np.argmax(ucb)
    next_x = x_grid[next_x_idx][0]
    axes[1].axvline(
        x=next_x, color="red", linestyle=":", linewidth=2,
        label=f"Next Sample (x={next_x:.2f})",
    )
    axes[1].set_xlabel("x", fontsize=12)
    axes[1].set_ylabel("UCB(x)", fontsize=12)
    axes[1].set_title("Acquisition Function: Upper Confidence Bound", fontsize=14)
    axes[1].legend(fontsize=10)
    axes[1].grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig("acquisition_function.png", dpi=150)
    plt.show()
    print("Saved: acquisition_function.png")


if __name__ == "__main__":
    main()

What does the bottom panel show?

  • The green curve is the UCB value at every point. It combines the GP mean (exploitation) with the uncertainty (exploration).
  • The red dotted vertical line marks where UCB peaks. This is where BO would sample next if given another evaluation.
  • Notice that the peak is often in a region where the GP is uncertain, not necessarily where the mean is highest. That is the exploration/exploitation tradeoff in action.
  • Save visualize_bo.py.

Before you run, think about this: where do you predict the red "Next Sample" line will appear? Will it be near an observed point, or in a gap where no observations exist?

  • Run the script to generate both plots:
python visualize_bo.py

Two matplotlib windows appear in sequence. After closing both, you should see:

Saved: gp_surrogate.png
Saved: acquisition_function.png

The "Next Sample" line likely points to a region with high uncertainty. That is BO choosing to explore where it knows the least.

Seeing a TypeError or shape mismatch?

  • Check that x_grid uses .reshape(-1, 1) and that fill_between uses x_grid.flatten(). The GP expects a 2D array for prediction but matplotlib needs a 1D array for the x-axis.
  • Make sure each code chunk is indented inside def main():. All the plotting code lives inside the function.

Help me fix a shape error in my visualize_bo.py script.

✔️ Awesome, I've got everything!

Great. Double check you have saved visualize_bo.py and that both gp_surrogate.png and acquisition_function.png appear in your project folder.

ⓧ I'd like to double check the full code

Your visualize_bo.py is built from the chunks above. Scroll up and confirm you have all of them in this order:

  • Imports and mystery_function definition.
  • def main(): with the BO run and observed point extraction.
  • GP fitting with GaussianProcessRegressor and prediction grid.
  • Plot 1: GP Surrogate Model saved as gp_surrogate.png.
  • Plot 2 top panel: GP model repeated.
  • Plot 2 bottom panel: UCB acquisition function saved as acquisition_function.png.
  • if __name__ == "__main__": main() at the very end (unindented).

You now have visual proof of how Bayesian optimization makes its decisions. The GP models the function from sparse observations, and the acquisition function uses that model to pick the most informative next point. Next up, you will apply this same technique to a real machine learning model and watch BO find better hyperparameters than random search.

Tune a Real ML Model with Bayesian Optimization

Your visualizations prove that Bayesian optimization works on a math function. The real payoff comes when you point it at something expensive to evaluate, like training a machine learning model.

In this step, you will apply BO to tune a Random Forest classifier on a real dataset. Each evaluation means training the model and running 5-fold cross-validation. Random search wastes those evaluations blindly. BO spends them strategically.

In this step, get ready to:
  • Define an objective function that trains and scores a Random Forest classifier.
  • Run random search as a baseline, then beat it with Bayesian optimization using the same evaluation budget.
  • Print the best hyperparameters and accuracy improvement BO achieved.
Create the tuning script and define the objective function

Hyperparameter tuning is expensive. Each "evaluation" trains a full model and cross-validates it. You need an objective function that takes hyperparameters as input and returns a score.

  • Create a new file called tune_model.py in your project folder and paste the following code:
"""Apply Bayesian optimization to tune a real ML model."""

import random
import numpy as np
from sklearn.datasets import load_wine
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
from bayes_opt import BayesianOptimization


def main():
    # Load dataset
    data = load_wine()
    X, y = data.data, data.target

    print(f"Dataset: Wine ({X.shape[0]} samples, {X.shape[1]} features, 3 classes)")
    print()

    # Define the objective function
    def rf_objective(n_estimators, max_depth, min_samples_split):
        clf = RandomForestClassifier(
            n_estimators=int(n_estimators),
            max_depth=int(max_depth),
            min_samples_split=int(min_samples_split),
            random_state=42,
        )
        scores = cross_val_score(clf, X, y, cv=5, scoring="accuracy")
        return scores.mean()

What does this code do?

The Wine dataset from scikit-learn has 178 samples across 3 wine classes, with 13 chemical features. It is small enough to train quickly but complex enough that hyperparameters matter.

The rf_objective function is what BO will call repeatedly. It takes three hyperparameters (n_estimators, max_depth, min_samples_split), trains a Random Forest, runs 5-fold cross-validation, and returns the mean accuracy.

Add the random search baseline

Before you let BO loose, you need a fair baseline. Random search will pick 20 random hyperparameter combinations from the same search space.

  • In tune_model.py, add the following code directly below the return scores.mean() line:
    # --- Random Search Baseline ---
    print("=" * 50)
    print("RANDOM SEARCH (20 evaluations)")
    print("=" * 50)

    random.seed(42)
    best_random_score = 0
    best_random_params = {}

    for i in range(20):
        n_est = random.randint(10, 300)
        max_d = random.randint(2, 20)
        min_s = random.randint(2, 20)
        score = rf_objective(n_est, max_d, min_s)
        if score > best_random_score:
            best_random_score = score
            best_random_params = {
                "n_estimators": n_est,
                "max_depth": max_d,
                "min_samples_split": min_s,
            }
            print(f"  New best at eval {i + 1}: accuracy = {score:.4f}")

    print(f"\nBest Random Search accuracy: {best_random_score:.4f}")
    print(f"Best params: {best_random_params}")


if __name__ == "__main__":
    main()

What does this code do?

The loop generates 20 random hyperparameter sets within the bounds (n_estimators 10-300, max_depth 2-20, min_samples_split 2-20). Each set trains a full model and cross-validates it. The best score is tracked as the baseline to beat.

The if __name__ block at the bottom makes the script runnable. You will add more code above it in the next chunk.

  • Save tune_model.py and run the random search baseline:
python tune_model.py

You should see the dataset info printed, followed by "New best" lines as random search stumbles onto better hyperparameter combos. Note the final accuracy it achieves.

Seeing an import error?

Make sure your virtual environment is still active. You should see (venv) at the start of your terminal prompt. If not, activate it again with venv\Scripts\activate (Windows) or source venv/bin/activate (macOS/Linux).

help me debug an import error when running tune_model.py

Add Bayesian optimization and compare the results

Random search gave you a baseline accuracy. Now you will point BO at the same search space with the same total budget of 20 evaluations and see if it can do better.

  • In tune_model.py, find the if __name__ == "__main__": line at the bottom of the file.
  • Add the following code above that line (between the random search and the if __name__ block):
    # --- Bayesian Optimization ---
    print("\n" + "=" * 50)
    print("BAYESIAN OPTIMIZATION (20 evaluations)")
    print("=" * 50)

    optimizer = BayesianOptimization(
        f=rf_objective,
        pbounds={
            "n_estimators": (10, 300),
            "max_depth": (2, 20),
            "min_samples_split": (2, 20),
        },
        random_state=42,
        verbose=2,
    )

    optimizer.maximize(init_points=5, n_iter=15)

    # --- Final Comparison ---
    print("\n" + "=" * 50)
    print("RESULTS COMPARISON")
    print("=" * 50)
    print(f"  Random Search best accuracy: {best_random_score:.4f}")
    print(f"  Bayesian Opt best accuracy:  {optimizer.max['target']:.4f}")
    print(f"\n  BO best parameters:")
    for param, value in optimizer.max["params"].items():
        print(f"    {param}: {int(value)}")
    improvement = optimizer.max["target"] - best_random_score
    print(f"\n  Accuracy improvement: {improvement:+.4f}")

What does this code do?

The optimizer uses 5 random initial points and 15 guided iterations for the same total budget of 20. It searches the same hyperparameter space as random search, but uses the GP surrogate to pick smarter combinations.

The final comparison prints both methods' best accuracy side by side, the winning hyperparameters, and the exact accuracy improvement BO achieved.

  • Save tune_model.py.

Before you run this, what do you think will happen? Will BO beat random search on accuracy with the same 20 evaluations? By how much?

  • Run the complete tuning script:
python tune_model.py

You will see the random search results first, then BO's iteration table showing it intelligently probing the hyperparameter space. The final RESULTS COMPARISON shows BO achieving higher cross-validation accuracy than random search.

Why does BO win here?

Each evaluation trains a Random Forest and runs 5-fold cross-validation. That is expensive. Random search spends those evaluations blindly. BO uses each result to update its surrogate model of the accuracy surface, then picks the next hyperparameter combination where it predicts the best return.

With only 20 evaluations over a 3-dimensional space, the difference between blind sampling and guided sampling becomes significant.

Script takes a long time to finish?

Each evaluation trains a Random Forest with up to 300 trees and cross-validates it 5 times. The full script should finish in under a minute on most machines. If it is significantly slower, check that no other heavy processes are running.

help me figure out why my tune_model.py is running slowly

✔️ Awesome, I've got everything!

Great. Double check you have saved tune_model.py and confirmed BO beats random search in the printed output.

ⓧ I'd like to double check the full code

Here is the complete tune_model.py file. Scroll up and confirm each section is present in this order:

  • Imports and def main(): with dataset loading and objective function definition.
  • Random Search Baseline section (the for i in range(20) loop).
  • Bayesian Optimization section (optimizer.maximize(init_points=5, n_iter=15)).
  • Final Comparison section printing both scores and the improvement.
  • The if __name__ == "__main__": block at the very end.

You have now seen Bayesian optimization outperform random search on a real ML task. With the same compute budget, BO found better hyperparameters by learning from each evaluation instead of guessing blindly. Next up, explore how the acquisition function's exploration parameter changes BO's strategy in the Secret Mission.

Secret mission

Explore Acquisition Function Strategies

The acquisition function controls whether Bayesian optimization explores unknown regions or exploits known good ones. In this mission, you'll run BO with two different kappa values and visualize how each strategy samples the search space differently.

Clean Up Your Resources

Clean Up Your Resources

Decide whether to keep your resources running, pause them to come back later, or delete them entirely. This project runs entirely locally, so there are no ongoing costs.

Resources you used:

  • A Python virtual environment (venv/ folder) on your local machine.
  • Several PNG visualization files (random_search_convergence.png, convergence_comparison.png, gp_surrogate.png, acquisition_function.png, acquisition_comparison.png).
  • A git repository with Python scripts (random_search.py, bayesian_opt.py, visualize_bo.py, tune_model.py).

Keep everything running

No action needed. Choose this if you plan to keep experimenting with Bayesian optimization or want to push the project to GitHub as a portfolio piece.

  • The project folder stays on your Desktop. There are no ongoing costs.
  • You can push to GitHub anytime. The repo is self-contained with requirements.txt and .gitignore already configured.

Pause - I'll come back to this later

Deactivate your virtual environment to free up your terminal session. Your files stay on disk and you can reactivate anytime.

  • Deactivate the virtual environment by running:
deactivate
  • When you come back, reactivate with:

macOS

source venv/bin/activate

Windows

venv\Scripts\activate

Delete - I don't want to use this again

Remove all project resources and start fresh. This deletes your scripts, visualizations, and virtual environment permanently.

  • Deactivate the virtual environment first:
deactivate
  • Move out of the project folder and delete it:

macOS

cd ~/Desktop
rm -rf bayesian-optimization

Windows

cd ~/Desktop
rmdir /s /q bayesian-optimization
  • Confirm the folder is gone by running:
ls
  • You should no longer see the bayesian-optimization folder in the output.

Nice Work!

Nice Work!

Nice work. You just built a complete Bayesian optimization project from scratch, starting with the pain of random search and ending with a production-ready hyperparameter tuning workflow.

You've learned how to:

  • Use Bayesian optimization to find better optima than random search with the same evaluation budget, proving that intelligent sampling beats brute force.
  • Visualize the Gaussian Process surrogate model, confidence bands, and Upper Confidence Bound acquisition function to reveal exactly how BO decides where to sample next.
  • Apply BO to tune a real Random Forest classifier on the Wine dataset, achieving higher cross-validation accuracy than random search with the same compute budget.
  • Secret Mission: Compared exploration-heavy vs exploitation-heavy acquisition strategies by tuning the kappa parameter and visualized the tradeoff.

Ready to quiz yourself?