<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def hsepd_pdf(sigma1, sigma2, xi, beta, sim=None, obs=None, node=None, skip_nan=False):
 """Calculate the probability densities based on the heteroskedastic skewed exponential power distribution. For convenience, the required parameters of the probability density function as well as the simulated and observed values are stored in a dictonary: The following test function allows the variation of one parameter and prints some and plots all of probability density values corresponding to different simulated values: When varying parameter `beta`, the resulting probabilities correspond to the Laplace distribution (1.0), normal distribution (0.0), and the uniform distribution (-1.0), respectively. Note that we use -0.99 instead of -1.0 for approximating the uniform distribution to prevent from running into numerical problems, which are not solved yet: 10.0, 0.002032, 0.000886, 0.0 15.0, 0.008359, 0.010798, 0.0 20.0, 0.034382, 0.048394, 0.057739 25.0, 0.141421, 0.079788, 0.057739 30.0, 0.034382, 0.048394, 0.057739 35.0, 0.008359, 0.010798, 0.0 40.0, 0.002032, 0.000886, 0.0 .. testsetup:: When varying parameter `xi`, the resulting density is negatively skewed (0.2), symmetric (1.0), and positively skewed (5.0), respectively: 10.0, 0.0, 0.000886, 0.003175 15.0, 0.0, 0.010798, 0.012957 20.0, 0.092845, 0.048394, 0.036341 25.0, 0.070063, 0.079788, 0.070063 30.0, 0.036341, 0.048394, 0.092845 35.0, 0.012957, 0.010798, 0.0 40.0, 0.003175, 0.000886, 0.0 .. testsetup:: In the above examples, the actual `sigma` (5.0) is calculated by multiplying `sigma1` (0.2) with the mean simulated value (25.0), internally. This can be done for modelling homoscedastic errors. Instead, `sigma2` is multiplied with the individual simulated values to account for heteroscedastic errors. With increasing values of `sigma2`, the resulting densities are modified as follows: 10.0, 0.000886, 0.002921, 0.005737 15.0, 0.010798, 0.018795, 0.022831 20.0, 0.048394, 0.044159, 0.037988 25.0, 0.079788, 0.053192, 0.039894 30.0, 0.048394, 0.04102, 0.032708 35.0, 0.010798, 0.023493, 0.023493 40.0, 0.000886, 0.011053, 0.015771 .. testsetup:: """

    sim, obs = prepare_arrays(sim, obs, node, skip_nan)
    sigmas = _pars_h(sigma1, sigma2, sim)
    mu_xi, sigma_xi, w_beta, c_beta = _pars_sepd(xi, beta)
    x, mu = obs, sim
    a = (x-mu)/sigmas
    a_xi = numpy.empty(a.shape)
    idxs = mu_xi+sigma_xi*a < 0.
    a_xi[idxs] = numpy.absolute(xi*(mu_xi+sigma_xi*a[idxs]))
    a_xi[~idxs] = numpy.absolute(1./xi*(mu_xi+sigma_xi*a[~idxs]))
    ps = (2.*sigma_xi/(xi+1./xi)*w_beta *
          numpy.exp(-c_beta*a_xi**(2./(1.+beta))))/sigmas
    return ps