7. Математика в NumPy

Когда дело доходит до математики и вычислений над какими-то наборами данных всегда появляется один итот же вопрос - что дороже, время потраченное на программную реализацию алгоритма вычислений или время компьютера, которое необходимо компьютеру для выполнения данной программы. Ответ очень прост, если за стеной вашего кабинета нет суперкомпьютера, то ваше время несомненно дороже. Да, знание ассемблера, C или Fortran позволит вам создавать очень быстрые программы, но на создание потребуется гораздо больше времени, чем в случае создания тех же программ на языке Python.

NumPy - это удачная попытка соединить скорость вычислений на языках C и Fortran со скоростью разработки программ на языке Python. Еще одним огромным достоинством NumPy является большой набор уже готовых математических функций, охватывающий многие разделы математики. Причём, большинство данных функций могут быть легко оптимизированы под конкретную задачу. А механизм транслирования NumPy позволяет избежать написания большого количества циклов, особенно вложенных, которые так часто встречаются при работе с многомерными массивами.

7.1. Базовые математические операции

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

>>> 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 попытается выполнить конкатенацию, а не
... #  сложение.
... 
>>> a + c
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: can only concatenate list (not "int") to list
>>> 
>>> 
>>> #  Теперь умножим массив 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]], dtype=int32)
>>> 
>>> a%3    #  остаток от деления
array([[2, 1],
       [2, 1]], dtype=int32)
>>> 
>>> a**3    #  возведение в степень
array([[ 125,  343],
       [1331, 2197]], dtype=int32)
>>>
>>> 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)    #  натуральный логарифм элементов массива
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
Traceback (most recent call last):
  File "", line 1, in 
TypeError: Cannot cast ufunc add output from 
dtype('float64') to dtype('int32') with casting rule 'same_kind'
>>>
>>>
>>> #  А вот преобразование целочисленного типа в вещественный возможно
>>> b += a
>>> b
array([  6.1,   8.2,  10.3,  12.4,  14.5])
>>>
>>> #  Операции не должны изменять тип массива
...
>>> a /= 2    #  Требуется прелбразование типа данных
Traceback (most recent call last):
  File "", line 1, in 
TypeError: ufunc 'true_divide' output (typecode 'd') could not be 
coerced to provided output parameter (typecode 'l') according to
the casting rule ''same_kind''
>>>
>>> a
array([ 6,  8, 10, 12, 14])
>>>
>>> a **=2    #  Преобразование типа данных не требуется
>>> a
array([ 36,  64, 100, 144, 196])
>>> 
>>> a //= 5
>>> a
array([ 7, 12, 20, 28, 39])
>>>
>>> a %= 5
>>> a
array([2, 2, 0, 3, 4])

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

>>> 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('int32')
>>> b.dtype
dtype('float64')
>>>
>>> c = a + b
>>> c
array([ 0.  ,  2.25,  4.5 ,  6.75,  9.  ])
>>> c.dtype
dtype('float64')
>>>
>>> d = np.array([1 +1j, 2 + 2j, 3 + 3j, 4 + 4j, 5 + 5j])
>>> d
array([ 1.+1.j,  2.+2.j,  3.+3.j,  4.+4.j,  5.+5.j])
>>> d.dtype
dtype('complex128')
>>>
>>> e = c*d
>>> e
array([  0.0 +0.j ,   4.5 +4.5j,  13.5+13.5j,  27.0+27.j ,  45.0+45.j ])
>>> e.dtype
dtype('complex128')

Применение логических операций к массивам, так же возможно и так же выполняется поэлементно. Результатом таких операций является массив булевых значений (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], dtype=int32)
>>> 
>>> 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]], dtype=int32)

Хотя, если честно, их размеры должны быть не равны, а должны быть совместимыми. Если их размеры совместимы, т.е. один массив может быть растянут до размеров другого, то в дело включается механизм транслирования массивов 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])
>>> 
>>> 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])

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

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

>>> np.log(0)
__main__:1: RuntimeWarning: divide by zero encountered in log
-inf

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

>>> a = np.array([0])
>>> 1/a
__main__:1: RuntimeWarning: divide by zero encountered in true_divide
array([inf])

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

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

>>> a = np.array([0, 1, np.inf])
>>> 
>>> np.cos(a)
__main__:1: RuntimeWarning: invalid value encountered in cos
array([1.        , 0.54030231,        nan])

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

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

