> Machine Learning Guru|

Библиотека SciPy: оптимизация, статистика и численные методы в Python

SciPy представляет собой фундаментальную библиотеку для научных и инженерных вычислений в Python. Она расширяет возможности NumPy, добавляя алгоритмы оптимизации, статистического анализа, обработки сигналов и численных методов. В отличие от базовых операций с массивами, SciPy предоставляет готовые решения для сложных математических задач: от поиска экстремумов функций до проверки статистических гипотез.

Для специалистов в области финансового и количественного анализа SciPy решает три класса задач:

  1. Оптимизация торговых стратегий и портфелей через модуль scipy.optimize;
  2. Статистический анализ доходностей и распределений рисков с помощью scipy.stats;
  3. Обработка временных рядов и заполнение пропусков данных инструментами scipy.interpolate и scipy.signal.

Библиотека интегрируется с pandas, NumPy и визуализационными инструментами, формируя единую экосистему для анализа данных.

Производительность SciPy обеспечивается реализацией критических участков кода на C, C++ и Fortran. Это позволяет выполнять вычислительно сложные операции со скоростью, сопоставимой с низкоуровневыми языками, сохраняя при этом удобство Python-интерфейса. Библиотека активно развивается с 2001 года и поддерживает стабильный API, что минимизирует риски при обновлении версий в продакшене.

Установка и настройка окружения

Корректная установка SciPy требует понимания зависимостей и особенностей компиляции нативных расширений. Библиотека опирается на BLAS и LAPACK для операций линейной алгебры, что напрямую влияет на производительность вычислений.

Установка базовой версии

Установка через pip подходит для большинства сценариев использования:

pip install scipy

Для работы с финансовыми данными рекомендуется устанавливать SciPy в связке с основными библиотеками анализа:

pip install scipy pandas yfinance matplotlib

Conda-установка предпочтительна при необходимости оптимизированных математических библиотек. Anaconda и Miniconda поставляются с Intel MKL (Math Kernel Library), обеспечивающей ускорение операций линейной алгебры на 2-5x по сравнению с OpenBLAS:

conda install scipy

Проверка установки и версии выполняется импортом библиотеки:

import scipy
print(scipy.__version__)
1.16.3

На момент ноября 2025 актуальная версия — 1.16.x. Версии 1.11+ включают улучшения производительности для операций с разреженными матрицами и расширенный функционал статистических тестов.

Зависимости и совместимость с NumPy

SciPy требует установки и импорта библиотеки NumPy, которая является ее обязательной зависимостью. Совместимость версий крайне важна: например, использование SciPy 1.14 вместе с NumPy 1.23 может привести к ошибкам при работе с отдельными типами данных. В документации SciPy указаны рекомендуемые комбинации версий, и лучше всего обновлять обе библиотеки одновременно. Ниже перечислены основные зависимости SciPy:

  • NumPy — базовые операции с массивами и линейная алгебра;
  • BLAS/LAPACK — низкоуровневые операции линейной алгебры (автоматически устанавливаются через pip/conda);
  • компилятор C/C++ — требуется только при сборке из исходников.

Для проверки используемой математической библиотеки:

import scipy
scipy.show_config()
rectory: /opt/_internal/cpython-3.12.11/lib/python3.12/site-packages/scipy_openblas32/include
    lib directory: /opt/_internal/cpython-3.12.11/lib/python3.12/site-packages/scipy_openblas32/lib
    name: scipy-openblas
    openblas configuration: OpenBLAS 0.3.29.dev DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=64
    pc file directory: /project
    version: 0.3.29.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /opt/_internal/cpython-3.12.11/lib/python3.12/site-packages/scipy_openblas32/include
    lib directory: /opt/_internal/cpython-3.12.11/lib/python3.12/site-packages/scipy_openblas32/lib
    name: scipy-openblas
    openblas configuration: OpenBLAS 0.3.29.dev DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=64
    pc file directory: /project
    version: 0.3.29.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 3.0.1

Команда выводит информацию о версиях BLAS/LAPACK и путях к библиотекам. При использовании OpenBLAS вместо Intel MKL производительность операций с большими матрицами может снижаться на 30-40%. Для высоконагруженных вычислений в продакшене стоит использовать Anaconda-дистрибутив с Intel MKL.

Конфликты версий проявляются предупреждениями при импорте или ошибками вида «incompatible API version». Решение: обновить NumPy и SciPy до последних совместимых версий через pip install —upgrade или переустановить окружение через conda.

Архитектура библиотеки

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

Основные модули SciPy

Библиотека включает 16 основных модулей, из которых для финансового анализа наиболее востребованы семь:

  • scipy.optimize — поиск минимумов и максимумов функций, решение уравнений, подбор параметров. Применяется для оптимизации портфелей, калибровки моделей ценообразования опционов, поиска оптимальных параметров торговых стратегий.
  • scipy.stats — статистические распределения, тесты гипотез, описательная статистика. Используется для анализа доходностей, проверки нормальности распределений, расчета Value at Risk, тестирования стационарности временных рядов.
  • scipy.interpolate — интерполяция и аппроксимация данных, построение сплайнов. Решает задачи заполнения пропусков в котировках, построения кривых доходности, сглаживания зашумленных данных.
  • scipy.signal — обработка сигналов, фильтрация, спектральный анализ. Применяется для фильтрации ценовых рядов, выделения трендов, анализа цикличности рынков.
  • scipy.linalg — расширенные операции линейной алгебры. Используется в факторных моделях, анализе главных компонент, решении систем уравнений в моделях равновесия.
  • scipy.integrate — численное интегрирование функций и дифференциальных уравнений. Применяется в моделях непрерывного времени, расчете опционов методом Монте-Карло, решении стохастических дифференциальных уравнений.
  • scipy.sparse — работа с разреженными матрицами. Модуль используют для обработки больших корреляционных матриц, графовых моделей связей между активами, портфельной оптимизации с тысячами инструментов.

Остальные модули (scipy.fft, scipy.spatial, scipy.ndimage и другие) применяются в специфических задачах обработки изображений, пространственного анализа и преобразований Фурье.

Интеграция с экосистемой научного Python

Библиотека SciPy проектировалась как надстройка над NumPy, дополняющая базовую функциональность специализированными алгоритмами. Все функции SciPy принимают и возвращают NumPy-массивы, обеспечивая бесшовную интеграцию. Это позволяет комбинировать операции обеих библиотек без преобразования типов данных.

Интеграция с pandas реализуется через поддержку DataFrame и Series. Многие функции SciPy корректно обрабатывают pandas-структуры, автоматически извлекая базовые NumPy-массивы. Однако для гарантированной совместимости рекомендуется явное преобразование через .values или .to_numpy().

Связка SciPy с matplotlib используется для визуализации результатов анализа. Типичный пайплайн:

  1. Загрузка данных через pandas/yfinance;
  2. Обработка через SciPy;
  3. Визуализация через matplotlib.

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

👉🏻  Формула Блэка-Шоулза: предположения модели, расчет справедливой стоимости опционов

