我尝试使用函数直方图将一个变量（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")
1条答案
按热度按时间lawou6xi1#
我不知道为什么会发生这种情况，除非这是一个舍入误差？
这实际上是一个舍入误差。
具体来说，当你在这里计算每个bin中的sst之和时：
与我尝试计算总和的两种替代方法相比，结果错误了0.1%。
我有两个想法来解决这个问题。
方法一
最简单的解决方法是更精确地进行计算。
如果不做此更改，sst的dtype为float 32。
方法二
出于性能考虑，您可能希望将计算保持在32位浮点数中。它们比64位浮点数快一些。另一种解决方案是在求和之前减去平均值，以提高数值稳定性。
这是从每个数据点中减去SST的总体平均值，然后在最后将其加回。由于浮点数在0附近具有更高的精度，这使得计算更精确。
对比
下面是方法#1的图：
方法#2的图看起来是相同的。这两种方法在1.32 * 10-5的范围内相等。