Minimal example below 


import numpy as np
import pandas as pd
import pymc as pm


def generate_data(num_samples: int) -> pd.DataFrame:
    rng = np.random.default_rng(seed=42)
    beta = 1.0
    sigma = 10.0
    x = rng.normal(loc=0.0, scale=1.0, size=num_samples)
    y = beta * x   sigma * rng.normal(size=num_samples)
    return pd.DataFrame({"x": x, "y": y})


def make_model(frame: pd.DataFrame) -> pm.Model:
    with pm.Model() as model:
        # Data
        x = pm.Data("x", frame["x"])
        y = pm.Data("y", frame["y"])

        # Prior
        beta = pm.Normal("beta", sigma=10.0)
        sigma = pm.HalfNormal("sigma", sigma=20.0)

        # Linear model
        mu = beta * x

        # Likelihood
        pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)
    return model

if __name__ == "__main__":
    frame = generate_data(10000)
    model = make_model(frame)
    with model:
        advi = pm.ADVI(
            start={"beta": 1.0, "sigma": 10.0},
            start_sigma={"beta": 0.5, "sigma_log__": 0.5}
        )
       tracker = Tracker(
            mean=advi.approx.mean.eval,
            std=advi.approx.std.eval
        )
        approx = pm.fit(
            n=1_000_000,
            method=advi,
            callbacks=[
                CheckParametersConvergence(diff="relative", tolerance=1e-3),
                CheckParametersConvergence(diff="absolute", tolerance=1e-3),
                tracker
            ],
        )

    def create_shared_params(self, start=None, start_sigma=None):
        # NOTE: `Group._prepare_start` uses `self.model.free_RVs` to identify free variables and
        # `DictToArrayBijection` to turn them into a flat array, while `Approximation.rslice` assumes that the free
        # variables are given by `self.group` and that the mapping between original variables and flat array is given
        # by `self.ordering`. In the cases I looked into these turn out to be the same, but there may be edge cases or
        # future code changes that break this assumption.
        start = self._prepare_start(start)
        rho1 = np.zeros((self.ddim,))

        if start_sigma is not None:
            for name, slice_, *_ in self.ordering.values():
                sigma = start_sigma.get(name)
                if sigma is not None:
                    rho1[slice_] = np.log(np.expm1(np.abs(sigma)))
        rho = rho1

        return {
            "mu": pytensor.shared(pm.floatX(start), "mu"),
            "rho": pytensor.shared(pm.floatX(rho), "rho"),
        }

 _select = dict(advi=ADVI, fullrank_advi=FullRankADVI, svgd=SVGD, asvgd=ASVGD)
    if isinstance(method, str):
        method = method.lower()
        if method in _select:
            inference = _select[method](model=model, **inf_kwargs)
        else:
            raise KeyError(f"method should be one of {set(_select.keys())} or Inference instance")
    elif isinstance(method, Inference):
        inference = method
...

import numpy as np
import pandas as pd
import pymc as pm


def generate_data(num_samples: int) -> pd.DataFrame:
    rng = np.random.default_rng(seed=42)
    beta = 1.0
    sigma = 10.0
    x = rng.normal(loc=0.0, scale=1.0, size=num_samples)
    y = beta * x   sigma * rng.normal(size=num_samples)
    return pd.DataFrame({"x": x, "y": y})


def make_model(frame: pd.DataFrame) -> pm.Model:
    with pm.Model() as model:
        # Data
        x = pm.Data("x", frame["x"])
        y = pm.Data("y", frame["y"])

        # Prior
        beta = pm.Normal("beta", sigma=10.0)
        sigma = pm.HalfNormal("sigma", sigma=20.0)

        # Linear model
        mu = beta * x

        # Likelihood
        pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)
    return model


if __name__ == "__main__":
    num_samples = 10_000
    frame = generate_data(num_samples=num_samples)
    model = make_model(frame)
    with model:
        advi = pm.ADVI(
        )
        tracker = Tracker(
            mean=advi.approx.mean.eval,
            std=advi.approx.std.eval
        )
        t0 = time.time()
        approx = pm.fit(
            n=1_000_000,
            method=advi,
            callbacks=[
                CheckParametersConvergence(diff="relative", tolerance=1e-3),
                CheckParametersConvergence(diff="absolute", tolerance=1e-3),
                tracker
            ],
        )
        t = time.time() - t0
        print(f"Time for fit is {t:.3f}s.")

    frame_new = pd.concat([frame, generate_data(100)], axis=0)
    model = make_model(frame_new)
    with model:
        advi = pm.ADVI(
            start={"beta_foo": tracker["mean"][-1][0], "sigma": np.exp(tracker["mean"][-1][1])},
            start_sigma={"beta": tracker["std"][-1][0], "sigma_log__": tracker["std"][-1][1]}
        )
        tracker = Tracker(
            mean=advi.approx.mean.eval,
            std=advi.approx.std.eval
        )
        t0 = time.time()
        approx = pm.fit(
            n=1_000_000,
            method=advi,
            callbacks=[
                CheckParametersConvergence(diff="relative", tolerance=1e-3),
                CheckParametersConvergence(diff="absolute", tolerance=1e-3),
                tracker
            ],
        )
        t = time.time() - t0
        print(f"Time for fit is {t:.3f}s.")