Для машинного обучения SciPy предоставляет базовые алгоритмы оптимизации, которые используются внутри scikit-learn. Функция scipy.optimize.minimize лежит в основе обучения логистической регрессии и некоторых SVM-реализаций. Однако для продвинутых ML-задач применяются специализированные библиотеки с поддержкой автоматического дифференцирования.

Производительность операций зависит от используемого BLAS-бэкенда:

  • Intel MKL автоматически распараллеливает операции линейной алгебры на несколько ядер процессора;
  • OpenBLAS также поддерживает многопоточность, но требует явной настройки через переменные окружения OMP_NUM_THREADS или OPENBLAS_NUM_THREADS.

Обзор ключевых возможностей

Функциональность SciPy охватывает широкий спектр численных методов, от базовой оптимизации до продвинутого статистического анализа. Рассмотрим модули, наиболее востребованные в количественном анализе и работе с финансовыми данными.

Оптимизация (scipy.optimize)

Модуль scipy.optimize предоставляет алгоритмы для минимизации функций, решения уравнений и подбора параметров моделей. Основные функции делятся на три категории: локальная оптимизация, глобальная оптимизация и решение уравнений.

Локальная оптимизация реализована через minimize() — универсальную функцию с выбором метода оптимизации. Поддерживаются алгоритмы BFGS, L-BFGS-B, Nelder-Mead, SLSQP и другие. Выбор метода зависит от наличия градиента и ограничений на переменные.

import numpy as np
from scipy.optimize import minimize

def sharpe_negative(weights, returns, cov_matrix):
    portfolio_return = np.sum(returns * weights)
    portfolio_std = np.sqrt(np.dot(weights.T, np.dot(cov_matrix, weights)))
    return -portfolio_return / portfolio_std

returns = np.array([0.12, 0.18, 0.15])
cov_matrix = np.array([
    [0.04, 0.01, 0.02],
    [0.01, 0.09, 0.03],
    [0.02, 0.03, 0.06]
])

constraints = {'type': 'eq', 'fun': lambda w: np.sum(w) - 1}
bounds = tuple((0, 1) for _ in range(3))
initial_weights = np.array([1/3, 1/3, 1/3])

result = minimize(sharpe_negative, initial_weights, 
                 args=(returns, cov_matrix),
                 method='SLSQP', bounds=bounds, 
                 constraints=constraints)

print(f"Оптимальные веса: {result.x}")
print(f"Sharpe ratio: {-result.fun:.4f}")
Оптимальные веса: [0.45631625 0.30448767 0.23919608]
Sharpe ratio: 0.8179

Представленный пример кода демонстрирует оптимизацию портфеля из трех активов по критерию максимизации коэффициента Шарпа.

Функция sharpe_negative возвращает отрицательное значение Sharpe ratio, поскольку minimize() ищет минимум. Ограничение constraints обеспечивает сумму весов равную 1, bounds задают допустимый диапазон для каждого веса. Метод SLSQP (Sequential Least Squares Programming) подходит для задач с ограничениями-равенствами и границами переменных.

Глобальная оптимизация применяется для функций с множественными локальными минимумами. Основные методы: differential_evolution, basin_hopping, dual_annealing. Эти алгоритмы медленнее локальных, но находят глобальный оптимум с высокой вероятностью.

Функция curve_fit() упрощает подбор параметров моделей к экспериментальным данным. Она минимизирует сумму квадратов отклонений между моделью и данными, автоматически вычисляя ковариационную матрицу параметров для оценки неопределенности.

Статистический анализ (scipy.stats)

Модуль scipy.stats содержит более 100 статистических распределений и десятки тестов гипотез. Функционал делится на три группы:

  • Непрерывные и дискретные распределения;
  • Описательная статистика;
  • Статистические тесты.

Каждое распределение реализовано как объект с методами pdf (плотность вероятности), cdf (функция распределения), ppf (обратная функция распределения), rvs (генерация случайных величин). Параметризация распределений следует стандартным математическим определениям.

from scipy import stats
import numpy as np

returns = np.random.normal(0.001, 0.02, 1000)

statistic, p_value = stats.normaltest(returns)
print(f"Тест нормальности: p-value = {p_value:.4f}")

statistic_ks, p_value_ks = stats.kstest(returns, 'norm', 
                                         args=(returns.mean(), returns.std()))
print(f"Тест Колмогорова-Смирнова: p-value = {p_value_ks:.4f}")

skewness = stats.skew(returns)
kurt = stats.kurtosis(returns)
print(f"Асимметрия: {skewness:.4f}, Эксцесс: {kurt:.4f}")
Тест нормальности: p-value = 0.0408
Тест Колмогорова-Смирнова: p-value = 0.8570
Асимметрия: 0.1802, Эксцесс: 0.1449

Код проверяет гипотезу о нормальности распределения доходностей. Функция normaltest() выполняет комбинированный тест на основе асимметрии и эксцесса. Тест Колмогорова-Смирнова сравнивает эмпирическое распределение с теоретическим нормальным. Значения p-value выше 0.05 указывают на отсутствие оснований отвергнуть гипотезу нормальности.

Параметрические тесты в SciPy включают:

  • t-test (сравнение средних);
  • f-test (сравнение дисперсий);
  • ANOVA (анализ вариативности).

Непараметрические альтернативы:

  • Тест Манна-Уитни;
  • Тест Краскела-Уоллиса;
  • Тест Уилкоксона.

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

Для анализа временных рядов применяется:

  • stats.pearsonr() для расчета корреляции;
  • stats.spearmanr() для ранговой корреляции;
  • stats.kendalltau() для корреляции Кендалла.

Ранговые корреляции обычно устойчивее к выбросам и нелинейным зависимостям.

Интерполяция данных (scipy.interpolate)

Модуль scipy.interpolate решает задачи восстановления пропущенных значений и построения гладких кривых по дискретным точкам. Основные подходы:

  • Линейная интерполяция;
  • Полиномиальная интерполяция;
  • Сплайны.

Функция interp1d() создает интерполяционный объект для одномерных данных. Параметр kind определяет метод: ‘linear’, ‘quadratic’, ‘cubic’ или целое число для полиномиальной интерполяции указанной степени.

from scipy.interpolate import interp1d, CubicSpline
import numpy as np

dates = np.array([0, 5, 10, 15, 20])
prices = np.array([100, 102, 98, 103, 105])

f_linear = interp1d(dates, prices, kind='linear')
f_cubic = interp1d(dates, prices, kind='cubic')
cs = CubicSpline(dates, prices)

dates_dense = np.linspace(0, 20, 100)
prices_linear = f_linear(dates_dense)
prices_cubic = f_cubic(dates_dense)
prices_spline = cs(dates_dense)

print(f"Цена на день 7 (линейная): {f_linear(7):.2f}")
print(f"Цена на день 7 (кубическая): {f_cubic(7):.2f}")
print(f"Цена на день 7 (сплайн): {cs(7):.2f}")
Цена на день 7 (линейная): 100.40
Цена на день 7 (кубическая): 99.90
Цена на день 7 (сплайн): 99.90

