如何实现直方图相除并绘制结果？获取随红移变化的卫星星系占比

我已经绘制了三幅红移直方图：中心星系、卫星星系、总星系的红移分布，但尝试绘制随红移变化的卫星星系占比直方图时得到的结果毫无意义。以下是我使用的Python代码：

from astropy.io import fits
import matplotlib.pyplot as plt
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)
import matplotlib.ticker as ticker

hdulist = fits.open('./GalCat_lean-1.fits')
tbdata = hdulist[1].data
satellite_mask = tbdata['mass_halo'] == -1.0
sat_data = tbdata[satellite_mask]
central_mask = tbdata['mass_halo'] != -1.0
cen_data = tbdata[central_mask]

sat_z = sat_data.field(3)
cen_z = cen_data.field(3)

sat_masshalo = sat_data.field(0)
cen_masshalo = cen_data.field(0)

sat_Nsat = sat_data.field(6)
cen_Nsat = cen_data.field(6)


mass_halo = tbdata.field(0) 
RA = tbdata.field(1)
Dec = tbdata.field(2)
z_spec = tbdata.field(3)
halo_id = tbdata.field(4)
dx = tbdata.field(5)
N_sat = tbdata.field(6)
mag = tbdata.field(7)

#Central Redshift histogram:
## set binwidth 
binwidth = 0.01
hist1, bins, patches = plt.hist(cen_z, bins=np.arange(0., max(cen_z) + binwidth, binwidth), facecolor = 'blue')
plt.title('Central Redshift Histogram')
plt.xlabel("Redshift 'z'")
plt.ylabel("Population Density")
plt.show()

#Satellite Redshift histogram:
## set binwidth 
binwidth = 0.01
hist2, bins, patches = plt.hist(sat_z, bins=np.arange(0., max(sat_z) + binwidth, binwidth), facecolor = 'red')
plt.title('Satellite Redshift Histogram')
plt.xlabel("Redshift 'z'")
plt.ylabel("Population Density")
plt.show()

#Satellite Total Histogram
binwidth = 0.01
hist3, bins, patches = plt.hist(z_spec, bins=np.arange(0., max(sat_z) + binwidth, binwidth), facecolor = 'green')
plt.title('Total Redshift Histogram')
plt.xlabel("Redshift 'z'")
plt.ylabel("Population Density")
plt.show()
plt.hist(hist2/hist3, bins, facecolor='green')

你当前的核心错误是直接对直方图频数的比值hist2/hist3调用plt.hist()——这是在对比例值做二次直方图统计，而非按红移区间展示比例。正确做法是用完全统一的红移区间，计算每个区间内卫星星系占总星系的比例，再用柱状图/折线图展示。

修正步骤 #

• 统一所有统计的红移区间边界，确保卫星星系和总星系的区间完全对应
• 计算每个区间的卫星星系占比，处理除数为0的异常情况
• 直接绘制比例与红移区间的对应关系

修正后的代码 #

from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np  # 补充缺失的numpy导入

hdulist = fits.open('./GalCat_lean-1.fits')
tbdata = hdulist[1].data
satellite_mask = tbdata['mass_halo'] == -1.0
sat_data = tbdata[satellite_mask]
central_mask = tbdata['mass_halo'] != -1.0
cen_data = tbdata[central_mask]

sat_z = sat_data.field(3)
cen_z = cen_data.field(3)
z_spec = tbdata.field(3)

# 统一设置红移区间：基于总数据的范围生成bins
binwidth = 0.01
z_min = 0.0
z_max = max(z_spec)
bins = np.arange(z_min, z_max + binwidth, binwidth)

# 计算各区间的星系数量
hist_sat, _ = np.histogram(sat_z, bins=bins)
hist_total, _ = np.histogram(z_spec, bins=bins)

# 计算卫星星系占比：处理无星系的区间，避免NaN
sat_fraction = np.zeros_like(hist_total, dtype=float)
valid_bins = hist_total > 0
sat_fraction[valid_bins] = hist_sat[valid_bins] / hist_total[valid_bins]

# 绘制卫星占比随红移变化的柱状图
plt.figure(figsize=(10,6))
bin_centers = (bins[:-1] + bins[1:]) / 2  # 取每个区间的中点作为x轴坐标
plt.bar(bin_centers, sat_fraction, width=binwidth, color='purple', alpha=0.7)

plt.title('Satellite Galaxy Fraction vs Redshift')
plt.xlabel("Redshift 'z'")
plt.ylabel("Satellite Fraction (Sat/Total)")
plt.xlim(z_min, z_max)
plt.ylim(0, 1)  # 比例范围固定在0-1
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.show()

# 可选：同时展示总星系与卫星星系的分布对比
plt.figure(figsize=(10,6))
plt.hist(z_spec, bins=bins, facecolor='green', alpha=0.5, label='Total Galaxies')
plt.hist(sat_z, bins=bins, facecolor='red', alpha=0.5, label='Satellite Galaxies')
plt.title('Redshift Distribution: Total vs Satellite')
plt.xlabel("Redshift 'z'")
plt.ylabel("Population Density")
plt.legend()
plt.show()

关键修正说明 #

1. 补充缺失的numpy导入，避免np.arange运行报错
2. 基于总星系的红移范围生成统一区间，保证卫星星系与总星系的统计区间完全匹配
3. 用np.histogram直接获取频数数组，比plt.hist的返回值更适合后续计算
4. 处理空区间：当某个红移区间没有星系时，将占比设为0，避免绘图出现异常值
5. 用区间中点作为x轴坐标绘制柱状图，直观展示每个红移区间的卫星占比

内容的提问来源于stack exchange，提问作者Kazarama