How to use the pymc3.find_MAP function in pymc3

approx = pm.fit(
            random_seed=random_seed,
            n=n_init,
            method="advi",
            model=model,
            callbacks=cb,
            progressbar=progressbar,
            obj_optimizer=pm.adagrad_window,
        )  # type: pm.MeanField
        start = approx.sample(draws=chains)
        start = list(start)
        stds = approx.bij.rmap(approx.std.eval())
        cov = model.dict_to_array(stds) ** 2
        potential = quadpotential.QuadPotentialDiag(cov)
    elif init == "advi_map":
        start = pm.find_MAP(include_transformed=True)
        approx = pm.MeanField(model=model, start=start)
        pm.fit(
            random_seed=random_seed,
            n=n_init,
            method=pm.KLqp(approx),
            callbacks=cb,
            progressbar=progressbar,
            obj_optimizer=pm.adagrad_window,
        )
        start = approx.sample(draws=chains)
        start = list(start)
        stds = approx.bij.rmap(approx.std.eval())
        cov = model.dict_to_array(stds) ** 2
        potential = quadpotential.QuadPotentialDiag(cov)
    elif init == "map":
        start = pm.find_MAP(include_transformed=True)

"z",
                mu=0,
                sd=(1.0 - smoothing_param) * sd,
                shape=y.shape,
            )

            nu = pymc3.Gamma("nu", alpha=2, beta=1)
            obs = pymc3.StudentT(
                "obs",
                mu=z,
                sd=sd * smoothing_param,
                nu=nu,
                observed=y,
            )

            res = pymc3.find_MAP(vars=[z], method="L-BFGS-B")
            return res['z']

X1 = np.random.randn(size)
X2 = np.random.randn(size) * 0.2
Y = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(size)*sigma

basic_model = Model()
with basic_model:
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=2)
    sigma = HalfNormal('sigma', sd=1)
    mu = alpha + beta[0]*X1 + beta[1]*X2
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
    
print(basic_model)
    
map_estimate = find_MAP(model=basic_model)
print(map_estimate)
#{'alpha': array(0.9065985664354854), 'beta': array([ 0.948486  ,  2.60705513]), 'sigma_log': array(-0.03278146854842092)}

map_estimate2 = find_MAP(model=basic_model, fmin=opt.fmin_powell)
print(map_estimate2)

map_estimate3 = find_MAP(model=basic_model, fmin=opt.fmin_l_bfgs_b)
print(map_estimate3)

with basic_model:
    #step = Slice(vars=[sigma]) 
    #trace = sample(2000, start=map_estimate, step=step) 
    trace = sample(2000, start=map_estimate) 

beta_samples = trace['beta'] # shape 2000, 2
traceplot(trace);

returns[:5]

plt.plot(returns)
plt.ylabel('daily returns in %');


with pymc3.Model() as sp500_model:
    nu = pymc3.Exponential('nu', 1./10, testval=5.)
    sigma = pymc3.Exponential('sigma', 1./.02, testval=.1)
    s = ts.GaussianRandomWalk('s', sigma**-2, shape=len(returns))
    volatility_process =  pymc3.Deterministic('volatility_process', pymc3.exp(-2*s))
    r = pymc3.StudentT('r', nu, lam=1/volatility_process, observed=returns)
    
with sp500_model:
    print('optimizing...')
    start = pymc3.find_MAP(vars=[s], fmin=scipy.optimize.fmin_l_bfgs_b)
    
    print('sampling... (slow!)')
    step = pymc3.NUTS(scaling=start)
    trace = pymc3.sample(100, step, progressbar=False)

    # Start next run at the last sampled position.
    step = pymc3.NUTS(scaling=trace[-1], gamma=.25)
    trace = pymc3.sample(1000, step, start=trace[-1], progressbar=False)
    

pymc3.traceplot(trace, [nu, sigma]);

fig, ax = plt.subplots(figsize=(15, 8))
#returns.plot(ax=ax)
ax.plot(returns)
ax.plot(1/np.exp(trace['s',::30].T), 'r', alpha=.03);

def run(n=2000):
    if n == "short":
        n = 50
    import matplotlib.pyplot as plt

    with model:
        start = find_MAP(fmin=opt.fmin_powell)
        trace = sample(n, Slice(), start=start)

    plt.plot(x, y, 'x')
    glm.plot_posterior_predictive(trace)
    # plt.show()

light_curve = op.get_light_curve(orbit=orbit, r=r, t=t, y=y) + mean
    
    # Here we track the value of the model light curve for plotting
    # purposes
    pm.Deterministic("light_curve", light_curve)
    
    # In this line, we simulate the dataset that we will fit
    flux = xo.eval_in_model(light_curve)
    flux += ferr * np.random.randn(len(flux))
    
    # The likelihood function assuming known Gaussian uncertainty
    pm.Normal("obs", mu=light_curve, sd=ferr, observed=flux)
    
    # Fit for the maximum a posteriori parameters given the simuated
    # dataset
    map_soln = pm.find_MAP(start=model.test_point)


np.random.seed(42)
sampler = xo.PyMC3Sampler(window=100, finish=200)
with model:
    burnin = sampler.tune(tune=2500, start=map_soln, step_kwargs=dict(target_accept=0.9))
    trace = sampler.sample(draws=3000)

with pm.Model() as model:
    S = np.eye(3)
    nu = 5
    mu = pm.Normal('mu', mu=0, sd=1, shape=3)

    # Use the transformed Wishart distribution
    # Under the hood this will do a Cholesky decomposition
    # of S and add two RVs to the sampler: c and z
    prec = pm.WishartBartlett('prec', S, nu)

    # To be able to compare it to truth, convert precision to covariance
    cov = pm.Deterministic('cov', tt.nlinalg.matrix_inverse(prec))

    lp = pm.MvNormal('likelihood', mu=mu, tau=prec, observed=data)

    start = pm.find_MAP()
    step = pm.NUTS(scaling=start)


def run(n=3000):
    if n == "short":
        n = 50
    with model:
        trace = pm.sample(n, step, start)

    pm.traceplot(trace)

if __name__ == '__main__':
    run()

cov = eta**2 * pymc3.gp.cov.Matern52(1, ell)
            gp = pymc3.gp.Latent(cov_func=cov)

            f = gp.prior("f", X=x.reshape(-1, 1))

            sigma = pymc3.HalfCauchy("sigma", beta=2)
            # sigma = pymc3.Normal("sigma")
            # sigma = 0.1
            nu = pymc3.Gamma("nu", alpha=2, beta=1)
            # sigma = 0.01
            # nu = 0.01
            # y_ = pymc3.StudentT("y", mu=f, lam=1.0/sigma, nu=nu, observed=y)
            y_ = pymc3.StudentT("y", mu=f, sd=sigma, nu=nu, observed=y)

            # res = pymc3.find_MAP(model=model, method="L-BFGS-B")
            res = pymc3.find_MAP(vars=[f], method="L-BFGS-B")
            return res["f"]

pm.summary(advi_trace)

basic_model = pm.Model()

with basic_model:
    # Priors
    mu = pm.Normal('mu', mu=0.0, sd=5.0)
    inv_softplus_sigma = pm.Normal('inv_softplus_sigma', mu=0.0, sd=1.0)

    # Model
    y = pm.Normal(
        'y', mu=mu, sd=tt.nnet.softplus(inv_softplus_sigma), observed=data
    )

    print('Sampling based approach')
    params = pm.find_MAP(model=basic_model)
    print(params)