Всё, что нужно знать о своих планах, случайностях и стохастическом программировании

Все мы прекрасно знаем, что очень часто наши планы идут не по плану именно из-за случайностей. В такие моменты очень трудно обойтись без жаргонизмов, нецензурной брани и отборного трехэтажного. Но все же есть способ сделать наши планы более устойчивыми и состоятельными — это стохастическое программирование (далее SP — stochastic programming).

Минимальнейший теоретический минимум

Как всегда, концептуально все достаточно просто. Допустим, мы хотим максимизировать вот такую функцию:

Что может быть банальнее параболы? Согласен, слишком просто. Но давайте представим, что по какой-то причине коэффициент при вдруг начал испытывать случайные отклонения: . Вообще, немного интересно взглянуть, что будет с параболой, если другие ее коэффициенты тоже будут случайными. Давайте нарисуем:

Код для отрисовки

import numpy as np
from scipy.stats import bernoulli, uniform, norm, expon
from scipy.integrate import quad
from scipy.optimize import linprog
from scipy.optimize import differential_evolution

import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")

fig, ax = plt.subplots(1, 3, figsize=(12, 4), dpi=150)

x = np.linspace(0.5, 5)
y_determ = -2 * x**2 + 12 * x - 14

iterations = 50

ax[0].plot(x, y_determ, 'C3', zorder=10)
for i in range(iterations):
    a = norm.rvs(loc=-2, scale = 0.1)
    y = a * x**2 + 12 * x - 14
    ax[0].plot(x, y, 'b', alpha=0.2)

ax[1].plot(x, y_determ, 'C3', zorder=10)
for i in range(iterations):
    b = uniform.rvs(loc=11, scale = 2)
    y = -2*x**2 + b*x - 14
    ax[1].plot(x, y, 'b', alpha=0.2)

ax[2].plot(x, y_determ, 'C3', zorder=10)
for i in range(iterations):
    c = expon.rvs(loc=-14, scale = 1)
    y = -2*x**2 + 12*x + c
    ax[2].plot(x, y, 'b', alpha=0.2)

ax[0].set_title(r'$a \sim N(-2, 0.1^{2})$' + '\n' + r'$b=12; c=-14$')
ax[1].set_title(r'$b \sim U(11, 13)$' + '\n' + r'$a=-2; c=-14$')
ax[2].set_title(r'$c \sim Exp(-14)$' + '\n' + r'$a=-2; b=12$')

ax[0].set_ylabel('Y')
ax[0].set_xlabel('x')
ax[1].set_xlabel('x')
ax[2].set_xlabel('x')

plt.show()

Красная линия — это детерминированная вселенная, очень хорошее место: с радугами и единорогами, там абсолютно все зависит только от наших желаний, решений и действий.

Реальность — это множество сценариев или даже множество возможных вселенных (синие линии), каждая из которых может реализоваться по воле случая, чаще всего с полным игнором наших желаний и надежд.

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

В этом месте может показаться, что теперь все, что нужно для счастья, это взять все параметры и переменные в наших моделях и заменить их на их pdf-ы! Но подождите, не все так просто. Давайте представим, что теперь все коэффициенты нашей параболы являются случайными. Если изобразить множество максимумов, то получим такое:

Код для отрисовки

x_vertex, y_vertex = [], []

N = 5000
a = norm.rvs(loc=-2, scale = 0.1, size=N)
b = uniform.rvs(loc=11, scale = 2, size=N)
c = expon.rvs(loc=-14, scale = 1, size=N)
x_vertex.append(- b / (2 * a))
y_vertex.append((4*a*c - b**2) / (4*a))

f, ax = plt.subplots(1, 3, figsize=(12, 4), dpi=150)

ax[0].scatter(x_vertex, y_vertex, s=1, alpha=0.2)
sns.kdeplot(x_vertex, ax=ax[1])
sns.kdeplot(y_vertex, ax=ax[2])

ax[0].set_title(r'KDE $Z^{*}$')
ax[0].set_xlabel(r'$Z^{*}_{x}$')
ax[0].set_ylabel(r'$Z^{*}_{y}$')