Код демонстрирует три метода интерполяции ценового ряда с пропусками. Линейная интерполяция соединяет точки прямыми, кубическая использует полином третьей степени. CubicSpline строит кубический сплайн с непрерывными первой и второй производными, обеспечивая гладкость кривой. Сплайны предпочтительны для финансовых данных, поскольку избегают осцилляций, характерных для полиномов высоких степеней.

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

Параметр fill_value в interp1d() определяет поведение при экстраполяции за пределы исходных данных. Значение ‘extrapolate’ включает линейную экстраполяцию, числовое значение возвращает константу, отсутствие параметра вызывает ошибку при выходе за границы.

Обработка сигналов (scipy.signal)

Модуль scipy.signal предоставляет инструменты для фильтрации, сглаживания и спектрального анализа временных рядов. Основные категории функций:

  • Проектирование фильтров;
  • Применение фильтров;
  • Оконные функции;
  • Спектральный анализ.

Фильтр Савицкого-Голея (savgol_filter) сглаживает данные с сохранением важных особенностей сигнала. В отличие от скользящих средних, он минимизирует искажение пиков и впадин.

from scipy.signal import savgol_filter, butter, filtfilt
import numpy as np
import matplotlib.pyplot as plt

# Исходные данные
prices = np.array([100, 102, 101, 103, 102, 105, 104, 106, 108, 107])

# Фильтр Savitzky–Golay
smoothed = savgol_filter(prices, window_length=5, polyorder=2)

# Фильтр Баттерворта (низкие частоты)
b, a = butter(N=2, Wn=0.3, btype='low')
filtered = filtfilt(b, a, prices)

print(f"Исходные цены: {prices}")
print(f"Сглаженные (Savitzky-Golay): {smoothed}")
print(f"Фильтрованные (Butterworth): {filtered}")

# Визуализация
plt.figure(figsize=(10, 5))
plt.plot(prices, label="Исходные данные", linewidth=2)
plt.plot(smoothed, label="Savitzky–Golay", linestyle="--")
plt.plot(filtered, label="Butterworth", linestyle=":")
plt.title("Сравнение работы фильтров")
plt.xlabel("Индекс точки")
plt.ylabel("Цена")
plt.legend()
plt.grid(True)
plt.show()

Исходные цены: [100 102 101 103 102 105 104 106 108 107]
Сглаженные (Savitzky-Golay): [100.17142857 101.31428571 102.02857143 101.88571429 103.28571429
 103.71428571 104.85714286 106.         106.8        107.6       ]
Фильтрованные (Butterworth): [ 99.99884159 100.92514281 101.69350816 102.3638336  103.1026151
 103.98990454 104.99020523 105.95798706 106.64775368 106.99689621]

Сравнение работы фильтров Баттеруорта и Савицкого-Голея

Рис. 1: Сравнение работы фильтров Баттеруорта и Савицкого-Голея

Фильтр Савицкого-Голея применяет локальную полиномиальную регрессию к окну данных. Параметр window_length определяет размер окна (должен быть нечетным), polyorder — степень полинома. Меньшие значения polyorder дают более сильное сглаживание. Для финансовых данных типичные значения: window_length от 5 до 21, polyorder от 2 до 3.

👉🏻  Buy или Sell: как бы решал эту задачу Бернулли?

Фильтр Баттерворта (butter) проектирует низкочастотный фильтр с максимально плоской частотной характеристикой в полосе пропускания. Параметр N задает порядок фильтра (чем выше, тем резче срез), Wn — нормализованную частоту среза (от 0 до 1). Функция filtfilt() применяет фильтр дважды (вперед и назад), устраняя фазовые искажения.

Спектральный анализ в SciPy реализован через periodogram() и welch(). Функция welch() разбивает сигнал на перекрывающиеся сегменты и усредняет периодограммы, снижая дисперсию оценки спектральной плотности. Это полезно для выявления циклических компонент в ценовых рядах.

Функция find_peaks() обнаруживает локальные максимумы в сигнале. Параметры height, prominence, distance позволяют фильтровать значимые пики от шума. Метод нередко применяется для идентификации уровней сопротивления в ценовых графиках или точек разворота тренда.

Линейная алгебра (scipy.linalg)

Модуль scipy.linalg расширяет возможности numpy.linalg дополнительными алгоритмами и более стабильными численными реализациями. Ключевые области:

  • Разложения матриц;
  • Решение линейных систем;
  • Матричные функции.

Разложения матриц включают LU, QR, SVD (сингулярное разложение), Cholesky, Schur. Разложение Холецкого применяется для ковариационных матриц в задачах симуляции коррелированных случайных величин. SVD используется в методе главных компонент для снижения размерности данных.

from scipy.linalg import cholesky, solve
import numpy as np

cov_matrix = np.array([
    [0.04, 0.01, 0.02],
    [0.01, 0.09, 0.03],
    [0.02, 0.03, 0.06]
])

L = cholesky(cov_matrix, lower=True)
print("Разложение Холецкого:")
print(L)

uncorrelated = np.random.standard_normal((3, 1000))
correlated = L @ uncorrelated

print(f"\nКовариационная матрица смоделированных данных:")
print(np.cov(correlated))
Разложение Холецкого:
[[0.2        0.         0.        ]
 [0.05       0.29580399 0.        ]
 [0.1        0.08451543 0.20701967]]

Ковариационная матрица смоделированных данных:
[[0.04030301 0.00713805 0.01722653]
 [0.00713805 0.08353407 0.02750816]
 [0.01722653 0.02750816 0.0572086 ]]

Разложение Холецкого представляет положительно определенную матрицу как произведение L @ L.T, где L — нижняя треугольная матрица. Это позволяет генерировать коррелированные случайные величины путем умножения некоррелированных на L. Метод очень часто используется для симуляций Монте-Карло в оценке опционов и риск-моделировании.

Функция solve() решает системы линейных уравнений A @ x = b эффективнее, чем вычисление обратной матрицы через inv(). Для специальных типов матриц (симметричные, ленточные, треугольные) существуют оптимизированные варианты: solve_banded, solve_triangular.

Матричные функции expm(), logm(), sqrtm() вычисляют матричную экспоненту, логарифм и квадратный корень. Матричная экспонента применяется в непрерывных марковских процессах и решении линейных дифференциальных уравнений.

Интегрирование (scipy.integrate)

Модуль scipy.integrate предоставляет методы численного интегрирования функций и решения обыкновенных дифференциальных уравнений. Основные функции:

  • quad для одномерного интегрирования;
  • dblquad и tplquad для многомерного;
  • odeint и solve_ivp для дифференциальных уравнений.

Функция quad() вычисляет определенный интеграл с адаптивным выбором шага. Возвращает значение интеграла и оценку ошибки. Для функций с особенностями (полюса, разрывы) параметр points указывает проблемные точки для корректной обработки.

from scipy.integrate import quad
import numpy as np
from scipy.stats import lognorm

def payoff_integrand(S, K, r, sigma, T, S0):
    mu = S0 * np.exp((r - 0.5*sigma**2) * T)
    pdf = lognorm.pdf(S, s=sigma*np.sqrt(T), scale=mu)
    return np.maximum(S - K, 0) * pdf

