python - 由于舍入错误，无法使用 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.

如您所见，问题在于 is 的正确限制F(xarr)而2.9999999999999996不是确切值3.0。

您能否建议任何解决此舍入错误问题的方法？我很感激任何帮助。

编辑：临时解决方案是使用mpmath.quad函数，mpmath.mp.dps = 20但这会使脚本相对较慢。

标签： pythonnumpyinterpolation

Answer 1

我使用任意精度算术模块解决了我的问题，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 ]