ax[1].set_title(r'KDE $Z^{*}_{x}$')
ax[2].set_title(r'KDE $Z^{*}_{y}$')
ax[1].set_xlabel('x')
ax[2].set_xlabel('y')

plt.show()

Что это за распределение? Я не знаю. Могу лишь сказать, что арифметика (и алгебра) для случайных величин работает не так, как для обычных чисел. Тут .

Даже в простейшем случае, например, когда и результатом арифметических операций будет вот такое:

Код для отрисовки

# Параметры равномерных распределений
a, b = 1e-10, 1  # X ~ U(0, 1)
c, d = 1e-10, 1  # Y ~ U(0, 1)

# Функции плотности
def f_X(x):
    return 1 / (b - a) if a <= x <= b else 0

def f_Y(y):
    return 1 / (d - c) if c <= y <= d else 0

# Численная свёртка для суммы Z = X + Y
def convolution_sum(z):
    # Пределы интегрирования: max(a, z - d) <= x <= min(b, z - c)
    lower = max(a, z - d)
    upper = min(b, z - c)
    if lower >= upper:
        return 0
    integral, _ = quad(lambda x: f_X(x) * f_Y(z - x), lower, upper)
    return integral

# Численная свёртка для произведения Z = X * Y
def convolution_prod(z):
    # Пределы интегрирования: x in [a, b] and [z/d, z/c] (для z > 0)
    if z == 0:
        return 0  # Особый случай
    lower = max(a, z / d)
    upper = min(b, z / c)
    if lower >= upper:
        return 0
    integral, _ = quad(lambda x: f_X(x) * f_Y(z / x) / abs(x), lower, upper)
    return integral

# Численная свёртка для частного Z = X / Y
def convolution_div(z):
    # Пределы интегрирования: y in [c, d] and [a/z, b/z] (для z > 0)
    if z == 0:
        return 0  # Особый случай
    lower = max(c, a / z)
    upper = min(d, b / z)
    if lower >= upper:
        return 0
    integral, _ = quad(lambda y: f_X(y * z) * f_Y(y) * abs(y), lower, upper)
    return integral

# Векторизация для массивов z
conv_sum_vec = np.vectorize(convolution_sum)
conv_prod_vec = np.vectorize(convolution_prod)
conv_div_vec = np.vectorize(convolution_div)

# Диапазон z для суммы (X + Y)
z_sum = np.linspace(1e-4, 2, 1000)
pdf_sum = conv_sum_vec(z_sum)

# Диапазон z для произведения (X * Y)
z_prod = np.linspace(5e-3, 1, 1000)
pdf_prod = conv_prod_vec(z_prod)

# Диапазон z для частного (X / Y)
z_div = np.linspace(1e-3, 5, 1000)
pdf_div = conv_div_vec(z_div)

# Монте-Карло для проверки
X = np.random.uniform(a, b, 100000)
Y = np.random.uniform(c, d, 100000)
Z_mc_sum = X + Y
Z_mc_prod = X * Y
Z_mc_div = X / Y
Z_mc_div = Z_mc_div[Z_mc_div < 100]

# Визуализация
f, ax = plt.subplots(1, 3, figsize=(12, 3), dpi=150)

ax[0].plot(z_sum, pdf_sum, c='C3', lw=2, label='PDF')
ax[0].hist(Z_mc_sum, bins=50, density=True, alpha=0.5, label='MC')
ax[0].set_title(r'$X + Y$')
ax[0].legend()

ax[1].plot(z_prod, pdf_prod, c='C3', lw=2, label='PDF')
ax[1].hist(Z_mc_prod, bins=50, density=True, alpha=0.5, label='MC')
ax[1].set_title(r'$X * Y$')
ax[1].legend()

ax[2].plot(z_div, pdf_div, c='C3', lw=2, label='PDF')
ax[2].hist(Z_mc_div, bins=1000, density=True, alpha=0.5, label='MC')
ax[2].set_xlim(0, 5)
ax[2].set_title(r'$X \; / \; Y$')
ax[2].legend()

plt.show()

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

и супермультипликативность, когда

