Рефакторинг оставил рекурсию бесконечного цикла

Question

При попытке реализовать эту формулу

простым способом, используя этот код

def P(n, p, limit):
    if n == limit:
        return 1

    if n == -limit:
        return 0

    return p*P(n+1, p, limit) + (1-p)*P(n-1, p, limit)

вы обязательно получите эту ошибку.

RecursionError: превышена максимальная глубина рекурсии

Итак, как правильно реализовать такой класс проблем с помощью кода? Есть ли метод, используемый для рефакторинга этого и получения решения? Или единственное, что нужно решить сначала математически, а потом только потом - в вычислениях?

Answer 1

Чтобы расширить предложения @Prune, это повторение можно записать как два других повторения, перенеся любой из терминов RHS в LHS:

Они нециклические, потому что каждый термин зависит только от терминов до или после него, а не от обоих.

---

Однако они не могут быть решены независимо, потому что это требует знания P _1-l и P _l-1 , где l = limit.

Выполнение некоторых замен для удобства:

Обратите внимание, что, поскольку отношения содержат только линейные члены (нет, например, P _n × P _n-1 и т. Д.), Все члены в их расширениях также будет линейным по переменным x, y:

Эти недавно введенные коэффициенты также подчиняются аналогичным рекуррентным соотношениям:

Их граничные условия:

Теперь, когда граничные условия известны, коэффициенты могут быть вычислены итеративно. Осталось только уравнять выражения для P _1-l и P _l-1 из обоих повторений:

Решение этой пары уравнений дает P _1-l и P _l-1 , из которых можно вычислить любой P _n (используя коэффициенты вычислено ранее):

---

Обновление : пример реализации Python

# solve a recurrence of the form
# Tn = c1 * Tn-1 + c2 * Tn-2, n >= 0
def solve_recurrence(c1, c2, t1, t2, n):
    if n <= 0: return t2
    if n == 1: return t1
    a, b = t2, t1
    for i in range(1, n):
        a, b = b, c1 * b + c2 * a
    return b

# solve the loop recurrence
def solve_loop(limit, p, n):
    assert(limit > 0 and abs(n) <= limit)
    assert(p > 0 and p < 1)

    # corner cases
    if n == -limit: return 0
    if n ==  limit: return 1

    # constants
    a, c = 1 / p, 1 / (1 - p)
    b, d = 1 - a, 1 - c

    # coefficients at 2l - 1
    e = 2 * limit - 1
    l = solve_recurrence(a, b, 1, 0, e)
    m = solve_recurrence(a, b, 0, 0, e)
    s = solve_recurrence(c, d, 1, 0, e)
    t = solve_recurrence(c, d, 0, 1, e)
    # note that m can just be 0 but was left in for generality

    # solving the simultaneous equations
    # to get P at 1 - l and l - 1
    x = (s * m + t) / (1 - s * l)
    y = (l * t + m) / (1 - l * s)

    # iterate from the closest side to improve precision
    if n < 0:
        return solve_recurrence(a, b, x, 0, limit + n)
    else:
        return solve_recurrence(c, d, y, 1, limit - n)

Экспериментальные результаты для limit = 10, p = 0.1:

n       | P(n)            abs. error
========================================
-10     | 0.0000000E+00   --
-9      | 6.5802107E-19   0.0000000E+00
-8      | 6.5802107E-18   1.7995889E-30
-7      | 5.9879917E-17   1.7995889E-29
-6      | 5.3957728E-16   5.0003921E-28
-5      | 4.8568535E-15   8.1994202E-27
-4      | 4.3712339E-14   8.9993252E-27
-3      | 3.9341171E-13   5.1002066E-25
-2      | 3.5407061E-12   5.9938283E-25
-1      | 3.1866355E-11   5.9339980E-18
 0      | 2.8679726E-10   5.3406100E-17
 1      | 2.5811749E-09   3.3000371E-21
 2      | 2.3230573E-08   2.6999175E-20
 3      | 2.0907516E-07   2.2999591E-19
 4      | 1.8816764E-06   7.9981086E-19
 5      | 1.6935088E-05   1.9989978E-18
 6      | 1.5241579E-04   3.5000757E-16
 7      | 1.3717421E-03   1.5998487E-15
 8      | 1.2345679E-02   3.2000444E-14
 9      | 1.1111111E-01   6.9999562E-14
 10     | 1.0000000E+00   --

---