>>> np.cos(0)/np.sin(0)
__main__:1: RuntimeWarning: divide by zero encountered in double_scalars
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
>>> #  Более точное значение числа pi
... #  3,141592653589793238462643383279

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

Тем не менее и на этом сюрпризы не заканчиваются. Если число 1.633123935319537e+16 является самым больши, которое может появиться при вычислениях, оно вполне ожидаемо должно появиться в самых разных ситуациях. Например:

>>> a = np.array([2, 2*10, 2**1000])
>>> a
array([2, 20,
       107150860718626732094842504906000181056140481170553360
       744375038837035105112493612249319837881569585812759467
       291755314682518714528569231404359845775746985748039345
       677748242309854210746050623711418779541821530464749835
       819412673987675591655439460770629145711964776865421676
       60429831652624386837205668069376],
      dtype=object)
>>> 
>>> 1/a
array([0.5, 0.05, 9.332636185032189e-302], dtype=object)
>>> 
>>> a%7
array([2, 6, 2], dtype=object)

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

>>> a = np.array([2])
>>> 
>>> a**(10**100)    # нажав enter, не пугайтесь, ваш компьютер зависнет

В заключение могу лишь сказать, что все предельные случаи требуют кардинальных решений. Некоторые решения имеются в самом NumPy, некоторые предоставляют другие пакеты. Если вам необходимы точные решения, то лучше обратиться к системам компьютерной алгебры и символьных вычислений, например пакету SymPy - маленький, но мощьный пакет Python для символьных вычислений. Если вы решили отправиться в самые дебри теории чисел, алгебры и криптографии, то лучшим решением окажется программа GAP. Программа GAP не является программой Python, но имеет Python интерфейс в замечательной программе Sage, которая определенно заслуживает вашего внмания.


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

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

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

>>> 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,5)
array([[1069, 1558],
       [2337, 3406]])
>>>
>>> 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.9999999999999982
>>>
>>> 1*3 - 3*4    #  Результат должен быть целым числом
-9

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

>>> np.linalg.det(a)
-8.9999999999999982
>>> round(np.linalg.det(a))
-9.0
>>>
>>> b = np.arange(1,48,3).reshape(4,4)
>>> np.linalg.det(b)
-1.0223637331664275e-27
>>> round(np.linalg.det(b))
-0.0

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

>>> a
matrix([[1, 3],
        [4, 3]])
>>>
>>> a.T
matrix([[1, 4],
        [3, 3]])
>>>
>>> b
array([[ 1,  4,  7, 10],
       [13, 16, 19, 22],
       [25, 28, 31, 34],
       [37, 40, 43, 46]])
>>>
>>> b.T
array([[ 1, 13, 25, 37],
       [ 4, 16, 28, 40],
       [ 7, 19, 31, 43],
       [10, 22, 34, 46]])

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

>>> a = np.array([[1,3],[4,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.])

7.3. Статистика

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

Элементарные статистические функции:

>>> a = np.random.randint(20, size = (5, 5))
>>> a
array([[ 9, 11,  7, 19, 10],
       [17,  4, 19, 14,  2],
       [ 1,  1,  6, 15,  9],
       [18,  6, 18, 18, 18],
       [ 2,  3, 19,  0,  9]])
>>> 
>>> np.amin(a)    #  Минимальный элемент массива
0
>>> np.amax(a)    #  максимальный элемент
19
>>> np.amin(a, axis = 0)    #  минимальный элементо вдоль первой оси (столбцы)
array([1, 1, 6, 0, 2])
>>> np.amin(a, axis = 1)    #  минимальный элементо вдоль второй оси (строки)
array([7, 2, 1, 6, 0])
>>> 
>>> 
>>> #  Процентили:
... 
>>> np.percentile(a, 25)
4.0
>>> np.percentile(a, 50)
9.0
>>> np.percentile(a, 75)
18.0

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

>>> a = np.random.randint(13, size = (5, 5))
>>> a
array([[ 3,  5,  0,  1,  5],
       [12,  0,  9,  5, 12],
       [ 0,  1,  1,  4,  4],
       [ 6,  5,  5,  3,  1],
       [ 4,  7,  4,  7,  8]])
>>> 
>>> np.median(a)    #  медиана элементов массива
4.0
>>> np.mean(a)    #  среднее арифметическое
4.48
>>> np.var(a)    #  дисперсия
11.0496
>>> np.std(a)    #  стандартное отклонение
3.324093861490677

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

>>> 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, ..., 21, 20, 23],
       [ 4,  1,  6, ..., 19, 15, 22]])