Вот только эти эффекты могут быть как положительными, так и отрицательными. По опыту могу сказать, что если модель сложная, и решение одной задачи становится исходными данными для другой, то лучше заморочиться и проверить. Мое знакомство с SP началось именно так. Иногда к небольшим улучшениям приводит простая замена средних значений используемых величин на их распределения.

Даже не знаю, что еще добавить? Случайность — это неопределенность, а вероятность — это мера неопределенности, с которой мы умеем хорошо работать. Тем не менее совладать с реальностью все равно бывает непросто.

Так какая она, РЕАЛЬНОСТЬ?

Представим, что нам необходимо подготовить часть фундамента для небоскреба. Основной компонент бетона (60-70% объема) — это щебень (небольшие камешки от 5 до 40 мм с шероховатой поверхностью), который иногда получают путем перемалывания камней из гравия. В тоже время гравий содержит не только камни, которые можно перемолоть для щебня, но еще и песок, а он тоже является основным компонентом бетона.

Допустим, у нас всего два месторождения гравия:

• Доля камней: — в первом месторождении, — во втором месторождении;
• Доля песка: — в первом месторождении, — во втором месторождении;
• Минимальные количества в тоннах: — камни, песок;
• Стоимости за тонну: — в первом месторождении, — во втором месторождении;
• Необходимые вероятности соблюдения ограничений: — для камней, — для песка.

Давайте решим задачу так, словно никаких случайностей нет и просто возьмем средние значения.

Код для решения и отрисовки

# Средние значения параметров
a11 = (0.6 + 0.7) / 2   # среднее для U[0.6, 0.7]
a12 = (0.75 + 0.9) / 2  # среднее для U[0.75, 0.9]
a21 = 1 - a11
a22 = 1 - a12
c1 = (2.2 + 2.5) / 2    # средняя стоимость 1-го месторождения
c2 = (1.7 + 2.1) / 2    # средняя стоимость 2-го месторождения

# Целевая функция (минимизация стоимости)
c = [c1, c2]

# Ограничения в виде A_ub * x <= b_ub
A = [
    [-a11, -a12],  # -x1*a11 - x2*a12 <= -b1 (камни)
    [-a21, -a22]   # -x1*a21 - x2*a22 <= -b2 (песок)
]
b = [-6500, -3000]

# Границы переменных (x1 >= 0, x2 >= 0)
bounds = [(0, None), (0, None)]

# Решение задачи линейного программирования
res = linprog(c, A_ub=A, b_ub=b, bounds=bounds, method='highs')

# Вывод решения:
print(f"Оптимальное количество из 1-го месторождения: {res.x[0]:.2f} тонн")
print(f"Оптимальное количество из 2-го месторождения: {res.x[1]:.2f} тонн")
print(f"Минимальная стоимость: {res.fun:.2f} у.е.")

# Построим область допустимых решений
x1 = np.linspace(6000, 10000, 500)

# Уравнения ограничений (преобразованные к x2)
# 1) a11*x1 + a12*x2 >= b1 => x2 >= (b1 - a11*x1)/a12
x2_1 = (6500 - a11*x1)/a12
# 2) a21*x1 + a22*x2 >= b2 => x2 >= (b2 - a21*x1)/a22
x2_2 = (3000 - a21*x1)/a22

# Визуализация
plt.figure(figsize=(10, 6), dpi=150)
plt.plot(x1, x2_1, label=r'Камни: $0.65x_{1} + 0.825x_{2} ≥ 6500$')
plt.plot(x1, x2_2, label=r'Песок: $0.35x_{1} + 0.175x_{2} ≥ 3000$')
plt.fill_between(x1, np.maximum(x2_1, x2_2), 10000, alpha=0.1, label='Область допустимых решений')

# Оптимальная точка
plt.scatter(res.x[0], res.x[1], color='red', zorder=100,
            label=f'Оптимум ({res.x[0]:.0f}, {res.x[1]:.0f})')

# Линии постоянной стоимости
for cost in [20000, 21489, 22000]:
    x2_cost = (cost - c1*x1)/c2
    plt.plot(x1, x2_cost, '--', alpha=0.3, label=f'Стоимость={cost}' if cost==21489 else None)

