It sounds like you are trying to do optimization rather than Bayesian inference. If so, then I would suggest using scipy’s built-in optimization routines.
Alternatively, it seems like the kind of thing one could do would be the following:
obs = np.array([1.726567e+06, 1.589836e+06, 1.643981e+06, 1.584314e+06])
with pm.Model() as model1:
kappa = pm.HalfNormal("kappa", 10)
mu = pm.Normal("mu", mu=np.mean(obs), sigma=np.std(obs))
b = pm.HalfNormal("b", 3 * np.std(obs))
# Based on distr of raw data
likelihood = pm.AsymetricLaplace("likelihood", kappa = kappa, mu = mu, b = b, observed=obs)
idata = pm.sample()