Примечания:

• Из кода и уравнений ясно, что пределы повторяемости / граничные условия могут быть произвольными.
• Числовая точность может быть улучшена с помощью методов суммирования, таких как Kahan-Neumaier.
• Существует также альтернативная явная формула для итерационного метода, используемого в solve_recurrence, которая может давать более точные результаты, если используется вместе с библиотечными функциями.

Answer 2

У меня есть другое решение, использующее символические вычисления.

class X(object):
    """Representing unknown value which will be resolved in future"""
    def __init__(self, value=1.0):
        self.value = value

    def __add__(self, other):
        cls = self.__class__
        return cls(value=self.value+other.value)

    def __sub__(self, other):
        cls = self.__class__
        return cls(value=self.value-other.value)

    def __mul__(self, other):
        cls = self.__class__
        if isinstance(other, cls):
            return cls(value=self.value*other.value)
        elif isinstance(other, numbers.Real):
            return cls(value=self.value*other)
        else:
            raise TypeError(f'Cannot multiply with type {type(other)}')

    def __rmul__(self, other):
        return self.__mul__(other)

Вместо немедленного получения числовых результатов мы полагаемся на тот факт, что полученное значение будет каким-то образом выражено как

и каждый промежуточный член будет представлен как

Именно поэтому мы определяем операции добавления / вычитания / умножения для нашего объекта.

Таким образом, имея представление такого символа в коде, мы можем решить это уравнение с помощью подхода, аналогичного принятому ответу (мы исходим из сторон -limit и limit и выражаем другие термины как еще неизвестный продукт). При достижении базовой точки (0 состояние) у нас будет отношение нашей вероятности и числовое значение для этого отношения, что позволяет нам легко вычислить результат

_memory = []
def fill_memory(n):
    global _memory
    _memory = [(0, X(1))]*(2*n+1)
    _memory[-n] = (0, X(0))
    _memory[n] = (1, X(0))


def resolve(n, p):
    q = 1-p
    for i in reversed(range(n)):
        # resolve from right to center
        prob_next = _memory[i+1]
        prob_prev = _memory[i-i]
        _memory[i] = (
            q*prob_prev[0] + p*prob_next[0],
            q*prob_prev[1] + p*prob_next[1],
        )

        # resolve from left to center
        prob_prev = _memory[-i-1]
        prob_next = _memory[-i+1]
        _memory[-i] = (
            q*prob_prev[0] + p*prob_next[0],
            q*prob_prev[1] + p*prob_next[1],
        )

    # calculate base probability
    solution = _memory[0]
    coef = (X(value=1) - solution[1]).value
    p0 = solution[0] / coef
    return p0

Я получил правильные и очень точные результаты с этим подходом (я использую p=0.6 и N=4 (предел))

def main():
    N = 4
    p = 0.6
    fill_memory(N)
    p0 = resolve(N, p)
    print(f'Probability for N={N} steps is {p0*100:.4}%')
    for i in range(2*N+1):
        print(i-N, '->', _memory[i-N])

4 - 83.505%
3 - 77.143%
2 - 69.231%
1 - 60.0%

Как вы можете видеть, это тот же результат, что и при математическом разрешении.

Answer 3

Рекурсия зависит от разбиения задачи на части, каждая из которых ближе к базовому варианту. Формула, которую вы даете, бесконечна, потому что она зависит от ряда терминов, которые одинаково не определены. Простая математическая запись не является гарантией вычислимости.

В этом случае вы должны переписать отношение рекурсии в некруговых терминах. Вы все еще можете использовать рекурсию, но должен быть определенный конец. В данной формуле P [n] зависит от P [n + 1] и P [n-1], каждый из которых напрямую зависит от P [n]. Это алгоритмически круглая.

Вам нужно извлечь алгоритм из базовых случаев, начиная с конца и работая в своем направлении. Одним из неудачных результатов этого рекурсивного отношения является то, что вам нужно решить весь диапазон [-limit, limit] для чтобы получить значение для желаемой позиции. Значения базового случая колеблются от концов до того места, где у вас есть точка зрения.

Нет никакого инструмента, чтобы символически разгадать это для вас. Прежде всего, обратите внимание, что ваши базовые случаи не отображаются в уравнении. Даже если вы предоставили их в качестве дополнительных уравнений, вывод ряда экспонент из этого отношения выходит за рамки возможностей анализа сегодняшних решателей.