Введение в NumPy (2)

Пожалуй, первое с чего стоит начать, так это с того, что массивы NumPy могут быть обычными операндами в математических выражениях:

import numpy as np

a = [1, 2, 3]    #  список Python
b = np.array([1, 2, 3])    #  массив NumPy
c = 7

#  Если мы умножим список на число, то...
a * c

[1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3]

#  Если попытаемся к списку это число прибавить, то Python попытается выполнить конкатенацию, а не сложение.

#  Теперь умножим массив NumPy на число:
b * c

array([ 7, 14, 21])

#  Прибавим к массиву число:
b + c

array([ 8,  9, 10])

Для выполнения таких операций на Python, мы были вынуждены писать циклы.

Писать такие циклы в NumPy, нет никакой необходимости, потому что все операции и так выполняются поэлементно:

a = np.array([[5, 7], [11, 13]])
a / 3 #  обычное деление

array([[1.66666667, 2.33333333],
       [3.66666667, 4.33333333]])

a // 3 #  целочисленное деление

array([[1, 2],
       [3, 4]])

a % 3 #  остаток от деления

array([[2, 1],
       [2, 1]])

a ** 3 #  возведение в степень

array([[ 125,  343],
       [1331, 2197]])

1 / a #  частное 1 и каждого элемента массива

array([[0.2       , 0.14285714],
       [0.09090909, 0.07692308]])

-a #  изменение знака элементов массива

array([[ -5,  -7],
       [-11, -13]])

Точно так же обстоят дела и с математическими функциями:

a = np.arange(6)
a

array([0, 1, 2, 3, 4, 5])

np.sin(a)    #  синус каждого элемента массива

array([ 0.        ,  0.84147098,  0.90929743,  0.14112001, -0.7568025 ,
       -0.95892427])

np.log(a)    #  натуральный логарифм элементов массива

/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:1: RuntimeWarning: divide by zero encountered in log
  """Entry point for launching an IPython kernel.

array([      -inf, 0.        , 0.69314718, 1.09861229, 1.38629436,
       1.60943791])

Такие операции как +=, -=, *=, /= и прочие подобные, могут применяться к массивам и так же выполняются поэлементно.

Они не создают новый массив, а изменяют старый:

a = np.array([1, 2, 3, 4, 5])
b = np.array([0.1, 0.2, 0.3, 0.4, 0.5])

a += 2
a

array([3, 4, 5, 6, 7])

a *= 2
a

array([ 6,  8, 10, 12, 14])

#  Вещественный тип ('float64') не может быть
#  преобразован в целочисленный ('int32'):
# a += b

#  А вот преобразование целочисленного типа в вещественный возможно
b += a
b

array([ 6.1,  8.2, 10.3, 12.4, 14.5])

При работе с массивами разного типа, тип результирующего массива приводится к более общему:

a = np.arange(5)
a

array([0, 1, 2, 3, 4])

b = np.linspace(0, 5, 5)
b

array([0.  , 1.25, 2.5 , 3.75, 5.  ])

a.dtype

dtype('int64')

b.dtype

dtype('float64')

c = a + b
c

array([0.  , 2.25, 4.5 , 6.75, 9.  ])

c.dtype

dtype('float64')

Применение логических операций к массивам, так же возможно и так же выполняется поэлементно.

Результатом таких операций является массив булевых значений (True и False):

a = np.array([2, 3, 5, 7, 11, 13])
a > 5

array([False, False, False,  True,  True,  True])

a == 7

array([False, False, False,  True, False, False])

b = np.array([2, 2, 5, 5, 11, 11])
a > b

array([False,  True, False,  True, False,  True])

a == b

array([ True, False,  True, False,  True, False])

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

a = np.array([1, 2, 3])
(a + 3) * 7

array([28, 35, 42])

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

a = np.array([1, 2, 3])
b = np.array([3, 2, 1])
a + b

array([4, 4, 4])

a ** b

array([1, 4, 3])

a = np.arange(9).reshape(3, 3)
a

array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

b = np.arange(9, 0, -1).reshape(3, 3)
b

array([[9, 8, 7],
       [6, 5, 4],
       [3, 2, 1]])

a + b

array([[9, 9, 9],
       [9, 9, 9],
       [9, 9, 9]])

a ** b

array([[   0,    1,  128],
       [ 729, 1024,  625],
       [ 216,   49,    8]])

Хотя, если честно, их размеры должны быть не равны, а должны быть совместимыми.

Если их размеры совместимы, т.е. один массив может быть растянут до размеров другого, то в дело включается механизм транслирования массивов NumPy.

Этот механизм очень прост, но имеет весьма специфичные нюансы.

Рассмотрим простой пример:

a = np.arange(9).reshape(3, 3)
a

array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

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

array([1, 2, 3])

c = np.array([[1, 2, 3], [1, 2, 3], [1, 2, 3]])
c

array([[1, 2, 3],
       [1, 2, 3],
       [1, 2, 3]])

a * b

array([[ 0,  2,  6],
       [ 3,  8, 15],
       [ 6, 14, 24]])

a * c

array([[ 0,  2,  6],
       [ 3,  8, 15],
       [ 6, 14, 24]])

В данном примере массив b может быть растянут до размеров массива a и станет абсолютно идентичен массиву c.

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

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

Вычисление суммы всех элементов в массиве и прочие унарные операции в NumPy реализованы как методы класса ndarray:

a = np.arange(10)
a

array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

a.sum()

45

a.min()

0

a.max()

9

По умолчанию, эти операции применяются к массиву, как к обычному списку чисел, без учета его ранга (размерности).

Но если указать в качестве параметра одну из осей axis, то вычисления будут производиться именно по ней:

b = np.arange(16).reshape(4, 4)
b

array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

b.sum(axis=0)    #  Сумма элементов каждого столбца

array([24, 28, 32, 36])

b.sum(axis=1)    #  Сумма элементов каждой строки

array([ 6, 22, 38, 54])

b.min(axis=1)    #  Минимальный элемент каждой строки

array([ 0,  4,  8, 12])

b.max(axis=0)    #  Максимальный элемент каждого столбца

array([12, 13, 14, 15])

Значения -inf, inf и nan

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

Убедимся в этом еще раз:

np.log(0)

/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:1: RuntimeWarning: divide by zero encountered in log
  """Entry point for launching an IPython kernel.

