Невозможно использовать interp1d из-за ошибки округления

Question

Я пытаюсь растянуть точки сетки, используя функции плотности. Например, учитывая следующее распределение точек (равномерно распределенное):

Приведенный ниже код изменит распределение на что-то вроде этого:

import numpy as np
from scipy.integrate import quad
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
from scipy.special import erf

# the x interval limits
a = 0.0
b = 3.0

# The density function, normalize it
_density_func = lambda x: 5*x
density_func = lambda x: _density_func(x) / quad(_density_func, a, b)[0]*(b-a)

# The xarr and yarr
Npts = 50
xarr = np.linspace(a, b, Npts)
yarr = np.zeros_like(xarr)

# Calculate the primitive function F = integral of density_func, from a, to t normalized by the int(density_func, a, b)
# F = np.vectorize(lambda t: quad(density_func, a, t)[0] / quad(density_func, a, b)[0])
F = np.vectorize(lambda t: quad(density_func, a, t)[0])

# if debug is true, print some info
debug = True
if debug:
    print('The range of xarr is: [', a, b, ']')
    print('The range of F(xarr) is: [', F(xarr).min(), F(xarr).max(), ']')

# Calculate the new x distribution of points using the inverse function of F.
# Use the trick of interpolation to calculate the inverse function,i.e: interp1d(y, x)(x)
xnew = interp1d(F(xarr), xarr)(xarr) 

# plot the result
plt.scatter(xnew, yarr, facecolors='none', edgecolors='black')
plt.show()

Когда я запускаю этот скрипт, я получаю следующую ошибку:

The range of xarr is: [ 0.0 3.0 ]
The range of F(xarr) is: [ 0.0 2.9999999999999996 ]
Traceback (most recent call last):
  File "meshDensity.py", line 38, in <module>
    xnew = interp1d(F(xarr), xarr)(xarr)
  File "C:\Users\navaro\Anaconda3\lib\site-packages\scipy\interpolate\polyint.py", line 79, in __call__
    y = self._evaluate(x)
  File "C:\Users\navaro\Anaconda3\lib\site-packages\scipy\interpolate\interpolate.py", line 664, in _evaluate
    below_bounds, above_bounds = self._check_bounds(x_new)
  File "C:\Users\navaro\Anaconda3\lib\site-packages\scipy\interpolate\interpolate.py", line 696, in _check_bounds
    raise ValueError("A value in x_new is above the interpolation "
ValueError: A value in x_new is above the interpolation range.

Как видите, проблема в правильном ограничении F(xarr)вместо точного значения 2.9999999999999996 3.0.

Не могли бы вы предложить какое-либо решение этой проблемы с ошибками округления? Я ценю любую помощь.

Редактировать : временное решение - использовать функцию mpmath.quad с mpmath.mp.dps = 20, но это делает скрипт относительно медленным.

Answer 1

На самом деле в вашем случае вам не нужно использовать mpmath модуль. Вам нужно ТОЛЬКО округлить верхний предел F(xarr) или просто заменить его верхним пределом xarr

a = 0.0
b = 0.0

F_xarr = F(xarr)
# replace the upper limit of F_xarr by b
arg_max = np.argmax(F_xarr)
# Or use round function for the upper limit
round(F_xarr[arg_max], 7)

Answer 2

Я решил свою проблему, используя модуль арифметики произвольной точности, mpmath.

import numpy as np
import mpmath as mp
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

mp.mp.dps = 18
# the x interval limits
a = 0.0
b = 3.0

# The density function, normalize it
_density_func = lambda x: 5*x
density_func = lambda x: _density_func(x) / mp.quad(_density_func, [a, b])*(b-a)

# The xarr and yarr
Npts = 50
xarr = np.linspace(a, b, Npts)
yarr = np.zeros_like(xarr)

# Calculate the primitive function F = integral of density_func, from a, to t.
F = np.vectorize(lambda t: mp.quad(density_func, [a, t]))

# if debug is true, print some info
debug = True
if debug:
    print('The range of xarr is: [', a, b, ']')
    print('The range of F(xarr) is: [', F(xarr).min(), F(xarr).max(), ']')

# Calculate the new x distribution of points using the inverse function of F.
# Use the trick of interpolation to calculate the inverse function,i.e: interp1d(y, x)(x)
xnew = interp1d(F(xarr), xarr)(xarr) 

# plot the result
plt.scatter(xnew, yarr, facecolors='none', edgecolors='black')
plt.show()

После запуска скрипта я получаю:

The range of xarr is: [ 0.0 3.0 ]
The range of F(xarr) is: [ 0.0 3.0 ]

Answer 3

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

Как видите, ошибка 4.440892098500626e-16, которая, вероятно, находится вокругТочность с плавающей запятой в вашей системе

Математический подход, вероятно, верен, но я бы посоветовал справиться с этим с помощью программного обеспечения, т.е. вам действительно нужны 16 десятичных знаков? Если вы добавите туда пару round(num,6), вы сделаете свой код более устойчивым к таким проблемам

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

density_func = lambda x: round(_density_func(x) / quad(_density_func, a, b)[0]*(b-a),6)

РЕДАКТИРОВАТЬ: как ни странно, я столкнулся с той же проблемой только сейчас: 3.0-2.9999999999999996 = 4.440892098500626e-16

, где мы оба знаем, что это должно было быть 0,0000000000000004

Забавные вещи случаются, когда вы имеете дело с числами на грани точности с плавающей запятой