plt.xlim(6000, 10000)
plt.ylim(0, 4000)
plt.xlabel(r'Количество из 1-го месторождения ($x_{1}$), тонн')
plt.ylabel(r'Количество из 2-го месторождения ($x_{2}$), тонн')
plt.title('Оптимизация закупки гравия для фундамента')
plt.legend()
plt.show()

Стоимость оптимального решения составляет 21489 у.е. — запомним это.

Теперь решим задачу в стохастической постановке.

Код для решения и отрисовки

# Параметры распределений
a11_range = [0.6, 0.7]
a12_range = [0.75, 0.9]
c1_range = [2.2, 2.5]
c2_range = [1.7, 2.1]
b1, b2 = 6500, 3000
n_samples = 10000  # Число сценариев Монте-Карло

# Генерация случайных параметров
np.random.seed(42)
a11_samples = np.random.uniform(a11_range[0], a11_range[1], n_samples)
a12_samples = np.random.uniform(a12_range[0], a12_range[1], n_samples)
c1_samples = np.random.uniform(c1_range[0], c1_range[1], n_samples)
c2_samples = np.random.uniform(c2_range[0], c2_range[1], n_samples)
a21_samples = 1 - a11_samples
a22_samples = 1 - a12_samples

# Целевая функция (ожидаемая стоимость)
def expected_cost(x, c1_samples, c2_samples):
    return np.mean(c1_samples * x[0] + c2_samples * x[1])

def objective(x):
    return expected_cost(x, c1_samples, c2_samples)

# Ограничения в виде штрафа (soft constraints)
def constraint1(x):
    return np.mean(a11_samples * x[0] + a12_samples * x[1] >= b1) - 0.7

def constraint2(x):
    return np.mean(a21_samples * x[0] + a22_samples * x[1] >= b2) - 0.9

# Штрафная функция
def penalized_objective(x):
    cost = objective(x)
    penalty = 0
    if constraint1(x) < 0:
        penalty += 1e6 * abs(constraint1(x))
    if constraint2(x) < 0:
        penalty += 1e6 * abs(constraint2(x))
    return cost + penalty

# Область поиска
bounds = [(0, 10000), (0, 10000)]

result = differential_evolution(penalized_objective, bounds)
print(f"Оптимальные количества: x1 = {result.x[0]:.2f}, x2 = {result.x[1]:.2f}")
print(f"Ожидаемая стоимость: {result.fun:.2f} у.е.")

# Генерация линий ограничений и стоимости для случайных параметров
plt.figure(figsize=(12, 8))
for i in range(100):  # Рисуем 100 случайных сценариев
    a11, a12 = a11_samples[i], a12_samples[i]
    a21, a22 = 1 - a11, 1 - a12
    c1, c2 = c1_samples[i], c2_samples[i]
    
    # Линии ограничений
    x1 = np.linspace(6000, 11000, 100)
    x2_stones = (b1 - a11 * x1) / a12
    x2_sand = (b2 - a21 * x1) / a22
    plt.plot(x1, x2_stones, 'C1-', alpha=0.3, label='Камни' if i==0 else None)
    plt.plot(x1, x2_sand, 'C0-', alpha=0.3, label='Песок' if i==0 else None)
    
    # Линии постоянной стоимости
    #x2_cost = (result.fun - c1 * x1) / c2
    #plt.plot(x1, x2_cost, 'g-', alpha=0.05)

# Оптимальная точка
plt.scatter(result.x[0], result.x[1], c='red', s=80, zorder=100,
            label=f'Оптимум ({res.x[0]:.0f}, {res.x[1]:.0f})')

plt.legend()
plt.ylim(0, 5000)
plt.xlabel(r'Количество из 1-го месторождения ($x_{1}$), тонн')
plt.ylabel(r'Количество из 2-го месторождения ($x_{2}$), тонн')
plt.title('Оптимизация закупки гравия для фундамента\n(100 случайных сценариев)')

plt.show()

Стоимость оптимального решения 23222 у.е., что аж на 1733 у.е. больше, чем в «детерминированном» случае. Почему так? Мы заплатили за гарантии того, что необходимого количества камней будет в 80%-х случаев, а песка в 95%-х. В чем легко убедиться, если прокрутить множество сценариев.