-inf

Более того, в NumPy мы даже можем делить на ноль:

a = np.array([0])
1 / a

/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:2: RuntimeWarning: divide by zero encountered in true_divide
  

array([inf])

NumPy предупредил нас о том, что встретил деление на ноль, но тем не менее выдал ответ inf (плюс бесконечность).

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

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

В NumPy есть еще одно специальное значение - nan. Данное значение выдается тогда, когда результат вычислений не удается определить:

a = np.array([0, 1, np.inf])
a

array([ 0.,  1., inf])

np.cos(a)

/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in cos
  """Entry point for launching an IPython kernel.

array([1.        , 0.54030231,        nan])

Заметьте, что NumPy нас просто предупредил о том, что ему попалось недопустимое значение, но ошибки не возникло.

Дело в том, что в реальных вычислениях значения nan, inf или -inf встречается очень часто, поэтому появление этого значения проще обрабатывать специальными методами (функции numpy.isnan() и numpy.isinf()), чем постоянно лицезреть сообщения об ошибках.

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

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

np.cos(0)/np.sin(0)

/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:1: RuntimeWarning: divide by zero encountered in double_scalars
  """Entry point for launching an IPython kernel.

inf

np.sin(np.pi/2)/np.cos(np.pi/2)    #  ожидаем значение 0, но...

1.633123935319537e+16

Число 1.633123935319537e+16 появилось потому что в NumPy выполняются арифметические, а не символьные вычисления, т. е. число π хранится в памяти компьютера не как знание о том, что это математическая константа с бесконечным количеством десятичных знаков после запятой, а как обычное число с десятичной точкой (десятичная дробь) равная числу π с очень маленькой, но все же, погрешностью:

np.pi    #  значение числа pi в NumPy

3.141592653589793

Если из-за самых незначительных погрешностей вычисления все же возможны, то NumPy их обязательно выполнит.

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

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

Если вы решили отправиться в самые дебри теории чисел, алгебры и криптографии, то лучшим решением окажется программа GAP.

Программа GAP не является программой Python, но имеет Python интерфейс в замечательной программе Sage, которая определенно заслуживает вашего внимания (см. книгу).

Линейная алгебра

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

a = np.array([1, 2])
b = np.array([3, 4])
np.dot(a, b)

11

Произведение двумерных массивов по правилам линейной алгебры также возможно:

a = np.arange(2, 6).reshape(2, 2)
a

array([[2, 3],
       [4, 5]])

b = np.arange(6, 10).reshape(2, 2)
b