>>> 
>>> np.corrcoef(XY)
array([[1.        , 0.93158099],
       [0.93158099, 1.        ]])
>>> 
>>> #  Линейный коэфициент корреляции Пирсона
... #  величин 'x' и 'z'
... 
>>> XZ = np.stack((x, z), axis = 0)
>>> XZ
array([[ 1,  4,  3, ..., 21, 20, 23],
       [29, 22, 24, ..., 11,  9, 10]])
>>> 
>>> np.corrcoef(XZ)
array([[ 1.        , -0.92342264],
       [-0.92342264,  1.        ]])
>>> 
>>> 
>>> #  Кросс-корреляции:
...
>>> np.correlate(x, y)
array([1736])
>>> np.correlate(x, z)
array([1486])
>>> np.correlate(y, z)
array([1670])
>>>
>>>
>>> #  Ковариационные матрицы:
... 
>>> np.cov(XY)
array([[63.65555556, 49.82222222],
       [49.82222222, 44.93333333]])
>>> np.cov(XZ)
array([[ 63.65555556, -48.25555556],
       [-48.25555556,  42.9       ]])
>>> 
>>> np.cov(x)
array(63.65555556)
>>> np.cov(y)
array(44.93333333)
>>> np.cov(z)
array(42.9)

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


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

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

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

>>> np.random.rand()    #  случайное число от 0 до 1
0.9129353336551852
>>>
>>> np.random.rand(10)    #  одномерный массив случайных значений
array([0.06480376, 0.62333573, 0.0363841 , 0.23728226, 0.58031617,
       0.64754324, 0.87082152, 0.56534312, 0.40102058, 0.79272714])
>>> 
>>> np.random.rand(3, 4)    #  двумерный массив случайных значений
array([[0.01464067, 0.92706517, 0.61968748, 0.57490629],
       [0.02952706, 0.07091393, 0.66877515, 0.20989394],
       [0.94344177, 0.71377153, 0.49056813, 0.99597689]])
>>>
>>> np.random.randn(10)    #  случайные значения с нормальным распределением
array([ 1.24693967,  0.79698643, -0.72896891,  0.28832121, -1.20637466,
        0.26280556, -0.30526188, -0.13138396,  1.17051473, -1.55996539])
>>>
>>> np.random.randint(10)    #  случайное целое число от 0 до 10
9
>>>
>>> np.random.randint(10, 100)    #  случайное целое число от 10 до 100
80
>>>
>>> np.random.randint(10, size = 7)    #  одномерный массив случайных целых чисел
array([0, 8, 1, 0, 0, 7, 3])
>>> 
>>> np.random.randint(10, size = (4, 4))    #  двумерный массив случайных целых чисел
array([[2, 1, 1, 3],
       [1, 5, 9, 6],
       [6, 7, 7, 6],
       [0, 9, 8, 3]])

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

>>> x = np.arange(7)
>>> x
array([0, 1, 2, 3, 4, 5, 6])
>>> 
>>> np.random.shuffle(x)    #  перетасовывает содержимое массива
>>> x
array([2, 5, 3, 4, 6, 1, 0])
>>> np.random.shuffle(x)
>>> x
array([4, 0, 2, 5, 3, 6, 1])
>>> 
>>> np.random.permutation(7)    #  выполняет тоже самое, не требуя входного массива
array([6, 0, 4, 3, 5, 1, 2])
>>> 
>>> y = np.arange(9).reshape(3, 3)
>>> y
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])
>>> 
>>> np.random.shuffle(y)    #  перестановки выполняются только по 1-й оси
>>> y
array([[0, 1, 2],
       [6, 7, 8],
       [3, 4, 5]])

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

>>> np.random.beta(0.1, 0.6, size = 5)    #  бета распределение
array([5.05340612e-05, 3.28413864e-01, 2.20400277e-13, 9.05029465e-23,
       1.63347756e-03])
>>> 
>>> np.random.gamma(shape = 0.8, scale = 1.7, size = 5)    # гамма распределение
array([1.42116185, 0.10924208, 0.1738796 , 2.74877293, 1.69352153])
>>> 
>>> np.random.pareto(3.5, size = 5)    #  Паретто распределение
array([0.07597227, 0.17833497, 0.03677746, 1.73382556, 0.70577262])
>>> 
>>> np.random.chisquare(2.2, size = 5)     #  хи-квадрат распределение
array([5.73476387, 9.01545569, 0.6511614 , 0.87320945, 2.99932556])

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