Код для отрисовки

x = result.x
f, ax = plt.subplots(1, 3, figsize=(12, 3), dpi=150)
sns.histplot(a11_samples * x[0] + a12_samples * x[1], stat='density', ax=ax[0])
ax[0].axvline(b1, c='C3')
sns.histplot(a21_samples * x[0] + a22_samples * x[1], stat='density', ax=ax[1])
ax[1].axvline(b2, c='C3')
sns.histplot(c1_samples * x[0] + c2_samples * x[1], stat='density', ax=ax[2]);
ax[2].axvline(expected_cost(x, c1_samples, c2_samples), c='C3')

ax[0].set_title('Камни')
ax[1].set_title('Песок')
ax[2].set_title('Стоимость')

ax[1].set_ylabel('')
ax[2].set_ylabel('')

plt.show()

Если есть риск несоблюдения тех.процесса приготовления бетона и его могут забраковать, то лучше заплатить за страховку и использовать SP. Если у нас есть какие-то дополнительные данные и разумные гипотезы, то их тоже можно использовать путем уточнения распределений.

Важно понимать, что SP — это действительно что-то вроде механизма самострахования, который тоже стоит денег. Но этот же самый механизм очень хорошо защищает от неблагоприятных сценариев. Хочешь, чтобы все было «ровно» — плати. Вот такая она, реальность.

Непрерывная оптимизация

Давайте представим, что мы играем в следующую игру. Есть 100 необычных монет — вероятность выпадения орла на каждой отдельной монете может настраиваться. Монеты бросаются последовательно, одна за другой. Если выпадает орел, то мы получаем награду, которая обратнопропорциональна вероятности его выпадения. Если выпадает решка, то мы ничего не получаем. Игра останавливается, если количество выпавших орлов стало равно 50 или если выполнено подбрасывание всех ста монет. Цель состоит в том, чтобы получить наибольшую суммарную награду.

Какая стратегия игры будет оптимальной?

Наивная стратегия

Самый наивный способ заключается в том, чтобы выбрать одинаковую вероятность для всех монет, пусть она будет равна 0.5, тогда мы увидим следующее:

Код

def naive_strategy():
    p = 0.5
    
    heads = bernoulli.rvs(p=p, size=100)
    count_heads = np.cumsum(heads)
    
    if count_heads[-1] < 50:
        total_revard = count_heads[-1] * (1 / 0.5)
        count_steps = 100
    else:
        total_revard = 50 * (1 / 0.5)
        count_steps = np.where(count_heads >= 50)[0][0] + 1
    return total_revard, count_steps

data_naive_strategy = np.array([naive_strategy() for i in range(3000)])

total_reward = np.mean(data_naive_strategy[:, 0])
total_flips = np.mean(data_naive_strategy[:, 1])

f, ax = plt.subplots(1, 2, figsize=(10, 4), dpi=150)
sns.histplot(data_naive_strategy[:, 0], discrete=True, ax=ax[0])
ax[0].axvline(total_reward, c='C3')
sns.histplot(data_naive_strategy[:, 1], discrete=True, ax=ax[1])
ax[1].axvline(total_flips, c='C3')

ax[0].set_title(f'Общая награда (в среднем: {total_reward:.2f})')
ax[1].set_title(f'Количество бросков (в среднем: {total_flips:.2f})')
plt.show()

Жадная стратегия

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

Код

import numpy as np

def greedy_strategy(n_coins=100, max_heads=50):
    remaining_heads = max_heads
    remaining_coins = n_coins
    rewards = []

    for t in range(n_coins):
        if remaining_heads == 0:
            break  # Игра окончена
        p = remaining_heads / remaining_coins
        flip = np.random.rand() < p
        if flip:
            reward = 1 / p
            rewards.append(reward)
            remaining_heads -= 1
        remaining_coins -= 1

    return sum(rewards), t

data_greedy_strategy= np.array([greedy_strategy() for _ in range(10000)])

total_reward = np.mean(data_greedy_strategy[:, 0])
total_flips = np.mean(data_greedy_strategy[:, 1])