array([[6, 7],
       [8, 9]])

np.dot(a, b)

array([[36, 41],
       [64, 73]])

При этом размеры матриц (массивов) должны быть либо равны, а сами матрицы квадратными, либо быть согласованными, т.е. если размеры матрицы А равны [m,k], то размеры матрицы В должны быть равны [k,n]:

a = np.arange(2, 8).reshape(2, 3)
a

array([[2, 3, 4],
       [5, 6, 7]])

b = np.arange(4, 10).reshape(3, 2)
b

array([[4, 5],
       [6, 7],
       [8, 9]])

np.dot(a,b)

array([[ 58,  67],
       [112, 130]])

Также по правилам умножения матриц, мы можем умножить матрицу на вектор (одномерный массив).

При этом в таком умножении вектор столбец должен находиться справа, а вектор строка слева:

a = np.array([1, 2, 3])
a

array([1, 2, 3])

b = np.arange(4, 10).reshape(3, 2)
b

array([[4, 5],
       [6, 7],
       [8, 9]])

np.dot(a, b)

array([40, 46])

a = np.arange(1, 3).reshape(2, 1)
a

array([[1],
       [2]])

b = np.arange(4, 10).reshape(3, 2)
b

array([[4, 5],
       [6, 7],
       [8, 9]])

np.dot(b, a)

array([[14],
       [20],
       [26]])

Квадратные матрицы можно возводить в степень n т.е. умнажать сами на себя n раз:

a = np.arange(1, 5).reshape(2, 2)
a

array([[1, 2],
       [3, 4]])

np.dot(a, a)    #  Равносильно a**2

array([[ 7, 10],
       [15, 22]])

np.linalg.matrix_power(a, 2)

array([[ 7, 10],
       [15, 22]])

np.linalg.matrix_power(a, 0)

array([[1, 0],
       [0, 1]])

Довольно часто приходится вычислять ранг матриц:

a = np.arange(1, 10).reshape(3, 3)
a

array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])

np.linalg.matrix_rank(a)

2

b = np.arange(1, 24, 2).reshape(3, 4)
b

array([[ 1,  3,  5,  7],
       [ 9, 11, 13, 15],
       [17, 19, 21, 23]])

np.linalg.matrix_rank(b)

2

Еще чаще приходится вычислять определитель матриц, хотя результат вас может немного удивить:

a = np.array([[1, 3], [4, 3]])
a

array([[1, 3],
       [4, 3]])

np.linalg.det(a)

-8.999999999999998

1*3 - 3*4    #  Результат должен быть целым числом

-9

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

Это связано с тем, что алгоритм вычисления определителя использует LU-разложение - это намного быстрее, чем обычный алгоритм, но за скорость все же приходится немного заплатить ручным округлением (конечно, если таковое требуется):

np.linalg.det(a)

-8.999999999999998

round(np.linalg.det(a))

-9

b = np.arange(1, 48, 3).reshape(4, 4)
b

array([[ 1,  4,  7, 10],
       [13, 16, 19, 22],
       [25, 28, 31, 34],
       [37, 40, 43, 46]])

np.linalg.det(b)

0.0

round(np.linalg.det(b))

0

Транспонирование матриц:

a

array([[1, 3],
       [4, 3]])

a.T

array([[1, 4],
       [3, 3]])

Вычисление обратных матриц:

a

array([[1, 3],
       [4, 3]])

b = np.linalg.inv(a)
b

array([[-0.33333333,  0.33333333],
       [ 0.44444444, -0.11111111]])

np.dot(a, b)

array([[1., 0.],
       [0., 1.]])

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

#  система из двух линейных уравнений:
#  1*x1 + 5*x2 = 11
#  2*x1 + 3*x2 = 8
a = np.array([[1, 5], [2, 3]])
b = np.array([11, 8])

x = np.linalg.solve(a, b)
x

array([1., 2.])

np.dot(a, x)

array([11.,  8.])

Статистика

a = np.random.randint(20, size=(5, 5))
a

array([[ 2, 17,  7,  6,  4],
       [ 2,  2,  1,  7, 13],
       [17, 13,  0,  4, 15],
       [11, 19, 12, 12, 14],
       [13,  1, 12, 19,  2]])

np.amin(a)    #  Минимальный элемент массива

0

np.amax(a)    #  максимальный элемент

19

np.amin(a, axis=0)  #  минимальный элемент вдоль первой оси (столбцы)