>>> np.random.get_state()    #  Вы можете узнать состояние генератора
('MT19937', array([2609792837, 1032626479, 4264830444, ..., 3281379045, 3584140162,
       3132090685], dtype=uint32), 362, 1, -1.598984180826123)
>>> 
>>> np.random.seed(123)    #  устанавливает состояние генератора
>>> np.random.rand(5)
array([0.69646919, 0.28613933, 0.22685145, 0.55131477, 0.71946897])
>>> 
>>> np.random.seed(123)    #  установим генератор в тоже состояние
>>> np.random.rand(5)     #  генератор выдаст теже случайные данные
array([0.69646919, 0.28613933, 0.22685145, 0.55131477, 0.71946897])
>>> 
>>> 
>>> np.random.rand(5)     #  после этого состояние генератора изменится
array([0.42310646, 0.9807642 , 0.68482974, 0.4809319 , 0.39211752])

7.5. Дискретное преобразование Фурье

Если данные в ваших массивах - это сигналы: звуки, изображения, радиоволны, котировки акций и т.д., то вам наверняка понадобится дискретное преобразование Фурье. В NumPy представлены методы быстрого дискретного преобразования Фурье для одномерных, двумерных и многомерных сигналов, а так же некоторые вспомогательные функции. Рассмотрим некоторые простые примеры.

Одномерное дискретное преобразование Фурье:

>>> import numpy as np
>>> x = np.arange(5)
>>> y = np.sin(x)
>>> y
array([ 0.        ,  0.84147098,  0.90929743,  0.14112001, -0.7568025 ])
>>> 
>>> np.fft.fft(y)
array([ 1.13508592+0.        j, -0.82364155-1.97157177j,
        0.25609859-0.20886144j,  0.25609859+0.20886144j,
       -0.82364155+1.97157177j])
>>> 
>>> np.fft.ifft(y)    #  одномерное обратное преобразование Фурье
array([ 0.22701718+0.        j, -0.16472831+0.39431435j,
        0.05121972+0.04177229j,  0.05121972-0.04177229j,
       -0.16472831-0.39431435j])

Двумерное дискретное преобразование Фурье:

>>> x = np.arange(4)
>>> y = np.arange(4)
>>> xy = np.meshgrid(x, y, sparse = True)
>>> 
>>> z = np.cos(xy[0]) + np.sin(xy[1])
>>> z
array([[ 1.        ,  0.54030231, -0.41614684, -0.9899925 ],
       [ 1.84147098,  1.38177329,  0.42532415, -0.14852151],
       [ 1.90929743,  1.44959973,  0.49315059, -0.08069507],
       [ 1.14112001,  0.68142231, -0.27502683, -0.84887249]])
>>> 
>>> np.fft.fft2(z)
array([[ 8.10420557e+00+0.00000000e+00j,  5.66458735e+00-6.12117921e+00j,
         4.13417342e+00+0.00000000e+00j,  5.66458735e+00+6.12117921e+00j],
       [-3.63718971e+00-2.80140391e+00j,  0.00000000e+00+0.00000000e+00j,
         0.00000000e+00-2.22044605e-16j,  0.00000000e+00+0.00000000e+00j],
       [-2.93174264e-01+0.00000000e+00j,  0.00000000e+00+0.00000000e+00j,
         2.22044605e-16+0.00000000e+00j,  0.00000000e+00+0.00000000e+00j],
       [-3.63718971e+00+2.80140391e+00j,  0.00000000e+00+0.00000000e+00j,
         0.00000000e+00+2.22044605e-16j,  0.00000000e+00+0.00000000e+00j]])
>>> 
>>> np.fft.ifft2(z)    #  двумерное обратное преобразование Фурье
array([[ 5.06512848e-01+0.00000000e+00j,  3.54036709e-01+3.82573701e-01j,
         2.58385839e-01+0.00000000e+00j,  3.54036709e-01-3.82573701e-01j],
       [-2.27324357e-01+1.75087744e-01j,  0.00000000e+00+0.00000000e+00j,
         0.00000000e+00+1.38777878e-17j,  0.00000000e+00+0.00000000e+00j],
       [-1.83233915e-02+0.00000000e+00j,  0.00000000e+00+0.00000000e+00j,
         1.38777878e-17+0.00000000e+00j,  0.00000000e+00+0.00000000e+00j],
       [-2.27324357e-01-1.75087744e-01j,  0.00000000e+00+0.00000000e+00j,
         0.00000000e+00-1.38777878e-17j,  0.00000000e+00+0.00000000e+00j]])

