Summary:
the current weighted-stacking implementation differs from the intended SNR-squared weighted stacking described in publications. Per confirmation from @ryanaloomis, the original codebase used an SNR-squared weighted stacking approach, but the function implemented in functions.py uses a different weighting scheme.
Details:
The intended formula for weighted stacking is:
$I_{\rm stacked} = \frac{\sum_i (I_i \times weight_i)}{\sum_i weight_i}$,
where the weights are defined as SNR squared, $weight_i = (\frac{pred\_ints_i}{stds_i})^2$, i.e. inverse variance weighting scaled by predicted signal intensity. This corresponds to the following portion in @ryanaloomis 's code:
total_weight = np.sum((pred_ints/stds)**2)
intensity = np.sum(interp_data*(pred_ints[:,np.newaxis]/stds[:,np.newaxis])**2, axis=0)/total_weight
In contrast, there are two stacking functions implemented in molsim,
velocity_stack in functions.py;velocity_stack in stack.py.
The version imported by ipython_quickstart.py (from functions.py) currently uses a different formula with inconsistent weights between numerator and denominator:
$I_{\rm stacked} = \frac{\sum_i (I_i \times \frac{max(I_i)}{max(max(I^{sim}_{i})} \times \frac{1}{rms^2_i})}{\sum_i rms^2_i}$,
corresponding to the code portion:
max_int = max(peak_ints)
for obs in obs_chunks:
if obs.flag is False:
obs.weight = obs.peak_int/max_int
if noise_arr is None:
obs.weight /= obs.rms**2
else:
obs.weight /= np.interp(obs.cfreq, obs.freq_obs, obs.rms)**2
obs.int_weighted = obs.int_obs * obs.weight
obs.int_sim_weighted = obs.int_sim * obs.weight
Here, obs.int_weighted is subsequently passed into obs.int_samp, which is then appended to interped_ints. The final averaging step is as follows:
rms_arr = ((~np.isnan(interped_ints)) * interped_rms**2).sum(axis=0)
int_avg = np.nansum(interped_ints,axis=0)/rms_arr
This introduces inconsistency between the code implementations and the description of “SNR-weighted stacking” used in published work. Proposed actions are to revise the formula implemented in functions.py to match the original SNR-squared weighted stacking definition.
Summary:
the current weighted-stacking implementation differs from the intended SNR-squared weighted stacking described in publications. Per confirmation from @ryanaloomis, the original codebase used an SNR-squared weighted stacking approach, but the function implemented in
functions.pyuses a different weighting scheme.Details:
$I_{\rm stacked} = \frac{\sum_i (I_i \times weight_i)}{\sum_i weight_i}$ ,$weight_i = (\frac{pred\_ints_i}{stds_i})^2$ , i.e. inverse variance weighting scaled by predicted signal intensity. This corresponds to the following portion in @ryanaloomis 's code:
The intended formula for weighted stacking is:
where the weights are defined as SNR squared,
In contrast, there are two stacking functions implemented in
molsim,velocity_stackinfunctions.py;velocity_stackinstack.py.The version imported by
$I_{\rm stacked} = \frac{\sum_i (I_i \times \frac{max(I_i)}{max(max(I^{sim}_{i})} \times \frac{1}{rms^2_i})}{\sum_i rms^2_i}$ ,
ipython_quickstart.py(fromfunctions.py) currently uses a different formula with inconsistent weights between numerator and denominator:corresponding to the code portion:
Here,
obs.int_weightedis subsequently passed intoobs.int_samp, which is then appended tointerped_ints. The final averaging step is as follows:This introduces inconsistency between the code implementations and the description of “SNR-weighted stacking” used in published work. Proposed actions are to revise the formula implemented in
functions.pyto match the original SNR-squared weighted stacking definition.