How to model the parameters of AssymetricLaplace and Pareto distribution

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()