Очень часто при спектральном анализе используются оконные функции (оконное преобразование Фурье), некоторые из которых так же представлены в NumPy

>>> np.bartlett(10)    #  окно Барлетта (треугольное окно)
array([0.        , 0.22222222, 0.44444444, 0.66666667, 0.88888889,
       0.88888889, 0.66666667, 0.44444444, 0.22222222, 0.        ])
>>> 
>>> np.blackman(10)    #  окно Блэкмана
array([-1.38777878e-17,  5.08696327e-02,  2.58000502e-01,  6.30000000e-01,
        9.51129866e-01,  9.51129866e-01,  6.30000000e-01,  2.58000502e-01,
        5.08696327e-02, -1.38777878e-17])
>>> 
>>> np.hamming(10)    #  окно Хэмминга
array([0.08      , 0.18761956, 0.46012184, 0.77      , 0.97225861,
       0.97225861, 0.77      , 0.46012184, 0.18761956, 0.08      ])
>>> 
>>> np.hanning(10)    #  окно Ханна (Хеннинга)
array([0.        , 0.11697778, 0.41317591, 0.75      , 0.96984631,
       0.96984631, 0.75      , 0.41317591, 0.11697778, 0.        ])
>>> 
>>> np.kaiser(10, 4)     #  окно Кайзера
array([0.08848053, 0.29425961, 0.56437221, 0.82160913, 0.97885093,
       0.97885093, 0.82160913, 0.56437221, 0.29425961, 0.08848053])

7.6. Прочие математические разделы

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

7.6.1. Множества

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

>>> 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]])

Функция unique() так же позволяет получить индексы входного массива с первым вхождением уникальных элементов, а так же индексы уникального массива которые восстанавливают исходный массив:

>>> x = np.array([1, 1, 5, 1, 4, 5, 2, 4])
>>> 
>>> np.unique(x, return_index=True)    #  индексы первого вхождения
(array([1, 2, 4, 5]), array([0, 6, 4, 2], dtype=int32))
>>> x[0], x[6], x[4], x[2]
(1, 2, 4, 5)
>>> 
>>> np.unique(x, return_inverse=True)    #  индексы уникальных элементов в 'x'
(array([1, 2, 4, 5]), array([0, 0, 3, 0, 2, 3, 1, 2], dtype=int32))
>>>
>>> uniq, index = np.unique(x, return_inverse=True)
>>> uniq
array([1, 2, 4, 5])
>>> index
array([0, 0, 3, 0, 2, 3, 1, 2], dtype=int32)
>>> 
>>> uniq[index]    #  восстановление исходного массива
array([1, 1, 5, 1, 4, 5, 2, 4])

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

>>> 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])

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

Логические функции 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])
>>> 
>>> np.isfortran(A)    #  Является ли A фортран-смежным массивом?
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)
>>> y = np.random.randint(4, size = 5)
>>> x
array([0, 0, 3, 0, 3])
>>> y
array([1, 0, 2, 3, 2])
>>> 
>>> np.greater(x, y)    #  поэлементное x > y 
array([False, False,  True, False,  True])
>>> 
>>> np.less(x, y)    #  поэлементное x < y
array([ True, False, False,  True, False])
>>> 
>>> np.equal(x, y)    #  поэлементное x == y
array([False,  True, False, False, False])

7.6.3. Многочлены

Многочлены, довольно часто, встречаются в самых разных областях: физика, инженерное дело, теория приближений, криптография, теория вероятностей, комбинаторика и т.д. Многочлены определенного вида выделяются в особые классы, каждый из которых имеет свою нишу применения. В NumPy представлены функции для работы с классами многочленов Чебышёва, Лежандра, Лагерра и Эрмита. Помимо прочего, NumPy предоставляет широкий набор функций для работы с полиномами в общем виде, т.е. полиномы любого вида (Бернулли, Гегенбауэра, Шапиро, Якоби, Роджерса в том числе) могут быть легко сконструированы пользователем самостоятельно.

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