array([2, 1, 0, 4, 2])

np.amin(a, axis=1)  #  минимальный элемент вдоль второй оси (строки)

array([ 2,  1,  0, 11,  1])

#  Процентили:
np.percentile(a, 25)

2.0

np.percentile(a, 50)

11.0

np.percentile(a, 75)

13.0

Средние значения элементов массива и их отклонения:

a = np.random.randint(13, size=(5, 5))
a

array([[ 2, 12,  0,  4, 12],
       [ 5,  0, 10,  1,  1],
       [ 5,  4,  3,  4, 12],
       [12,  3, 10,  7,  5],
       [ 6,  5,  3,  0,  4]])

np.median(a)    #  медиана элементов массива

4.0

np.mean(a)    #  среднее арифметическое

5.2

np.var(a)    #  дисперсия

15.28

np.std(a)    #  стандартное отклонение

3.9089640571384128

Корреляционные коэфициенты и ковариационные матрицы величин:

x = np.array([1, 4, 3, 7, 10, 8, 14, 21, 20, 23])
y = np.array([4, 1, 6, 9, 13, 11, 16, 19, 15, 22])
z = np.array([29, 22, 24, 20, 18, 14, 16, 11, 9, 10])

#  Линейный коэфициент корреляции Пирсона
#  величин 'x' и 'y'

XY = np.stack((x, y), axis=0)
XY

array([[ 1,  4,  3,  7, 10,  8, 14, 21, 20, 23],
       [ 4,  1,  6,  9, 13, 11, 16, 19, 15, 22]])

np.corrcoef(XY)

array([[1.        , 0.93158099],
       [0.93158099, 1.        ]])

#  Кросс-корреляции:
np.correlate(x, y)

array([1736])

np.correlate(x, z)

array([1486])

#  Ковариационные матрицы:
np.cov(XY)

array([[63.65555556, 49.82222222],
       [49.82222222, 44.93333333]])

np.cov(x)

array(63.65555556)

Так же NumPy предоставляет функции для вычисления гистограмм наборов данных различной размерности и некоторые другие статистичские функции.

Генерация случайных значений

Получение простых случайных данных:

np.random.rand()    #  случайное число от 0 до 1

0.17448219339797233

np.random.rand(10)    #  одномерный массив случайных значений

array([0.70252469, 0.67976413, 0.33093102, 0.71896991, 0.50738282,
       0.85770526, 0.54996816, 0.20724896, 0.69961605, 0.72648606])

np.random.rand(3, 4)    #  двумерный массив случайных значений

array([[0.09784915, 0.92849649, 0.57973627, 0.10194581],
       [0.21016446, 0.24861077, 0.3003258 , 0.1775901 ],
       [0.98115457, 0.64760134, 0.83332208, 0.09634793]])

np.random.randn(10)    #  случайные значения с нормальным распределением

array([-0.18961292, -0.92561458,  0.04817766, -0.12435168,  0.76708842,
       -0.14950132,  1.56869699, -0.19216652,  0.73780653, -0.90962697])

np.random.randint(10)    #  случайное целое число от 0 до 10

0

np.random.randint(10, 100)    #  случайное целое число от 10 до 100

79

np.random.randint(10, size=7)    #  одномерный массив случайных целых чисел

array([9, 0, 3, 4, 7, 3, 4])

np.random.randint(10, size=(4, 4))    #  двумерный массив случайных целых чисел

array([[9, 8, 4, 1],
       [2, 5, 2, 8],
       [6, 7, 4, 6],
       [4, 0, 3, 2]])

Перестановки:

x = np.arange(7)
x

array([0, 1, 2, 3, 4, 5, 6])

np.random.shuffle(x)    #  перетасовывает содержимое массива
x

array([3, 4, 2, 6, 1, 0, 5])

np.random.permutation(7)    #  выполняет тоже самое, не требуя входного массива

array([2, 5, 0, 6, 3, 4, 1])

y = np.arange(9).reshape(3, 3)
y

array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

np.random.shuffle(y)    #  перестановки выполняются только по 1-й оси

NumPy предоставляет порядка 30 функций, позволяющих генерировать случайные числа с самыми разными вероятностными распределениями:

np.random.beta(0.1, 0.6, size=5)    #  бета распределение

array([2.84742109e-02, 1.01824328e-17, 9.25048785e-04, 3.64882719e-01,
       3.25493988e-01])