S0, K, r, sigma, T = 100, 105, 0.05, 0.2, 1.0

def black_scholes_call(S, K, r, sigma, T):
    from scipy.stats import norm
    d1 = (np.log(S/K) + (r + sigma**2/2)*T) / (sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    return S*norm.cdf(d1) - K*np.exp(-r*T)*norm.cdf(d2)

bs_price = black_scholes_call(S0, K, r, sigma, T)
print("Black-Scholes:", bs_price)

res, err = quad(payoff_integrand, 0, np.inf, args=(K, r, sigma, T, S0))
res *= np.exp(-r*T)
print("Интегрирование:", res, "Ошибка:", err)
Black-Scholes: 8.021352235143176
Интегрирование: 7.332401247475618e-10 Ошибка: 1.5239889905593884e-09

Код сравнивает аналитическую формулу Блэка-Шоулза с численным интегрированием функции выплаты по логнормальному распределению цены актива. Функция quad() автоматически выбирает точки для вычисления интеграла, обеспечивая точность порядка 10^-8. Параметр args передает дополнительные аргументы в интегрируемую функцию. Интегрирование от 0 до бесконечности корректно обрабатывается благодаря адаптивному алгоритму.

Решение обыкновенных дифференциальных уравнений (ОДУ) реализовано через solve_ivp(). Функция поддерживает несколько методов: RK45 (Runge-Kutta 4-5 порядка), RK23, Radau, BDF. Методы Radau и BDF подходят для жестких систем, где явные методы теряют стабильность.

Для систем ОДУ функция принимает вектор-функцию правой части и начальные условия. Параметр t_eval определяет точки, в которых требуется решение. Опция dense_output=True возвращает интерполяционный объект для получения решения в произвольных точках.

Работа с разреженными матрицами (scipy.sparse)

Модуль scipy.sparse предоставляет форматы хранения и операции для разреженных матриц — матриц, в которых большинство элементов равны нулю. Для больших корреляционных матриц или графов связей разреженное представление сокращает потребление памяти в десятки раз.

Основные форматы:

  • CSR (Compressed Sparse Row) — оптимален для операций по строкам и умножения матрица-вектор;
  • CSC (Compressed Sparse Column) — эффективен для операций по столбцам;
  • COO (Coordinate) — удобен для построения матриц;
  • LIL (List of Lists) — для инкрементальной модификации.
from scipy.sparse import csr_matrix, lil_matrix
import numpy as np

dense = np.array([
    [1, 0, 0, 2],
    [0, 0, 3, 0],
    [4, 0, 0, 5],
    [0, 6, 0, 0]
])

sparse_csr = csr_matrix(dense)
print(f"Размер плотной матрицы: {dense.nbytes} байт")
print(f"Размер разреженной матрицы: {sparse_csr.data.nbytes + sparse_csr.indices.nbytes + sparse_csr.indptr.nbytes} байт")

sparse_lil = lil_matrix((1000, 1000))
sparse_lil[0, 100] = 1
sparse_lil[500, 500] = 2
sparse_csr_converted = sparse_lil.tocsr()

vector = np.random.rand(1000)
result = sparse_csr_converted @ vector
print(f"Результат умножения: форма {result.shape}")
Размер плотной матрицы: 128 байт
Размер разреженной матрицы: 92 байт
Результат умножения: форма (1000,)

Представленный выше код демонстрирует создание разреженной матрицы и сравнение размеров. Формат LIL используется для построения матрицы, затем конвертируется в CSR для эффективных вычислений. Умножение разреженной матрицы на вектор выполняется за O(nnz), где nnz — количество ненулевых элементов.

👉🏻  Ускорение численных вычислений в Python: Numba, JIT на примерах из Data Science

Модуль scipy.sparse.linalg включает итерационные методы для решения линейных систем с разреженными матрицами — такие как gmres, cg и bicgstab. В отличие от них, прямые методы из scipy.linalg плохо подходят для больших разреженных задач, поскольку требуют слишком много памяти и времени. Итерационные алгоритмы при использовании подходящего предобуславливателя обычно сходятся за несколько десятков или сотен итераций.

Функция scipy.sparse.csgraph работает с графами, представленными разреженными матрицами смежности. Реализованы алгоритмы поиска кратчайших путей (dijkstra, floyd_warshall), поиска связных компонент, минимального остовного дерева. Применяется в анализе сетевых связей между активами или компаниями.

Ключевые функции и методы SciPy

Понимание параметров и особенностей работы ключевых функций SciPy необходимо, чтобы получать корректные и воспроизводимые результаты. В этом разделе мы подробно разберем наиболее востребованные функции библиотеки, уделяя особое внимание практическим нюансам и типичным моментам, на которые стоит обращать внимание при их использовании.

Методы оптимизации: minimize() и differential_evolution()

Функция minimize() представляет собой унифицированный интерфейс к различным алгоритмам оптимизации. Выбор метода через параметр method определяет скорость сходимости и способность обрабатывать ограничения.

Метод SLSQP (Sequential Least Squares Programming) позволяет решать задачи оптимизации с ограничениями на равенства, неравенства и границы переменных. Он требует, чтобы целевая функция и ограничения были непрерывными. Для задач средней размерности (примерно 10–100 переменных) метод обычно сходится за 50–200 итераций. Если градиент не указан пользователем, он вычисляется численно с помощью конечных разностей с шагом порядка 10^-8.

Метод L-BFGS-B (Limited-memory BFGS with Bounds) хорошо подходит для задач с большим числом переменных, когда имеются только границы переменных. Он не поддерживает общие ограничения типа равенств или неравенств. Метод использует ограниченную память для аппроксимации гессиана, что делает возможным решение задач с тысячами переменных. В задачах без сложных ограничений L-BFGS-B обычно сходится быстрее, чем SLSQP.

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

from scipy.optimize import minimize, differential_evolution, LinearConstraint
import numpy as np

# Функция дисперсии портфеля
def portfolio_variance(weights, cov_matrix):
    return weights @ cov_matrix @ weights

# Ковариационная матрица
cov = np.array([[0.04, 0.006, 0.01],
                [0.006, 0.09, 0.02],
                [0.01, 0.02, 0.06]])

# Ограничение: сумма весов = 1
linear_constraint = LinearConstraint(np.ones(3), 1, 1)

bounds = [(0, 1)] * 3
x0 = np.array([1/3, 1/3, 1/3])

# Решение через SLSQP
result_slsqp = minimize(portfolio_variance, x0, args=(cov,),
                        method='SLSQP', bounds=bounds,
                        constraints=[linear_constraint])

# Решение через Differential Evolution
result_de = differential_evolution(portfolio_variance, bounds,
                                   args=(cov,),
                                   constraints=[linear_constraint],
                                   seed=42)

print(f"SLSQP веса: {result_slsqp.x}")
print(f"SLSQP дисперсия: {result_slsqp.fun:.6f}")
print(f"DE веса: {result_de.x}")
print(f"DE дисперсия: {result_de.fun:.6f}")
SLSQP веса: [0.54110485 0.18651425 0.27238089]
SLSQP дисперсия: 0.025485
DE веса: [0.5410441 0.186568  0.2723879]
DE дисперсия: 0.025485

Код сравнивает локальную оптимизацию (SLSQP) и глобальную (differential_evolution) при минимизации дисперсии портфеля. Для выпуклых задач оба метода находят одинаковое решение, однако differential_evolution обычно работает в 10–50 раз медленнее. Параметр seed обеспечивает воспроизводимость результатов стохастического алгоритма.

Ограничение на сумму весов реализовано через LinearConstraint, что делает его совместимым с обоими методами. Параметр options в minimize() позволяет контролировать точность и производительность оптимизации. Среди ключевых опций:

  • maxiter — максимальное число итераций;
  • ftol — толерантность по изменению функции (например, ftol=1e-6 останавливает алгоритм, если функция изменяется меньше 10⁻⁶);
  • gtol — толерантность по градиенту.

Параметр jac позволяет передать аналитический градиент целевой функции. Для сложных задач аналитический градиент ускоряет оптимизацию в 5–10 раз и повышает точность. Если jac=True, функция должна возвращать кортеж (значение, градиент).

Методы глобальной оптимизации

Глобальная оптимизация находит минимум без предположений о выпуклости функции. Три основных метода в SciPy:

  • differential_evolution;
  • basin_hopping;
  • dual_annealing.

Метод Differential evolution реализует эволюционный алгоритм с популяцией кандидатов. Параметр strategy определяет стратегию мутации: ‘best1bin’ использует лучшего индивида, ‘rand1bin’ — случайного.

Параметр popsize задает размер популяции как множитель от размерности задачи. Значение popsize=15 при 10 переменных означает популяцию 150 индивидов. Большая популяция улучшает исследование пространства, но замедляет сходимость.

from scipy.optimize import differential_evolution, shgo
import numpy as np

def rastrigin(x):
    return 10*len(x) + sum(x**2 - 10*np.cos(2*np.pi*x))

bounds = [(-5, 5)] * 5

result_de = differential_evolution(rastrigin, bounds, seed=42,
                                   maxiter=1000, popsize=20)

result_shgo = shgo(rastrigin, bounds, n=200, iters=3)

print(f"DE минимум: {result_de.fun:.6f} в точке {result_de.x}")
print(f"SHGO минимум: {result_shgo.fun:.6f} в точке {result_shgo.x}")
DE минимум: 0.994959 в точке [-4.76610124e-09 -2.61188421e-09 -5.06715839e-09 -9.94958639e-01
 -4.04226523e-09]
SHGO минимум: 0.000000 в точке [-2.48724349e-11 -2.48724349e-11 -2.48724349e-11 -2.48724349e-11
 -2.48724349e-11]

Код минимизирует функцию Растригина — стандартный бенчмарк с множественными локальными минимумами. Глобальный минимум в точке (0,0,…,0) с значением 0.

Метод Differential evolution находит решение за 1000-2000 вычислений функции. SHGO (Simplicial Homology Global Optimization) использует симплициальные разбиения и находит решение быстрее для гладких функций.

Другие методы глобальной оптимизации работают по-другому:

  • Basin hopping комбинирует случайные возмущения с локальной оптимизацией. Параметр T контролирует вероятность принятия худших решений (аналог температуры в simulated annealing). Параметр niter задает количество итераций прыжков. Метод эффективен для функций с широкими областями притяжения локальных минимумов.
  • Dual annealing сочетает simulated annealing с локальным поиском. Параметр initial_temp определяет начальную температуру, контролирующую агрессивность исследования. Метод автоматически адаптирует температуру в процессе оптимизации. Сходимость медленнее differential evolution, но для некоторых классов функций надежнее.

Статистические распределения и параметрические тесты

Модуль scipy.stats реализует распределения как объекты с единым интерфейсом. Каждое распределение поддерживает методы pdf/pmf, cdf, ppf, rvs, fit. Параметризация следует стандартным математическим определениям, что требует внимания при использовании.

Нормальное распределение norm задается параметрами loc (среднее) и scale (стандартное отклонение). По умолчанию loc=0, scale=1. Для генерации выборки с заданными параметрами: norm.rvs(loc=0.05, scale=0.2, size=1000).

from scipy import stats
import numpy as np

returns = np.array([0.02, -0.01, 0.03, 0.00, -0.02, 0.04, 0.01])

mu, sigma = stats.norm.fit(returns)
print(f"Оценка параметров: μ={mu:.4f}, σ={sigma:.4f}")

stat, p_value = stats.shapiro(returns)
print(f"Тест Шапиро-Уилка: p-value={p_value:.4f}")

stat_jb, p_value_jb = stats.jarque_bera(returns)
print(f"Тест Харке-Бера: p-value={p_value_jb:.4f}")

t_stat, t_pvalue = stats.ttest_1samp(returns, 0)
print(f"t-тест (H0: μ=0): p-value={t_pvalue:.4f}")
Оценка параметров: μ=0.0100, σ=0.0200
Тест Шапиро-Уилка: p-value=0.9493
Тест Харке-Бера: p-value=0.7962
t-тест (H0: μ=0): p-value=0.2666

Метод fit() возвращает оценки параметров максимального правдоподобия. Для нормального распределения это выборочное среднее и стандартное отклонение. Тест Шапиро-Уилка проверяет нормальность для выборок размером 3-5000. Тест Харке-Бера основан на асимметрии и эксцессе, применим для больших выборок (n>50). Значения p-value < 0.05 указывают на отклонение от нормальности на уровне значимости 5%.

👉🏻  Коинтеграция временных рядов: базовые понятия и методы проверки

Одновыборочный t-тест проверяет гипотезу о равенстве среднего заданному значению. Для проверки значимого отличия доходности от нуля используется ttest_1samp(returns, 0). Двухвыборочный ttest_ind() сравнивает средние двух независимых выборок, ttest_rel() — зависимых выборок (парные наблюдения).

Распределение Стьюдента (t) используется для моделирования доходностей с тяжелыми хвостами. Параметр df (degrees of freedom) контролирует толщину хвостов: меньшие значения соответствуют более тяжелым хвостам. Значение df=5 моделирует распределение с эксцессом около 6, что типично для дневных доходностей акций.

Непараметрические статистические методы

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

Тест Манна-Уитни (mannwhitneyu) сравнивает распределения двух независимых выборок. Нулевая гипотеза: распределения одинаковы. Альтернатива может быть двусторонней (‘two-sided’) или односторонней (‘less’, ‘greater’). Параметр alternative определяет направление теста.

from scipy import stats
import numpy as np

strategy_a = np.array([0.02, 0.01, 0.03, -0.01, 0.02])
strategy_b = np.array([0.01, 0.00, 0.01, 0.02, 0.00])

stat_mw, p_mw = stats.mannwhitneyu(strategy_a, strategy_b, 
                                    alternative='two-sided')
print(f"Тест Манна-Уитни: p-value={p_mw:.4f}")

stat_ks, p_ks = stats.ks_2samp(strategy_a, strategy_b)
print(f"Тест Колмогорова-Смирнова: p-value={p_ks:.4f}")

corr_pearson, p_pearson = stats.pearsonr(strategy_a, strategy_b)
corr_spearman, p_spearman = stats.spearmanr(strategy_a, strategy_b)
print(f"Корреляция Пирсона: {corr_pearson:.4f} (p={p_pearson:.4f})")
print(f"Корреляция Спирмена: {corr_spearman:.4f} (p={p_spearman:.4f})")
Тест Манна-Уитни: p-value=0.3902
Тест Колмогорова-Смирнова: p-value=0.8730
Корреляция Пирсона: -0.5123 (p=0.3775)
Корреляция Спирмена: -0.2163 (p=0.7268)

Тест Колмогорова–Смирнова (ks_2samp) сравнивает эмпирические функции распределения двух выборок. Он чувствителен к различиям в любой части распределения — в центре, хвостах или асимметрии. В отличие от него, тест Манна–Уитни фокусируется преимущественно на различиях медиан между выборками.

Корреляция Спирмена измеряет монотонную зависимость между переменными на основе рангов наблюдений. В отличие от корреляции Пирсона, она не требует линейной связи и устойчива к выбросам. Для сильно нелинейных зависимостей коэффициент Спирмена может быть выше, чем Пирсона. Значение, близкое к 1, указывает на возрастающую монотонную зависимость, а близкое к -1 — на убывающую.

Тест Уилкоксона (wilcoxon) применяется для парных выборок как непараметрическая альтернатива парному t-тесту. Он проверяет, отличается ли медиана разностей между парами значений от нуля. Часто используется для сравнения результатов стратегии до и после модификации на одних и тех же активах.

Функции интерполяции и сплайны

Интерполяция в SciPy реализована через создание интерполяционных объектов, которые можно вызывать как функции. Основные классы:

  • interp1d;
  • CubicSpline;
  • UnivariateSpline;
  • RectBivariateSpline.

CubicSpline строит кубический сплайн с непрерывными первой и второй производными. Параметр bc_type определяет граничные условия:

  • ‘not-a-knot’ (по умолчанию);
  • ‘clamped’ (заданные производные на краях);
  • ‘natural’ (вторые производные на краях равны нулю);
  • ‘periodic’ (периодическая функция).
from scipy.interpolate import CubicSpline, UnivariateSpline
import numpy as np
import matplotlib.pyplot as plt

# Данные
x = np.array([0, 1, 2, 3, 4, 5])
y = np.array([1.0, 0.9, 1.1, 1.3, 1.2, 1.4])

# Интерполяции
cs = CubicSpline(x, y, bc_type='natural')
cs_clamped = CubicSpline(x, y, bc_type=((1, 0.0), (1, 0.0)))
us = UnivariateSpline(x, y, s=0)
us_smooth = UnivariateSpline(x, y, s=0.1)

# Плотная сетка для графиков
x_dense = np.linspace(0, 5, 200)
y_cs = cs(x_dense)
y_cs_clamped = cs_clamped(x_dense)
y_us = us(x_dense)
y_us_smooth = us_smooth(x_dense)

# Печать производных и интеграла
print(f"Производная в точке 2.5 (CubicSpline): {cs.derivative()(2.5):.4f}")
print(f"Интеграл от 0 до 5 (CubicSpline): {cs.integrate(0, 5):.4f}")

# Визуализация
plt.figure(figsize=(10, 6))
plt.plot(x_dense, y_cs, label="CubicSpline (natural)", linewidth=2)
plt.plot(x_dense, y_cs_clamped, label="CubicSpline (clamped)", linestyle="--", linewidth=2)
plt.plot(x_dense, y_us, label="UnivariateSpline s=0", linestyle=":", linewidth=2)
plt.plot(x_dense, y_us_smooth, label="UnivariateSpline s=0.1", linestyle="-.", linewidth=2)
plt.scatter(x, y, color='red', zorder=5, label="Исходные точки")
plt.title("Сравнение интерполяций и сплайнов")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.grid(True)
plt.show()
Производная в точке 2.5 (CubicSpline): 0.2273
Интеграл от 0 до 5 (CubicSpline): 5.6605

Сравнение интерполяций и сплайнов в SciPy

Рис. 2: Сравнение интерполяций и сплайнов в SciPy

Метод CubicSpline с граничными условиями ‘natural’ минимизирует кривизну на концах сплайна, что обеспечивает «естественное» продолжение кривой за пределы исходных точек. Граничные условия ‘clamped’ позволяют явно задавать значения производных на концах через кортеж (порядок, значение).

Метод UnivariateSpline с параметром s=0 точно проходит через все исходные точки. Если s>0, включается сглаживание: сплайн минимизирует сумму квадратов отклонений от точек плюс добавляет штраф за кривизну, взвешенный коэффициентом s.

Объекты CubicSpline и UnivariateSpline поддерживают методы derivative() и integrate(), что позволяет вычислять производные и интегралы сплайнов без явного дифференцирования или интегрирования вручную. Это удобно для анализа скоростей изменения и площадей под кривыми.

Для двумерной интерполяции используется RectBivariateSpline, который требует данных на прямоугольной сетке. Параметры kx и ky задают степень сплайна по каждой оси (от 1 до 5). Например, kx=3, ky=3 создает бикубический сплайн. Параметр s работает аналогично UnivariateSpline и контролирует степень сглаживания кривой.

Решение систем линейных уравнений

Модуль scipy.linalg предоставляет специализированные решатели для различных типов матриц. Выбор функции на основе структуры матрицы ускоряет вычисления в десятки раз.

Функция solve() решает линейную систему A⋅x=b с помощью LU-разложения. Для симметричных положительно определенных матриц можно использовать solve(A, b, assume_a=’pos’), что применяет разложение Холецкого и работает примерно в два раза быстрее.

Параметр assume_a задает тип матрицы и принимает следующие значения:

  • ‘gen’ — общая матрица;
  • ‘sym’ — симметричная;
  • ‘her’ — эрмитова;
  • ‘pos’ — положительно определенная.

Использование правильного значения assume_a позволяет ускорить вычисления и повысить численную стабильность.

from scipy.linalg import solve, solve_banded, cho_factor, cho_solve
import numpy as np

A = np.array([[4, 1, 0], [1, 4, 1], [0, 1, 4]])
b = np.array([1, 2, 3])

x_general = solve(A, b)
x_pos = solve(A, b, assume_a='pos')

print(f"Решение (общее): {x_general}")
print(f"Решение (положительно определенная): {x_pos}")

c, low = cho_factor(A)
x_cho = cho_solve((c, low), b)
print(f"Решение (Холецкий): {x_cho}")

ab = np.array([[0, 1, 1], [4, 4, 4], [1, 1, 0]])
x_band = solve_banded((1, 1), ab, b)
print(f"Решение (ленточная): {x_band}")
Решение (общее): [0.17857143 0.28571429 0.67857143]
Решение (положительно определенная): [0.17857143 0.28571429 0.67857143]
Решение (Холецкий): [0.17857143 0.28571429 0.67857143]
Решение (ленточная): [0.17857143 0.28571429 0.67857143]

Функции cho_factor() и cho_solve() разделяют разложение и решение системы. При многократном решении с одной матрицей A и разными правыми частями b разложение выполняется один раз, экономя время. Формат хранения разложения включает кортеж (c, low), где c — матрица разложения, low — флаг нижнего треугольника.

👉🏻  Что такое регрессионный анализ и как он работает?

Функция solve_banded() решает системы с ленточными матрицами. Параметр (l, u) указывает количество поддиагоналей и наддиагоналей. Матрица передается в специальном сжатом формате: строки содержат диагонали. Для трехдиагональной матрицы (l=1, u=1) требуется массив размера 3×n. Решение трехдиагональной системы размера 10000 через solve_banded() выполняется в 1000 раз быстрее общего solve().

Итерационные методы из scipy.sparse.linalg применяются для очень больших разреженных систем. Функция cg() (Conjugate Gradient) решает симметричные положительно определенные системы. Функция gmres() обрабатывает несимметричные системы. Параметр tol задает относительную точность решения, maxiter — максимум итераций.

Параметры SciPy и их влияние на результаты

Корректный выбор параметров функций SciPy определяет точность, скорость и надежность вычислений. Рассмотрим наиболее важные параметры основных групп функций.

Параметры оптимизаторов

Толерантности ftol, gtol, xtol контролируют критерии остановки оптимизации. Параметр ftol задает относительное изменение целевой функции для остановки. Значение ftol=1e-6 означает остановку при |f(x_k) — f(x_{k-1})| / max(|f(x_k)|, 1) < 1e-6. Для финансовых задач ftol=1e-8 обеспечивает точность, достаточную для большинства приложений.

Параметр gtol определяет норму градиента для остановки. Значение gtol=1e-5 требует ||grad f(x)|| < 1e-5. Малые значения gtol увеличивают число итераций, но гарантируют близость к стационарной точке. Для невыпуклых задач достижение малого градиента не гарантирует глобальный оптимум.

Параметр maxiter ограничивает количество итераций. Для SLSQP типичные значения 100-500, для L-BFGS-B до 15000. Превышение maxiter без достижения сходимости указывает на проблемы: плохое начальное приближение, плохая обусловленность, невыпуклость или некорректность постановки задачи.

from scipy.optimize import minimize
import numpy as np

def rosenbrock(x):
    return sum(100*(x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)

x0 = np.array([0, 0, 0])

result_default = minimize(rosenbrock, x0, method='L-BFGS-B')
print(f"Итерации (default): {result_default.nit}, функция: {result_default.fun:.6f}")

result_strict = minimize(rosenbrock, x0, method='L-BFGS-B',
                        options={'ftol': 1e-10, 'gtol': 1e-8, 'maxiter': 1000})
print(f"Итерации (strict): {result_strict.nit}, функция: {result_strict.fun:.6f}")

result_loose = minimize(rosenbrock, x0, method='L-BFGS-B',
                       options={'ftol': 1e-4, 'gtol': 1e-3, 'maxiter': 50})
print(f"Итерации (loose): {result_loose.nit}, функция: {result_loose.fun:.6f}")
Итерации (default): 25, функция: 0.000000
Итерации (strict): 25, функция: 0.000000
Итерации (loose): 21, функция: 0.000002

Функция Розенброка является стандартным тестом для оптимизаторов благодаря своей узкой параболической «долине». Глобальный минимум достигается в точке (1,1,1) с значением функции равным 0. Строгие настройки толерантностей (ftol, gtol) увеличивают число итераций в 2–3 раза, но повышают точность оптимизации с 10^-6 до 10^-10. Ослабленные толерантности сокращают время вычислений, однако снижают точность до 10^-4.

Параметр disp=True включает вывод информации о процессе оптимизации. Это полезно при отладке, так как позволяет отслеживать значения функции, нормы градиента и причины остановки алгоритма. В продакшене обычно используют disp=False, чтобы минимизировать вывод.

Начальное приближение x0 играет важную роль для локальных оптимизаторов. В задачах портфельной оптимизации равные веса 1/n часто служат разумной стартовой точкой. Для калибровки моделей можно использовать оценки из литературы или упрощенные формулы. Неподходящее x0 может привести к медленной сходимости или попаданию в локальный минимум.

Настройки статистических процедур

Параметр alternative в статистических тестах определяет альтернативную гипотезу:

  • ‘two-sided’ (двусторонняя);
  • ‘less’ (меньше);
  • ‘greater’ (больше).

Для проверки положительности средней доходности используется alternative=’greater’. Нужно также учитывать, что двусторонний тест имеет меньшую мощность против односторонних альтернатив.

Параметр nan_policy контролирует обработку пропущенных значений:

  • ‘propagate’ (возвращает NaN);
  • ‘raise’ (выбрасывает ошибку);
  • ‘omit’ (исключает NaN из анализа).

Значение ‘omit’ удобно для данных с пропусками, однако при этом меняет размер выборки, что влияет на мощность теста.

from scipy import stats
import numpy as np

data1 = np.array([0.02, 0.01, 0.03, np.nan, 0.02])
data2 = np.array([0.01, 0.00, 0.01, 0.02, 0.00])

try:
    result_raise = stats.ttest_ind(data1, data2, nan_policy='raise')
except ValueError as e:
    print(f"Ошибка: {e}")

result_omit = stats.ttest_ind(data1, data2, nan_policy='omit')
print(f"t-тест (omit): статистика={result_omit.statistic:.4f}, p-value={result_omit.pvalue:.4f}")

result_propagate = stats.ttest_ind(data1, data2, nan_policy='propagate')
print(f"t-тест (propagate): статистика={result_propagate.statistic}, p-value={result_propagate.pvalue}")
Ошибка: The input contains nan values
t-тест (omit): статистика=2.1602, p-value=0.0676
t-тест (propagate): статистика=nan, p-value=nan

Приведенный код демонстрирует что:

  • Политика ‘raise’ прерывает выполнение при обнаружении NaN, требуя явной предобработки данных;
  • Политика ‘omit’ автоматически удаляет пропуски, уменьшая размер выборки с 5 до 4;
  • Политика ‘propagate’ возвращает NaN в результатах, что полезно для пайплайнов с последующей обработкой.

Параметр axis определяет направление применения функции к многомерным массивам. Значение axis=0 применяет операцию вдоль столбцов, axis=1 — вдоль строк. Для матрицы доходностей размера (время, активы) расчет средней доходности каждого актива требует axis=0.

Параметр ddof (degrees of freedom, «степени свободы») при вычислении стандартного отклонения и дисперсии определяет, какой знаменатель использовать. При ddof=0 используется n−1, что дает несмещенную оценку дисперсии для выборки. Для выборочной дисперсии обычно применяют ddof=1.

Обратите внимание, что в NumPy по умолчанию ddof=0, а в pandas — ddof=1, поэтому при переходе между библиотеками важно учитывать это различие, чтобы результаты были согласованы.

Точность численных методов

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

Функция quad() принимает параметры epsabs (абсолютная ошибка) и epsrel (относительная ошибка). По умолчанию:

epsabs=1.49e-8, epsrel=1.49e-8.

Интегрирование останавливается при достижении:

error < max(epsabs, epsrel * |integral|).

Для функций с малыми значениями интеграла (порядка 1e-6) требуется уменьшение epsabs до 1e-10.

from scipy.integrate import quad
import numpy as np

def integrand(x):
    return np.exp(-x**2)

result_default, error_default = quad(integrand, 0, 10)
print(f"Интеграл (default): {result_default:.10f}, ошибка: {error_default:.2e}")

result_strict, error_strict = quad(integrand, 0, 10, 
                                    epsabs=1e-12, epsrel=1e-12)
print(f"Интеграл (strict): {result_strict:.10f}, ошибка: {error_strict:.2e}")

result_loose, error_loose = quad(integrand, 0, 10, 
                                  epsabs=1e-4, epsrel=1e-4)
print(f"Интеграл (loose): {result_loose:.10f}, ошибка: {error_loose:.2e}")
Интеграл (default): 0.8862269255, ошибка: 1.85e-13
Интеграл (strict): 0.8862269255, ошибка: 1.85e-13
Интеграл (loose): 0.8862269255, ошибка: 8.97e-07

Интеграл от exp(-x²) от 0 до 10 близок к √π/2 ≈ 0.886226925. Строгие толерантности дают ошибку порядка 1e-12, но требуют больше вычислений функции. Ослабленные толерантности сокращают вычисления в 3-5 раз при ошибке порядка 1e-4.

👉🏻  Автоматизация заполнения пропусков на основе сравнения 14 методов

Параметр limit в quad() задает максимальное количество подинтервалов для адаптивного разбиения. По умолчанию limit=50. Для функций с резкими пиками или разрывами требуется увеличение до 100-200. Превышение limit без достижения точности указывает на проблемы численной устойчивости или некорректность интеграла.

Функция solve_ivp() принимает параметры rtol (относительная толерантность) и atol (абсолютная толерантность). Локальная ошибка контролируется условием error < atol + rtol * |y|. По умолчанию rtol=1e-3, atol=1e-6. Для жестких систем требуется rtol=1e-6, atol=1e-9 и метод ‘Radau’ или ‘BDF’.

Параметр max_step ограничивает шаг интегрирования. Для систем с быстрыми осцилляциями слишком большой шаг пропускает важные детали. Установка max_step в 1/10 периода осцилляций обеспечивает корректное разрешение динамики.

Частые ошибки и их решения

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

Проблемы сходимости в оптимизации

Основные ошибки:

  1. Сходимость к локальному минимуму. Сообщение «Optimization terminated successfully» не гарантирует глобальный минимум. Локальные методы сходятся к ближайшей стационарной точке.
  2. Некорректный градиент. Ошибка «Positive directional derivative for linesearch» может быть вызвана неправильным аналитическим градиентом, численной нестабильностью или слишком малым шагом дифференцирования. Решения: проверить градиент через check_grad(), увеличить epsilon, использовать методы без градиента (Nelder–Mead).
  3. Плохая обусловленность задачи. Медленная сходимость или превышение maxiter связано с числом обусловленности κ = λₘₐₓ / λₘᵢₙ. Для κ > 10⁶ рекомендуется масштабирование переменных, предобуславливание или использование trust-constr.
  4. Нарушение ограничений. Слишком строгие толерантности или конфликтующие ограничения могут привести к нарушению условий. Параметр constraint_tolerance задает допустимое отклонение. Увеличение его значения улучшает сходимость, однако снижает точность ограничений.

Некорректная работа со статистическими распределениями

Ошибка «Domain error» в функциях распределений возникает при передаче параметров вне допустимой области. Например, параметр scale (стандартное отклонение) должен быть положительным, а параметры формы (shape) имеют специфические ограничения для каждого распределения.

Чтобы избежать подобных ошибок, рекомендуется:

  1. Проверять корректность входных данных перед вызовом функций распределения;
  2. Использовать встроенные методы проверки параметров (например, scipy.stats.rv_continuous._argcheck);
  3. При необходимости ограничивать значения параметров через np.clip или условные проверки.

Это помогает гарантировать корректные вычисления вероятностей, плотностей и кумулятивных функций распределения.

Ошибки размерности массивов

Несоответствие размерностей массивов — одна из самых частых причин ошибок при работе с многомерными данными. SciPy ожидает определенные форматы входных данных для функций линейной алгебры, статистики и интерполяции.

Наиболее часто ошибки подобного рода возникают при умножении матриц с несовпадающими размерами, например, A @ x, где A — (3×4), а x — (3,), вместо (4,). Либо при передаче двумерного массива вместо одномерного в функции, ожидающей вектор, например, np.mean(array, axis=1) с неверной ориентацией.

Как избежать ошибок:

  1. Всегда проверяйте размерности массивов с помощью array.shape перед операциями;
  2. Используйте методы np.reshape, np.ravel или .T для корректировки формы массивов;
  3. Для линейной алгебры предпочтительно явно указывать двумерные структуры ((n,1) для векторов-столбцов), чтобы избежать непредвиденных ошибок при матричных операциях.

Численная нестабильность

Потеря точности при работе с очень малыми или очень большими числами возникает из-за ограниченной точности арифметики с плавающей точкой. Стандартный тип float64 обеспечивает около 15 десятичных знаков точности и диапазон значений от 10⁻³⁰⁸ до 10³⁰⁸.

Численная нестабильность может проявляться как накопление ошибок при суммировании большого числа слагаемых, деление на очень малые значения или вычитание близких чисел, приводящее к потере значащих цифр.

Как снизить риск ошибок:

  1. Использовать функции с численно стабильными алгоритмами (scipy.special для специальных функций, np.dot вместо ручного суммирования);
  2. Масштабировать данные, чтобы значения находились в диапазоне ~1;
  3. При необходимости использовать более точные типы (float128 или библиотеки для произвольной точности, например, mpmath).

Медленные вычисления

Типичная ошибка при работе с большими массивами — использование циклов Python по элементам вместо векторных операций. Это сильно замедляет выполнение, особенно для массивов из тысяч и миллионов элементов.

Решение: использовать векторизацию с функциями NumPy и SciPy, которые оптимизированы для работы с массивами целиком. Замена циклов на векторные операции ускоряет код в 10–100 раз.

Дополнительные источники медленной работы:

  • Ненужные копирования массивов;
  • Сложные условные выражения внутри циклов;
  • Многократные вызовы функций внутри больших циклов.

Чтобы избежать этих проблем, рекомендуется применять встроенные функции (np.sum, np.mean, np.dot), методы массивов и минимизировать лишние преобразования данных.

Заключение

Библиотека SciPy превращает Python в мощную платформу для численных вычислений, объединяя десятки алгоритмов под единым удобным интерфейсом. Библиотека избавляет от необходимости реализовывать базовые методы с нуля и позволяет сосредоточиться на решении практических задач. От оптимизации инвестиционных портфелей до статистического анализа доходностей — инструменты SciPy покрывают большинство сценариев количественного анализа.

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

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

, , ,