>>> from numpy.polynomial import Polynomial as P
>>> 
>>> p = P([3, 2, 1])    #  создание многочлена 3 + 2*x + x**2
>>> p
Polynomial([3., 2., 1.], domain=[-1,  1], window=[-1,  1])
>>> 
>>> p + p    #  сложение многочленов
Polynomial([6., 4., 2.], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> p - p    #  вычитание
Polynomial([0.], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> p*p    #  умножение
Polynomial([ 9., 12., 10.,  4.,  1.], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> p**2    #  возведение в степень
Polynomial([ 9., 12., 10.,  4.,  1.], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> p//P([1, 2])    #  частное полиномов 3 + 2*x + x**2 и 1 + 2*x
Polynomial([0.75, 0.5 ], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> p%P([1, 2])    #  остаток от деления 3 + 2*x + x**2 на 1 + 2*x
Polynomial([2.25], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> 
>>> #  Получение значений полинома при различных 
... #  значениях переменной:
... 
>>> x = np.arange(7)
>>> x
array([0, 1, 2, 3, 4, 5, 6])
>>> 
>>> p(x)
array([ 3.,  6., 11., 18., 27., 38., 51.])
>>> 
>>>
>>> #  Рассматривая полином, как математическую функцию
... #  мы можем применять ее к другим полиномам
... 
>>> p(p)
Polynomial([18., 16., 12.,  4.,  1.], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> p.roots()    #  корни многочлена
array([-1.-1.41421356j, -1.+1.41421356j])
>>> 
>>> 
>>> #  Полиномы могут быть проинтегрированны и продифференцированны:
... 
>>> p.integ()    #  интегрирование полинома
Polynomial([0.        , 3.        , 1.        , 0.33333333], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> p.deriv(1)    #  первая производня полинома
Polynomial([2., 2.], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> p.deriv(2)    #  вторая производня полинома
Polynomial([2.], domain=[-1.,  1.], window=[-1.,  1.])

Иногда бывает крайне удобно создавать полиномы не указывая их коэфициенты, а указывая их корни:

>>> p1 = P.fromroots([-1, 1])
>>> p1
Polynomial([-1.,  0.,  1.], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> p2 = P.fromroots([-3.14, 3.14])
>>> p2
Polynomial([-9.8596,  0.    ,  1.    ], domain=[-1.,  1.], window=[-1.,  1.])
>>> 
>>> p3 = P.fromroots([1 + 1j, 2 + 2j, 3 + 3j])
>>> p3
Polynomial([12.-12.j,  0.+22.j, -6. -6.j,  1. +0.j], domain=[-1.,  1.], window=[-1.,  1.])

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

>>> x = np.arange(10)
>>> 
>>> p = np.polynomial.Polynomial([3, 2, 1])    #  обычный полином
>>> p
Polynomial([3., 2., 1.], domain=[-1,  1], window=[-1,  1])
>>> p(x)
array([  3.,   6.,  11.,  18.,  27.,  38.,  51.,  66.,  83., 102.])
>>> 
>>> 
>>> Cheb = np.polynomial.Chebyshev([3, 2, 1])    #  полином Чебышёва
>>> Cheb
Chebyshev([3., 2., 1.], domain=[-1,  1], window=[-1,  1])
>>> Cheb(x)
array([  2.,   6.,  14.,  26.,  42.,  62.,  86., 114., 146., 182.])
>>> 
>>> 
>>> Leg = np.polynomial.Legendre([3, 2, 1])    #  полином Лежандра
>>> Leg
Legendre([3., 2., 1.], domain=[-1,  1], window=[-1,  1])
>>> Leg(x)
array([  2.5,   6. ,  12.5,  22. ,  34.5,  50. ,  68.5,  90. , 114.5,
       142. ])
>>> 
>>> 
>>> Lag = np.polynomial.Laguerre([3, 2, 1])    #  полином Лагерра
>>> Lag
Laguerre([3., 2., 1.], domain=[0, 1], window=[0, 1])
>>> Lag(x)
array([ 6. ,  2.5,  0. , -1.5, -2. , -1.5,  0. ,  2.5,  6. , 10.5])
>>> 
>>> 
>>> Herm = np.polynomial.Hermite([3, 2, 1])    #  полином Эрмита
>>> Herm
Hermite([3., 2., 1.], domain=[-1,  1], window=[-1,  1])
>>> Herm(x)
array([  1.,   9.,  25.,  49.,  81., 121., 169., 225., 289., 361.])