f, ax = plt.subplots(1, 2, figsize=(10, 4), dpi=150)
sns.histplot(data_greedy_strategy[:, 0], ax=ax[0])
ax[0].axvline(np.mean(data_greedy_strategy[:, 0]), c='C3')
sns.histplot(data_greedy_strategy[:, 1], discrete=True, ax=ax[1])
ax[1].axvline(np.mean(data_greedy_strategy[:, 1]), c='C3')

ax[0].set_title(f'Общая награда (в среднем: {total_reward:.2f})')
ax[1].set_title(f'Количество бросков (в среднем: {total_flips:.2f})')
plt.show()

И правда стало гораздо лучше. Но нет ли тут какого-то подвоха — почему средняя общаяя награда не равна 100, а всего 99 с копейками?

Стохастическая динамическая стратегия

Жадные стратегии просты, интуитивны, но не учитывают неопределенность будущего. Давайте попробуем динамическую стратегию.

Код

# Параметры игры
n_coins = 100
max_heads = 50

# Таблицы политик и значений
V = np.zeros((n_coins + 2, max_heads + 1))  # V[t][s]
policy = np.zeros((n_coins + 2, max_heads + 1))

# Сетка вероятностей
p_values = np.linspace(0.01, 1.0, 100)

# Обратная прогонка
for t in range(n_coins, 0, -1):
    for s in range(max_heads, -1, -1):
        if s >= max_heads:
            policy[t][s] = 0.0
            continue

        best_value = -np.inf
        best_p = 0.01

        for p in p_values:
            reward_hit = 1.0 / p
            next_s = min(s + 1, max_heads)

            value = p * (reward_hit + V[t + 1][next_s]) + (1 - p) * V[t + 1][s]

            # Штраф за маленькие p
            penalty = 0.01 * (1.0 / p)
            value -= penalty

            if value > best_value:
                best_value = value
                best_p = p

        # Минимальная вероятность
        best_p = max(best_p, 0.01)

        V[t][s] = best_value
        policy[t][s] = best_p

def dynamic_strategy(policy, n_coins=100, max_heads=50):
    total_reward = 0.0
    used_heads = 0
    flips_made = 0

    for t in range(1, n_coins + 1):
        if used_heads >= max_heads:
            break  # Лимит орлов исчерпан

        # Получаем оптимальную вероятность из политики
        p = policy[t][used_heads]

        # Бросок монеты
        flip = np.random.rand() < p
        flips_made += 1

        if flip:
            reward = 1.0 / p
            total_reward += reward
            used_heads += 1

    return total_reward, flips_made

data_dynamic_strategy = np.array([dynamic_strategy(policy, n_coins=100, max_heads=50) for _ in range(10000)])

total_reward = np.mean(data_dynamic_strategy[:, 0])
total_flips = np.mean(data_dynamic_strategy[:, 1])

f, ax = plt.subplots(1, 2, figsize=(10, 4), dpi=150)
sns.histplot(data_dynamic_strategy[:, 0], ax=ax[0])
ax[0].axvline(np.mean(data_dynamic_strategy[:, 0]), c='C3')
sns.histplot(data_dynamic_strategy[:, 1], discrete=True, ax=ax[1])
ax[1].axvline(np.mean(data_dynamic_strategy[:, 1]), c='C3')

ax[0].set_title(f'Общая награда (в среднем: {total_reward:.2f})')
ax[1].set_title(f'Количество бросков (в среднем: {total_flips:.2f})')
plt.show()

Вот теперь средние значения общего выигрыша и количества бросков монеты очень близки к 100.

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

В заключение

Лучшая стратегия — это воспринимать мир таким, какой он есть, и SP очень хорошо вписывается в эту парадигму.

Слишком много всего осталось за кадром. Мало теории, формул, строгих обоснований, но просто оглянитесь вокруг, задумайтесь: все, что нас окружает — это результат чьих-то планов. Сколько МЫ теряем, игнорируя случайность. Жизнь в детерминированном мире — это иллюзия, которая обходится слишком дорого абсолютно всегда и везде от здравоохранения до экономики.