提问人：ClimateUnboxed 提问时间：8/16/2023 最后编辑：ClimateUnboxed 更新时间：8/17/2023 访问量：67
使用直方图将变量 Y 采样为 X 的函数时，python 舍入错误？
python rounding error when sampling variable Y as a function of X with histogram?
问：
我正在尝试使用函数直方图将变量 （SST） 采样为另一个变量 （TCWV） 的函数，并将权重设置为采样变量，如下所示：
# average sst over bins
num, _ = np.histogram(tcwv, bins=bins)
sstsum, _ = np.histogram(tcwv, bins=bins,weights=sst)
out=np.zeros_like(sstsum)
out[:]=np.nan
sstav = np.divide(sstsum,num,out=out, where=num>100)
可重复性的完整代码如下。我的问题是，当我绘制原始数据的散点图，然后绘制计算出的平均值时，平均值位于数据“云”之外，如下所示（见右图）：
我想不出为什么会发生这种情况，除非可能是四舍五入错误？
这是我的全部代码：
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
# if you have a recent netcdf libraries you can access it directly here
url = ('http://clima-dods.ictp.it/Users/tompkins/CRM/data/WRF_1min_mem3_grid4.nc#mode=bytes')
ds=Dataset(url)
### otherwise need to download, and use this:
###ifile="WRF_1min_mem3_grid4.nc"
###ds=Dataset(idir+ifile)
# axis bins
bins=np.linspace(40,80,21)
iran1,iran2=40,60
# can put in dict and loop
sst=ds.variables["sst"][iran1:iran2+1,:,:]
tcwv=ds.variables["tcwv"][iran1:iran2+1,:,:]
# don't need to flatten, just tried it to see if helps (it doesn't)
sst=sst.flatten()
tcwv=tcwv.flatten()
# average sst over bins
num, _ = np.histogram(tcwv, bins=bins)
sstsum, _ = np.histogram(tcwv, bins=bins,weights=sst)
out=np.zeros_like(sstsum)
out[:]=np.nan
sstav = np.divide(sstsum,num,out=out,where=num>100)
# bins centroids
avbins=(np.array(bins[1:])+np.array(bins[:-1]))/2
#plot
subsam=2
fig,(ax)=plt.subplots()
plt.scatter(tcwv.flatten()[::subsam],sst.flatten()[::subsam],s=0.05,marker=".")
plt.scatter(avbins,sstav,s=3,color="red")
plt.ylim(299,303)
plt.savefig("scatter.png")
答：
4赞
Nick ODell
8/17/2023
#1
我想不出为什么会发生这种情况，除非可能是四舍五入错误？
这实际上是一个舍入错误。
具体来说，当您在此处计算每个 bin 中的 sst 总和时：
sstsum, _ = np.histogram(tcwv, bins=bins,weights=sst)
结果错误了 0.1%，而我尝试了两种替代方法来计算总和。
对于解决此问题的方法，我有两个想法。
方法#1
最简单的解决方法是更精确地进行计算。
sstsum, _ = np.histogram(tcwv, bins=bins,weights=sst.astype('float64'))
如果不进行此更改，则 sst 的 dtype 为 float32。
方法#2
出于性能原因，您可能希望将计算保留在 32 位浮点数中。它们比 64 位浮点数快一些。另一种解决方案是在求和之前减去平均值以提高数值稳定性。
sst_mean = sst.mean()
num, _ = np.histogram(tcwv, bins=bins)
sstsum, _ = np.histogram(tcwv, bins=bins,weights=sst - sst_mean)
out=np.zeros_like(sstsum)
out[:]=np.nan
sstav = np.divide(sstsum,num,out=out,where=num>100)
sstav += sst_mean
这将从每个数据点中减去 sst 的总体平均值，然后在末尾将其加回。由于浮点数在 0 附近具有更高的精度，这使得计算更加精确。
比较
以下是方法 #1 的图：
方法 #2 的图看起来相同。这两种方法相等，相距 1.32 * 10-5 以内。
评论