Никита Волков
Пакет numpy предоставляет $n$-мерные однородные массивы (все элементы одного типа); в них нельзя вставить или удалить элемент в произвольном месте. В numpy реализовано много операций над массивами в целом. Если задачу можно решить, произведя некоторую последовательность операций над массивами, то это будет столь же эффективно, как в C или matlab - львиная доля времени тратится в библиотечных функциях, написанных на C.
In [1]:
import numpy as np
Можно преобразовать список в массив.
In [2]:
a = np.array([0, 2, 1])
a, type(a)
Out[2]:
print печатает массивы в удобной форме.
In [3]:
print(a)
Класс ndarray имеет много методов.
In [4]:
set(dir(a)) - set(dir(object))
Out[4]:
Наш массив одномерный.
In [5]:
a.ndim
Out[5]:
В $n$-мерном случае возвращается кортеж размеров по каждой координате.
In [6]:
a.shape
Out[6]:
size - это полное число элементов в массиве; len - размер по первой координате (в 1-мерном случае это то же самое).
In [7]:
len(a), a.size
Out[7]:
numpy предоставляет несколько типов для целых (int16, int32, int64) и чисел с плавающей точкой (float32, float64).
In [8]:
a.dtype, a.dtype.name, a.itemsize
Out[8]:
Индексировать массив можно обычным образом.
In [9]:
a[1]
Out[9]:
Массивы - изменяемые объекты.
In [10]:
a[1] = 3
print(a)
Массивы, разумеется, можно использовать в for циклах. Но при этом теряется главное преимущество numpy - быстродействие. Всегда, когда это возможно, лучше использовать операции над массивами как едиными целыми.
In [11]:
for i in a:
print(i)
Массив чисел с плавающей точкой.
In [12]:
b = np.array([0., 2, 1])
b.dtype
Out[12]:
Точно такой же массив.
In [13]:
c = np.array([0, 2, 1], dtype=np.float64)
print(c)
Преобразование данных
In [14]:
print(c.dtype)
print(c.astype(int))
print(c.astype(str))
Массив, значения которого вычисляются функцией. Функции передаётся массив. Так что в ней можно использовать только такие операции, которые применимы к массивам.
In [15]:
def f(i):
print(i)
return i ** 2
a = np.fromfunction(f, (5,), dtype=np.int64)
print(a)
In [16]:
a = np.fromfunction(f, (5,), dtype=np.float64)
print(a)
Массивы, заполненные нулями или единицами. Часто лучше сначала создать такой массив, а потом присваивать значения его элементам.
In [17]:
a = np.zeros(3)
print(a)
In [18]:
b = np.ones(3, dtype=np.int64)
print(b)
Если нужно создать массив, заполненный нулями, длины другого массива, то можно использовать конструкцию
In [19]:
np.zeros_like(b)
Out[19]:
Функция arange подобна range. Аргументы могут быть с плавающей точкой. Следует избегать ситуаций, когда $(конец-начало)/шаг$ - целое число, потому что в этом случае включение последнего элемента зависит от ошибок округления. Лучше, чтобы конец диапазона был где-то посредине шага.
In [20]:
a = np.arange(0, 9, 2)
print(a)
In [21]:
b = np.arange(0., 9, 2)
print(b)
Последовательности чисел с постоянным шагом можно также создавать функцией linspace. Начало и конец диапазона включаются; последний аргумент - число точек.
In [22]:
a = np.linspace(0, 8, 5)
print(a)
Последовательность чисел с постоянным шагом по логарифмической шкале от $10^0$ до $10^1$.
In [23]:
b = np.logspace(0, 1, 5)
print(b)
Массив случайных чисел.
In [24]:
print(np.random.random(5))
Случайные числа с нормальным (гауссовым) распределением (среднее 0, среднеквадратичное отклонение 1).
In [25]:
print(np.random.normal(size=5))
In [26]:
print(a + b)
In [27]:
print(a - b)
In [28]:
print(a * b)
In [29]:
print(a / b)
In [30]:
print(a ** 2)
Когда операнды разных типов, они пиводятся к большему типу.
In [31]:
i = np.ones(5, dtype=np.int64)
print(a + i)
numpy содержит элементарные функции, которые тоже применяются к массивам поэлементно. Они называются универсальными функциями (ufunc).
In [32]:
np.sin, type(np.sin)
Out[32]:
In [33]:
print(np.sin(a))
Один из операндов может быть скаляром, а не массивом.
In [34]:
print(a + 1)
In [35]:
print(2 * a)
Сравнения дают булевы массивы.
In [36]:
print(a > b)
In [37]:
print(a == b)
In [38]:
c = a > 5
print(c)
Кванторы "существует" и "для всех".
In [39]:
np.any(c), np.all(c)
Out[39]:
Модификация на месте.
In [40]:
a += 1
print(a)
In [41]:
b *= 2
print(b)
In [42]:
b /= a
print(b)
При выполнении операций над массивами деление на 0 не возбуждает исключения, а даёт значения np.nan или np.inf.
In [43]:
print(np.array([0.0, 0.0, 1.0, -1.0]) / np.array([1.0, 0.0, 0.0, 0.0]))
In [44]:
np.nan + 1, np.inf + 1, np.inf * 0, 1. / np.inf
Out[44]:
Сумма и произведение всех элементов массива; максимальный и минимальный элемент; среднее и среднеквадратичное отклонение.
In [45]:
b.sum(), b.prod(), b.max(), b.min(), b.mean(), b.std()
Out[45]:
In [46]:
x = np.random.normal(size=1000)
x.mean(), x.std()
Out[46]:
Имеются встроенные функции
In [47]:
print(np.sqrt(b))
print(np.exp(b))
print(np.log(b))
print(np.sin(b))
print(np.e, np.pi)
Иногда бывает нужно использовать частичные (кумулятивные) суммы. В нашем курсе такое пригодится.
In [48]:
print(b.cumsum())
Функция sort возвращает отсортированную копию, метод sort сортирует на месте.
In [49]:
print(np.sort(b))
print(b)
In [50]:
b.sort()
print(b)
Объединение массивов.
In [51]:
a = np.hstack((a, b))
print(a)
Расщепление массива в позициях 3 и 6.
In [52]:
np.hsplit(a, [3, 6])
Out[52]:
Функции delete, insert и append не меняют массив на месте, а возвращают новый массив, в котором удалены, вставлены в середину или добавлены в конец какие-то элементы.
In [53]:
a = np.delete(a, [5, 7])
print(a)
In [54]:
a = np.insert(a, 2, [0, 0])
print(a)
In [55]:
a = np.append(a, [1, 2, 3])
print(a)
Есть несколько способов индексации массива. Вот обычный индекс.
In [56]:
a = np.linspace(0, 1, 11)
print(a)
In [57]:
b = a[2]
print(b)
Диапазон индексов. Создаётся новый заголовок массива, указывающий на те же данные. Изменения, сделанные через такой массив, видны и в исходном массиве.
In [58]:
b = a[2:6]
print(b)
In [59]:
b[0] = -0.2
print(b)
In [60]:
print(a)
Диапазон с шагом 2.
In [61]:
b = a[1:10:2]
print(b)
In [62]:
b[0] = -0.1
print(a)
Массив в обратном порядке.
In [63]:
b = a[len(a):0:-1]
print(b)
Подмассиву можно присвоить значение - массив правильного размера или скаляр.
In [64]:
a[1:10:3] = 0
print(a)
Тут опять создаётся только новый заголовок, указывающий на те же данные.
In [65]:
b = a[:]
b[1] = 0.1
print(a)
Чтобы скопировать и данные массива, нужно использовать метод copy.
In [66]:
b = a.copy()
b[2] = 0
print(b)
print(a)
Можно задать список индексов.
In [67]:
print(a[[2, 3, 5]])
Можно задать булев массив той же величины.
In [68]:
b = a > 0
print(b)
In [69]:
print(a[b])
In [70]:
a = np.array([[0.0, 1.0], [-1.0, 0.0]])
print(a)
In [71]:
a.ndim
Out[71]:
In [72]:
a.shape
Out[72]:
In [73]:
len(a), a.size
Out[73]:
In [74]:
a[1, 0]
Out[74]:
Атрибуту shape можно присвоить новое значение - кортеж размеров по всем координатам. Получится новый заголовок массива; его данные не изменятся.
In [75]:
b = np.linspace(0, 3, 4)
print(b)
In [76]:
b.shape
Out[76]:
In [77]:
b.shape = 2, 2
print(b)
Можно растянуть в одномерный массив
In [78]:
print(b.ravel())
Арифметические операции поэлементные
In [79]:
print(a + 1)
print(a * 2)
print(a + [0, 1]) # второе слагаемое дополняется до матрицы копированием строк
print(a + np.array([[0, 2]]).T) # .T - транспонирование
print(a + b)
Поэлементное и матричное (только в Python 3.5) умножение.
In [80]:
print(a * b)
In [81]:
print(a @ b)
In [82]:
print(b @ a)
Умножение матрицы на вектор.
In [83]:
v = np.array([1, -1], dtype=np.float64)
print(b @ v)
In [84]:
print(v @ b)
Если у вас Питон более ранней версии, то для работы с матрицами можно использовать класс np.matrix, в котором операция умножения реализуется как матричное умножение.
In [85]:
np.matrix(a) * np.matrix(b)
Out[85]:
Внешнее произведение $a_{ij}=u_i v_j$
In [86]:
u = np.linspace(1, 2, 2)
v = np.linspace(2, 4, 3)
print(u)
print(v)
In [87]:
a = np.outer(u, v)
print(a)
Двумерные массивы, зависящие только от одного индекса: $x_{ij}=u_j$, $y_{ij}=v_i$
In [88]:
x, y = np.meshgrid(u, v)
print(x)
print(y)
Единичная матрица.
In [89]:
I = np.eye(4)
print(I)
Метод reshape делает то же самое, что присваивание атрибуту shape.
In [90]:
print(I.reshape(16))
In [91]:
print(I.reshape(2, 8))
Строка.
In [92]:
print(I[1])
Цикл по строкам.
In [93]:
for row in I:
print(row)
Столбец.
In [94]:
print(I[:, 2])
Подматрица.
In [95]:
print(I[0:2, 1:3])
Можно построить двумерный массив из функции.
In [96]:
def f(i, j):
print(i)
print(j)
return 10 * i + j
print(np.fromfunction(f, (4, 4), dtype=np.int64))
Транспонированная матрица.
In [97]:
print(b.T)
Соединение матриц по горизонтали и по вертикали.
In [98]:
a = np.array([[0, 1], [2, 3]])
b = np.array([[4, 5, 6], [7, 8, 9]])
c = np.array([[4, 5], [6, 7], [8, 9]])
print(a)
print(b)
print(c)
In [99]:
print(np.hstack((a, b)))
In [100]:
print(np.vstack((a, c)))
Сумма всех элементов; суммы столбцов; суммы строк.
In [101]:
print(b.sum())
print(b.sum(axis=0))
print(b.sum(axis=1))
Аналогично работают prod, max, min и т.д.
In [102]:
print(b.max())
print(b.max(axis=0))
print(b.min(axis=1))
След - сумма диагональных элементов.
In [103]:
np.trace(a)
Out[103]:
In [104]:
X = np.arange(24).reshape(2, 3, 4)
print(X)
Суммирование (аналогично остальные операции)
In [105]:
# суммируем только по нулевой оси, то есть для фиксированных j и k суммируем только элементы с индексами (*, j, k)
print(X.sum(axis=0))
# суммируем сразу по двум осям, то есть для фиксированной i суммируем только элементы с индексами (i, *, *)
print(X.sum(axis=(1, 2)))
In [106]:
np.linalg.det(a)
Out[106]:
Обратная матрица.
In [107]:
a1 = np.linalg.inv(a)
print(a1)
In [108]:
print(a @ a1)
print(a1 @ a)
Решение линейной системы $au=v$.
In [109]:
v = np.array([0, 1], dtype=np.float64)
print(a1 @ v)
In [110]:
u = np.linalg.solve(a, v)
print(u)
Проверим.
In [111]:
print(a @ u - v)
Собственные значения и собственные векторы: $a u_i = \lambda_i u_i$. l - одномерный массив собственных значений $\lambda_i$, столбцы матрицы $u$ - собственные векторы $u_i$.
In [112]:
l, u = np.linalg.eig(a)
print(l)
In [113]:
print(u)
Проверим.
In [114]:
for i in range(2):
print(a @ u[:, i] - l[i] * u[:, i])
Функция diag от одномерного массива строит диагональную матрицу; от квадратной матрицы - возвращает одномерный массив её диагональных элементов.
In [115]:
L = np.diag(l)
print(L)
print(np.diag(L))
Все уравнения $a u_i = \lambda_i u_i$ можно собрать в одно матричное уравнение $a u = u \Lambda$, где $\Lambda$ - диагональная матрица с собственными значениями $\lambda_i$ по диагонали.
In [116]:
print(a @ u - u @ L)
Поэтому $u^{-1} a u = \Lambda$.
In [117]:
print(np.linalg.inv(u) @ a @ u)
Найдём теперь левые собственные векторы $v_i a = \lambda_i v_i$ (собственные значения $\lambda_i$ те же самые).
In [118]:
l, v = np.linalg.eig(a.T)
print(l)
print(v)
Собственные векторы нормированы на 1.
In [119]:
print(u.T @ u)
print(v.T @ v)
Левые и правые собственные векторы, соответствующие разным собственным значениям, ортогональны, потому что $v_i a u_j = \lambda_i v_i u_j = \lambda_j v_i u_j$.
In [120]:
print(v.T @ u)
In [121]:
from scipy.integrate import quad, odeint
from scipy.special import erf
In [122]:
def f(x):
return np.exp(-x ** 2)
Адаптивное численное интегрирование (может быть до бесконечности). err - оценка ошибки.
In [123]:
res, err = quad(f, 0, np.inf)
print(np.sqrt(np.pi) / 2, res, err)
In [124]:
res, err = quad(f, 0, 1)
print(np.sqrt(np.pi) / 2 * erf(1), res, err)
In [125]:
x = np.arange(0, 25, 0.5).reshape((5, 10))
# Сохраняем в файл example.txt данные x в формате с двумя точками после запятой и разделителем ';'
np.savetxt('example.txt', x, fmt='%.2f', delimiter=';')
Получится такой файл
In [126]:
! cat example.txt
Теперь его можно прочитать
In [127]:
x = np.loadtxt('example.txt', delimiter=';')
print(x)
Нам пригодится только модуль scipy.stats.
Полное описание http://docs.scipy.org/doc/scipy/reference/stats.html
In [128]:
import scipy.stats as sps
Общий принцип:
$X$ — некоторое распределение с параметрами params
Кроме того для непрерывных распределений определены функции
А для дискретных
Параметры могут быть следующими:
Для примера сгенерируем выборку размера $N = 200$ из распределения $\mathscr{N}(1, 9)$ и посчитаем некоторые статистики.
В терминах выше описанных функций у нас $X$ = sps.norm, а params = (loc=1, scale=3).
In [129]:
sample = sps.norm.rvs(size=200, loc=1, scale=3)
print('Первые 10 значений выборки:\n', sample[:10])
print('Выборочное среденее: %.3f' % sample.mean())
print('Выборочная дисперсия: %.3f' % sample.var())
In [130]:
print('Плотность:\t\t', sps.norm.pdf([-1, 0, 1, 2, 3], loc=1, scale=3))
print('Функция распределения:\t', sps.norm.cdf([-1, 0, 1, 2, 3], loc=1, scale=3))
In [131]:
print('Квантили:', sps.norm.ppf([0.05, 0.1, 0.5, 0.9, 0.95], loc=1, scale=3))
Cгенерируем выборку размера $N = 200$ из распределения $Bin(10, 0.6)$ и посчитаем некоторые статистики.
В терминах выше описанных функций у нас $X$ = sps.binom, а params = (n=10, p=0.6).
In [132]:
sample = sps.binom.rvs(size=200, n=10, p=0.6)
print('Первые 10 значений выборки:\n', sample[:10])
print('Выборочное среденее: %.3f' % sample.mean())
print('Выборочная дисперсия: %.3f' % sample.var())
In [133]:
print('Дискретная плотность:\t', sps.binom.pmf([-1, 0, 5, 5.5, 10], n=10, p=0.6))
print('Функция распределения:\t', sps.binom.cdf([-1, 0, 5, 5.5, 10], n=10, p=0.6))
In [134]:
print('Квантили:', sps.binom.ppf([0.05, 0.1, 0.5, 0.9, 0.95], n=10, p=0.6))
Отдельно есть класс для многомерного нормального распределения. Для примера сгенерируем выборку размера $N=200$ из распределения $\mathscr{N} \left( \begin{pmatrix} 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \right)$.
In [135]:
sample = sps.multivariate_normal.rvs(mean=[1, 1], cov=[[2, 1], [1, 2]], size=200)
print('Первые 10 значений выборки:\n', sample[:10])
print('Выборочное среденее:', sample.mean(axis=0))
print('Выборочная матрица ковариаций:\n', np.cov(sample.T))
Некоторая хитрость :)
In [136]:
sample = sps.norm.rvs(size=10, loc=np.arange(10), scale=0.1)
print(sample)
Бывает так, что надо сгенерировать выборку из распределения, которого нет в `scipy.stats`.
Для этого надо создать класс, который будет наследоваться от класса rv_continuous для непрерывных случайных величин и от класса rv_discrete для дискретных случайных величин.
Пример есть на странице http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.html#scipy.stats.rv_continuous
Для примера сгенерируем выборку из распределения с плотностью $f(x) = \frac{4}{15} x^3 I\{x \in [1, 2] = [a, b]\}$.
In [137]:
class cubic_gen(sps.rv_continuous):
def _pdf(self, x):
return 4 * x ** 3 / 15
cubic = cubic_gen(a=1, b=2, name='cubic')
sample = cubic.rvs(size=200)
print('Первые 10 значений выборки:\n', sample[:10])
print('Выборочное среденее: %.3f' % sample.mean())
print('Выборочная дисперсия: %.3f' % sample.var())
Если дискретная случайная величина может принимать небольшое число значений, то можно не создавать новый класс, как показано выше, а явно указать эти значения и из вероятности.
In [138]:
some_distribution = sps.rv_discrete(name='some_distribution', values=([1, 2, 3], [0.6, 0.1, 0.3]))
sample = some_distribution.rvs(size=200)
print('Первые 10 значений выборки:\n', sample[:10])
print('Выборочное среденее: %.3f' % sample.mean())
print('Частота значений по выборке:', (sample == 1).mean(), (sample == 2).mean(), (sample == 3).mean())