np.random.gamma(shape=0.8, scale=1.7, size=5)    # гамма распределение

array([0.10021983, 0.16943748, 0.05614188, 0.55054814, 0.86767598])

np.random.pareto(3.5, size=5)    #  Паретто распределение

array([1.46231788, 0.46808183, 0.1089563 , 1.08172899, 0.04363823])

np.random.chisquare(2.2, size=5)     #  хи-квадрат распределение

array([ 0.21156785,  4.38347729, 10.16896532,  0.12034952,  4.05578084])

Вы так же имеете доступ к состоянию генератора случайных чисел, а так же можете управлять им:

#np.random.get_state()    #  Вы можете узнать состояние генератора

np.random.seed(123)    #  устанавливает состояние генератора

np.random.rand(5)

array([0.69646919, 0.28613933, 0.22685145, 0.55131477, 0.71946897])

Множества

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

x = np.array([1, 2, 1, 2, 3, 1])
np.unique(x)

array([1, 2, 3])

np.unique(x, return_counts=True)    #  количество вхождений

(array([1, 2, 3]), array([3, 2, 1]))

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

x = np.array([[0, 1, 1, 2], [0, 1, 1, 2], [9, 1, 1, 2]])
x

array([[0, 1, 1, 2],
       [0, 1, 1, 2],
       [9, 1, 1, 2]])

np.unique(x, axis=0)    #  множество уникальных строк

array([[0, 1, 1, 2],
       [9, 1, 1, 2]])

np.unique(x, axis=1)    #  множество уникальных столбцов

array([[0, 1, 2],
       [0, 1, 2],
       [9, 1, 2]])

Так же имеется ряд других полезных функций:

X = np.array([0, 2, 4, 6, 8])
Y = np.array([0, 3, 4, 6, 7])

np.in1d(X, Y)    #  наличие элементов из X среди элементов Y

array([ True, False,  True,  True, False])

np.intersect1d(X, Y)    #  пересечение множеств элементов массивов

array([0, 4, 6])

np.setdiff1d(X, Y)    #  разность множеств

array([2, 8])

np.setxor1d(X, Y)    #  симметрическая разность (исключающее или)

array([2, 3, 7, 8])

np.union1d(X, Y)     #  объединение множеств

array([0, 2, 3, 4, 6, 7, 8])

Логические операции

Логические функции NumPy, условно, можно разделить на два множества: первое - позволяет определять специфику элементов массива и самого массива; второе - обычные логические операции, которые действуют над массивами поэлементно.

Иногда, возникает потребность определить тип элементов:

A = np.array([-np.inf, np.inf, np.nan, -1, 0, 1, 1.47, 2, 3 + 2j])

np.isfinite(A)    #  Все ли элементы в A числа?

array([False, False, False,  True,  True,  True,  True,  True,  True])

np.isinf(A)    #  Есть ли в A бесконечности?

array([ True,  True, False, False, False, False, False, False, False])

np.isnan(A)    #  Есть ли в A значения nan?

array([False, False,  True, False, False, False, False, False, False])

np.iscomplex(A)    #  есть ли в A комплексные числа?

array([False, False, False, False, False, False, False, False,  True])

np.isreal(A)    #  Есть ли в A вещественные числа?

array([ True,  True,  True,  True,  True,  True,  True,  True, False])

Привычные нам логические операции выполняются над массивами булевых значений (массивы из значений True и False):

X = np.array([True, False, True, False])
Y = np.array([True, True, False, False])

np.logical_and(X, Y)    #  логическое "И"

array([ True, False, False, False])

np.logical_or(X, Y)    #  логическое "ИЛИ"

array([ True,  True,  True, False])

np.logical_not(X)    #  логическое "НЕ"

array([False,  True, False,  True])

np.logical_xor(X, Y)    # исключающее "ИЛИ"

array([False,  True,  True, False])

Помимо всего прочего, NumPy позволяет производить различные сравнения:

np.allclose([1, 2, 3], [1, 2, 2.99999])    #  Являются ли значения массивов близкими?

True

x = np.random.randint(4, size=5)
x

array([3, 1, 2, 1, 0])

y = np.random.randint(4, size=5)
y

array([1, 2, 3, 1, 0])

np.greater(x, y)    #  поэлементное x > y

array([ True, False, False, False, False])

np.less(x, y)    #  поэлементное x < y

array([False,  True,  True, False, False])

np.equal(x, y)    #  поэлементное x == y

array([False, False, False,  True,  True])