def METHOD_NAME(
    num_iters: int, batch_size: int, n_data: int, verbose: bool
) -> None:
    # Create dataset:
    key = jr.PRNGKey(123)
    x = jnp.sort(
        jr.uniform(key=key, minval=-2.0, maxval=2.0, shape=(n_data, 1)), axis=0
    )
    y = jnp.sin(x) + jr.normal(key=key, shape=x.shape) * 0.1
    D = Dataset(X=x, y=y)
    # Define GP model:
    prior = Prior(kernel=RBF(), mean_function=Constant())
    likelihood = Gaussian(num_datapoints=n_data)
    posterior = prior * likelihood
    # Define variational family:
    z = jnp.linspace(-2.0, 2.0, 10).reshape(-1, 1)
    q = VariationalGaussian(posterior=posterior, inducing_inputs=z)
    # Define loss function:
    elbo = ELBO(negative=True)
    # Train!
    trained_model, history = fit(
        model=q,
        objective=elbo,
        train_data=D,
        optim=ox.adam(0.1),
        num_iters=num_iters,
        batch_size=batch_size,
        verbose=verbose,
        key=jr.PRNGKey(123),
    )
    # Ensure the trained model is a Gaussian process posterior
    assert isinstance(trained_model, VariationalGaussian)
    # Ensure we return a history of the correct length
    assert len(history) == num_iters
    # Ensure we reduce the loss
    assert elbo(trained_model, D) < elbo(q, D)