In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
В статистике часто приходится считать выборочное среднее, т.е. по данной выборке значений $x_k$, $k=1..N$, нужно вычислить $$\bar x=\frac1N\sum_{k=1}^N x_k.$$ С точки зрения математики не имеет значения, как считать указанную сумму, так как результат сложения всегда будет один и тот же. Однако при вычислениях с плавающей запятой ответ будет зависеть от порядка выполнения операций, хотя бы потому, что сложения чисел с плавающей запятой не ассоциативно. Но будет ли зависеть точность вычислений от порядка операций? Давайте это проверим.
Сконструируем выборку таким образом, что сумма всех элементов равна $1$, и порядок элементов меняется в широком диапазоне. Для этого разобьем единицу на $K$ частей, и $k$-ую часть разобьем на $7^k$ равных значений. Полученные элементы перемешаем.
In [2]:
base=10 # параметр, может принимать любые целые значения > 1
def exact_sum(K):
"""Точное значение суммы всех элементов."""
return 1.
def samples(K):
""""Элементы выборки"."""
# создаем K частей из base^k одинаковых значений
parts=[np.full((base**k,), float(base)**(-k)/K) for k in range(0, K)]
# создаем выборку объединяя части
samples=np.concatenate(parts)
# перемешиваем элементы выборки и возвращаем
return np.random.permutation(samples)
def direct_sum(x):
"""Последовательная сумма всех элементов вектора x"""
s=0.
for e in x:
s+=e
return s
In [3]:
def number_of_samples(K):
"""Число элементов в выборке"""
return np.sum([base**k for k in range(0, K)])
def exact_mean(K):
"""Значение среднего арифметического по выборке с близкой к машинной точностью."""
return 1./number_of_samples(K)
def exact_variance(K):
"""Значение оценки дисперсии с близкой к машинной точностью."""
# разные значения элементов выборки
values=np.asarray([float(base)**(-k)/K for k in range(0, K)], dtype=np.double)
# сколько раз значение встречается в выборке
count=np.asarray([base**k for k in range(0, K)])
return np.sum(count*(values-exact_mean(K))**2)/number_of_samples(K)
Создадим выборку из значений, отличающихся на 6 порядков, и просуммируем элементы выборки.
In [4]:
K=7 # число слагаемых
x=samples(K) # сохраняем выборку в массив
print("Число элементов:", len(x))
print("Самое маленькое и большое значения:", np.min(x), np.max(x))
exact_sum_for_x=exact_sum(K) # значение суммы с близкой к машинной погрешностью
direct_sum_for_x=direct_sum(x) # сумма всех элементов по порядку
def relative_error(x0, x):
"""Погрешность x при точном значении x0"""
return np.abs(x0-x)/np.abs(x)
print("Погрешность прямого суммирования:", relative_error(exact_sum_for_x, direct_sum_for_x))
Попробуем теперь просуммировать элементы в порядке возрастания.
In [5]:
sorted_x=x[np.argsort(x)]
sorted_sum_for_x=direct_sum(sorted_x)
print("Погрешность суммирования по возрастанию:", relative_error(exact_sum_for_x, sorted_sum_for_x))
Попробуем просуммировать в порядке убывания.
In [6]:
sorted_x=x[np.argsort(x)[::-1]]
sorted_sum_for_x=direct_sum(sorted_x)
print("Погрешность суммирования по убыванию:", relative_error(exact_sum_for_x, sorted_sum_for_x))
Таким образом погрешность результата зависит от порядка суммирования. Как можно объяснить этот эффект?
На практике суммирование предпочтительно проводить не наивным способом, а компенсационным суммированием (см. алгоритм Кэхэна.
In [7]:
def Kahan_sum(x):
s=0.0 # частичная сумма
c=0.0 # сумма погрешностей
for i in x:
y=i-c # первоначально y равно следующему элементу последовательности
t=s+y # сумма s может быть велика, поэтому младшие биты y будут потеряны
c=(t-s)-y # (t-s) отбрасывает старшие биты, вычитание y восстанавливает младшие биты
s=t # новое значение старших битов суммы
return s
Kahan_sum_for_x=Kahan_sum(x) # сумма всех элементов по порядку
print("Погрешность суммирования по Кэхэну:", relative_error(exact_sum_for_x, Kahan_sum_for_x))
Сумма первых $N$ элементов последовательности из задания 4 может быть найдена явна: $$\sum_{k=1}^N\sin k=\frac{1}{2}\bigg(\sin n-\mathrm{ctg}\frac{1}{2}\cos n+\mathrm{ctg}\frac{1}{2}\bigg).$$
Кроме вычисления оценки математического ожидания, часто требуется вычислить оценку среднеквадратического отклонения или его квадрата - дисперсии. Дисперсия $D[X]$ случайной величины $X$ определена через математическое ожидание $E[X]$ следующим образом: $$D[X]=E[(X-E[X])^2].$$ Для оценки дисперсии мы можем воспользоваться формулой для оценки математического ожидания через выборочное среднее: $$E[X]\approx\frac1N\sum_{n=1}^N x_n,$$ т.е. можно предложить следующую формулу для оценки дисперсии (первая формула): $$D[X]\approx\frac1N\sum_{n=1}^N\left(x_n-\frac1N\sum_{n=1}^Nx_n\right)^2.$$ Полученная оценка является смещенной, т.е. ее мат. ожидание не совпадает с верным значением дисперсии, поэтому на практике нужно использовать следующую несмещенную оценку: $$D[X]\approx\frac1{N-1}\sum_{n=1}^N\left(x_n-\frac1N\sum_{n=1}^Nx_n\right)^2,$$ однако в этой работе мы удовлетворимся смещенной оценкой. К сожалению, наша формула не позволяет обновлять значения дисперсии при добавлении значения в выборку, так как требует двух проходов по выборке: сначала считается среднее, затем считается дисперсия. Однако в учебниках теории вероятности можно встретить и другую эквивалентную формулу для дисперсии, получим ее, опираясь на свойства мат. ожидания: $$D[X]=E[(X-E[X])^2]=E[X^2-2E[X]X+E[X]^2]=E[X^2]-2E[X]E[X]+E[E[X]^2]=E[X^2]-E[X]^2.$$ Снова заменяя мат. ожидание на выборочное среднее, получаем новую оценку для дисперсии (вторая формула): $$D[X]\approx \frac1N\sum_{n=1}^N x_n^2-\left(\frac1N\sum_{n=1}^Nx_n\right)^2.$$ Вторая формулы для вычисления дисперсии более привлекательна, так как обе суммы могут вычисляться одновременно, а значения мат. ожидания и дисперсии вычислить, последовательно добавляя значения. Действительно, введем обозначения для оценок мат. ожидания и дисперсии по первым $n$ членам выборки: $$E_n=\frac1n\sum_{k=1}^n x_n,\quad D_n=\frac1n\sum_{k=1}^n x_n^2-E_n^2.$$ Отсюда легко вывести рекуррентные формулы: $$E_{n}=\frac{x_{n}+(n-1)E_{n-1}}{n},\quad D_{n}=\frac{x_{n}^2+(n-1)D_{n-1}}{n}-E_{n}^2.$$ Хотя эти формулы и просты, погрешность вычислений по второй формуле может быть значительно выше, чем по первой. Проверим это.
Рассмотрим выборку, среднее для которой на порядки больше среднеквадратического отклонения. Пусть ровно половина значений больше среднего на $delta$, а половина меньше на $delta$. Оценка дисперсии и мат. ожидания в этом случае легко вычисляются явно.
In [8]:
# параметры выборки
mean=1e6 # среднее
delta=1e-5 # величина отклонения от среднего
def samples(N_over_two):
"""Генерирует выборку из 2*N_over_two значений с данным средним и среднеквадратическим
отклонением."""
x=np.full((2*N_over_two,), mean, dtype=np.double)
x[:N_over_two]+=delta
x[N_over_two:]-=delta
return np.random.permutation(x)
def exact_mean():
"""Значение среднего арифметического по выборке с близкой к машинной точностью."""
return mean
def exact_variance():
"""Значение оценки дисперсии с близкой к машинной точностью."""
return delta**2
x=samples(1000000)
In [9]:
print("Размер выборки:", len(x))
print("Среднее значение:", exact_mean())
print("Оценка дисперсии:", exact_variance())
print("Ошибка среднего для встроенной функции:",relative_error(exact_mean(),np.mean(x)))
print("Ошибка дисперсии для встроенной функции:",relative_error(exact_variance(),np.var(x)))
In [10]:
def direct_mean(x):
"""Среднее через последовательное суммирование."""
return direct_sum(x)/len(x)
print("Ошибка среднего для последовательного суммирования:",relative_error(exact_mean(),direct_mean(x)))
In [11]:
def direct_second_var(x):
"""Вторая оценка дисперсии через последовательное суммирование."""
return direct_mean(x**2)-direct_mean(x)**2
def online_second_var(x):
"""Вторая оценка дисперсии через один проход по выборке"""
m=x[0] # накопленное среднее
m2=x[0]**2 # накопленное среднее квадратов
for n in range(1,len(x)):
m=(m*(n-1)+x[n])/n
m2=(m2*(n-1)+x[n]**2)/n
return m2-m**2
print("Ошибка второй оценки дисперсии для последовательного суммирования:",relative_error(exact_variance(),direct_second_var(x)))
print("Ошибка второй оценки дисперсии для однопроходного суммирования:",relative_error(exact_variance(),online_second_var(x)))
In [12]:
def direct_first_var(x):
"""Первая оценка дисперсии через последовательное суммирование."""
return direct_mean((x-direct_mean(x))**2)
print("Ошибка первой оценки дисперсии для последовательного суммирования:",relative_error(exact_variance(),direct_first_var(x)))
Как мы видим, суммирование по первой формуле дает наиболее точный результат, суммирование по второй формуле менее точно, а однопроходная формула наименее точна.
Показательная функция имеет одно из самых простых разложений в ряд Тейлора:
$$e^x = \sum_{k=0}^\infty \frac{x^k}{k!}.$$Естественным желанием при решении задачи вычисления показательной функции является воспользоваться этим рядом. В данном разделе мы рассмотрим результативность этого подхода.
Так как на практике мы не можем суммировать бесконечное число слагаемых, то будем приближать ряд его частичной суммой:
$$e^x \approx \sum_{k=0}^N \frac{x^k}{k!}.$$Так как частичная сумма является многочленом, то для практического счета удобно воспользоваться (схемой Горнера)[ru.wikipedia.org/wiki/Схема_Горнера]:
$$e^x \approx 1+x\bigg(1+\frac{x}{2}\bigg(1+\frac{x}{3}\bigg(1+\frac{x}{4}\bigg(\ldots+\frac{x}{N}\bigg(1\bigg)\ldots\bigg)\bigg)\bigg)\bigg).$$Проведем эксперимент по оценки точности этого разложения. Сравнивать будем с библиотечной функцией numpy.exp, которая не дает совершенно точный ответ. Оценим погрешность библитечной функции, предполагая, что она вычисляется с максимальной возможной точностью. Число обусловленности показательной функции для относительной погрешности равно $\kappa_{exp}(x)=|x|$, тогда учитывая погрешности округления до числа с плавающей запятой, мы ожидаем предельную погрешность результата не менее $|x|\epsilon/2+\epsilon$.
In [13]:
def exp_taylor(x, N=None):
"""N-ая частичная сумма ряда Тейлора для экспоненты."""
acc = 1 # k-ая частичная сумму. Начинаем с k=0.
xk = 1 # Степени x^k.
inv_fact = 1 # 1/k!.
for k in range(1, N+1):
xk = xk*x
inv_fact /= k
acc += xk*inv_fact
return acc
def exp_horner(x, N=None):
"""N-ая частичная сумма ряда Тейлора для экспоненты методом Горнера."""
if N<=0: return 1 # Избегаем деления на ноль.
acc = 1 # Выражение во вложенных скобках в схеме Горнера
for k in range(N, 0, -1):
acc = acc/k*x + 1
return acc
def make_exp_test(fns, args={}, xmin=-1, xmax=1):
"""Проводит тест приближения fn показательной функции."""
x = np.linspace(xmin, xmax, 1000)
standard = np.exp(x)
theoretical_relative_error = (np.abs(x)/2+1)*np.finfo(float).eps
theoretical_absolute_error = theoretical_relative_error * standard
fig, ax1 = plt.subplots(1,1,figsize=(10,5))
ax2 = plt.twinx(ax1)
ax1.set_xlabel("Argument")
ax1.set_ylabel("Absolute error")
ax2.set_ylabel("Relative error")
ax1.semilogy(x, theoretical_absolute_error, '-r')
line, = ax2.semilogy(x, theoretical_relative_error, '--r')
line.set_label("theory (relative)")
for fn in fns:
subject = fn(x, **args)
absolute_error = np.abs(standard-subject)
relative_error = absolute_error/standard
ax1.semilogy(x, absolute_error, '-')
line, = ax2.semilogy(x, relative_error, '--')
line.set_label("{} (relative)".format(fn.__name__))
plt.legend()
plt.show()
make_exp_test([exp_taylor, exp_horner], args={"N": 3}, xmin=-0.001, xmax=0.001)
make_exp_test([exp_taylor, exp_horner], args={"N": 3}, xmin=-1, xmax=1)
make_exp_test([exp_taylor, exp_horner], args={"N": 3}, xmin=-10, xmax=10)
Ясно, что 4-x слагаемых слишком мало, чтобы хорошо приблизить ряд. Попробуем взять больше.
In [14]:
make_exp_test([exp_taylor, exp_horner], args={"N": 15}, xmin=-0.001, xmax=0.001)
make_exp_test([exp_taylor, exp_horner], args={"N": 15}, xmin=-1, xmax=1)
make_exp_test([exp_taylor, exp_horner], args={"N": 15}, xmin=-10, xmax=10)
Точность приближения растет с увеличением числа слагаемых, однако даже для умеренно больших аргументов ни одного верного знака в ответе не получается. Посмотрим, как погрешность изменяется в зависимости от числа слагаемых.
In [15]:
def cum_exp_taylor(x, N=None):
"""Вычисляет все частичные суммы ряда Тейлора для экспоненты по N-ую включительно."""
acc = np.empty(N+1, dtype=float)
acc[0] = 1 # k-ая частичная сумму. Начинаем с k=0.
xk = 1 # Степени x^k.
inv_fact = 1 # 1/k!.
for k in range(1, N+1):
xk = xk*x
inv_fact /= k
acc[k] = acc[k-1]+xk*inv_fact
return acc
x = -10
standard = np.exp(x)
theoretical_relative_error = (np.abs(x)/2+1)*np.finfo(float).eps
theoretical_absolute_error = theoretical_relative_error * standard
Ns = np.arange(100)
partial_sums = cum_exp_taylor(x, N=Ns[-1])
absolute_error = np.abs(partial_sums-standard)
relative_error = absolute_error/standard
fig, ax1 = plt.subplots(1,1,figsize=(10,5))
ax2 = plt.twinx(ax1)
ax1.set_xlabel("Argument")
ax1.set_ylabel("Absolute error")
ax2.set_ylabel("Relative error")
ax1.semilogy(Ns, Ns*0+theoretical_absolute_error, '-r')
line, = ax2.semilogy(Ns, Ns*0+theoretical_relative_error, '--r')
line.set_label("theory (relative)")
ax1.semilogy(Ns, absolute_error, '-')
line, = ax2.semilogy(Ns, relative_error, '--')
line.set_label("experiment (relative)")
plt.legend()
plt.show()
Оказывается, что даже суммируя очень большое число слагаемых, мы не можем достигнуть максимальной точности.