def arrival_times(n=1):
    fun1 = lambda x: (2 * x) / 5 + 15
    fun2 = lambda x: (56 - x) / 3
    fun3 = lambda x: (88 + 3 * x) / 7
    fun4 = lambda x: (83 - 3 * x) / 2
    fun5 = lambda x: (129 - x) / 7

    integrals = (integrate.quad(fun1, 0, 5), integrate.quad(fun2, 5, 8), integrate.quad(fun3, 8, 15),
                 integrate.quad(fun4, 15, 17), integrate.quad(fun5, 17, 24))

    integrals_res = [round(x[0], 4) for x in integrals]
    print(integrals_res)

    xk = (0, 1, 2, 3, 4)
    temp = []

    for x in integrals_res:
        temp.append(np.exp(-x) * x ** n / np.math.factorial(n))

    pk = tuple(temp)

    print(pk)

    custom = stats.rv_discrete(name='custom', values=(xk, pk))
    return custom.rvs(size=n)