python - 使用 numpy 广播/矢量化从其他数组构建新数组

Question

我正在研究 Quantopian 模型的股票排名因素。他们建议避免在自定义因子中使用循环。但是，我不确定在这种情况下如何避免循环。

def GainPctInd(offset=0, nbars=2):  
    class GainPctIndFact(CustomFactor):  
        window_length = nbars + offset  
        inputs = [USEquityPricing.close, ms.asset_classification.morningstar_industry_code]  
        def compute(self, today, assets, out, close, industries):
            # Compute the gain percents for all stocks
            asset_gainpct = (close[-1] - close[offset]) / close[offset] * 100  

            # For each industry, build a list of the per-stock gains over the given window  
            gains_by_industry = {}  
            for i in range(0, len(industries)):  
                industry = industries[0,i]  
                if industry in gains_by_industry:  
                    gains_by_industry[industry].append(asset_gainpct[i])  
                else:  
                    gains_by_industry[industry] = [asset_gainpct[i]]

            # Loop through each stock's industry and compute a mean value for that  
            # industry (caching it for reuse) and return that industry mean for  
            # that stock  
            mean_cache = {}  
            for i in range(0, len(industries)):  
                industry = industries[0,i]  
                if not industry in mean_cache:  
                    mean_cache[industry] = np.mean(gains_by_industry[industry])  
                out[i] = mean_cache[industry]  
    return GainPctIndFact()

调用计算函数时，assets 是 Assets 名称的一维数组，close 是多维 numpy 数组，其中 assets 中列出的每个 Assets 都有 window_length 收盘价(使用相同的索引号)，而行业是与一维数组中的每个 Assets 相关联的行业代码列表。我知道 numpy 在这一行中矢量化了 Gainpct 的计算:
asset_gainpct = (close[-1] - close[offset]) / close[offset] * 100
结果是 asset_gainpct 是每只股票的所有计算 yield 的一维数组。我不清楚的部分是我将如何使用 numpy 来完成计算，而无需我手动循环遍历数组。

基本上，我需要做的是根据所有股票所在的行业汇总所有股票的所有 yield ，然后计算这些值的平均值，然后将平均值分解为完整的 Assets 列表。

现在，我正在遍历所有行业并将 yield 百分比推送到行业索引字典中，该字典存储每个行业的 yield 列表。然后我计算这些列表的平均值并执行反向行业查找，以根据行业将行业 yield 映射到每个 Assets 。

在我看来，这应该可以在 numpy 中使用一些高度优化的数组遍历来完成，但我似乎无法弄清楚。我在今天之前从未使用过 numpy，而且我对 Python 还很陌生，所以这可能无济于事。

更新:

我修改了我的行业代码循环以尝试使用行业数组来使用屏蔽数组处理计算来屏蔽 asset_gainpct 数组，如下所示:

    # For each industry, build a list of the per-stock gains over the given window 
    gains_by_industry = {}
    for industry in industries.T:
        masked = ma.masked_where(industries != industry[0], asset_gainpct)
        np.nanmean(masked, out=out)

它给了我以下错误:

IndexError: Inconsistant shape between the condition and the input (got (20, 8412) and (8412,))

此外，作为旁注，行业以 20x8412 数组的形式出现，因为 window_length 设置为 20。额外值是前几天股票的行业代码，除非它们通常不会改变，因此它们可以是忽略。我现在正在迭代industries.T(行业的转置)，这意味着行业是一个 20 元素的数组，每个元素都具有相同的行业代码。因此，我只需要元素 0。

上面的错误来自 ma.masked_where() 调用。行业数组是 20x8412，所以我认为 asset_gainpct 是列为 (8412,) 的数组。我如何使这些与此调用兼容？

更新 2:

我再次修改了代码，解决了我遇到的其他几个问题。现在看起来像这样:

    # For each industry, build a list of the per-stock gains over the given window 
    unique_ind = np.unique(industries[0,])
    for industry in unique_ind:
        masked = ma.masked_where(industries[0,] != industry, asset_gainpct)
        mean = np.full_like(masked, np.nanmean(masked), dtype=np.float64, subok=False)
        np.copyto(out, mean, where=masked)

基本上，这里的新前提是我必须构建一个与输入数据中股票数量相同大小的均值填充数组，然后在应用我之前的掩码时将值复制到我的目标变量 (out) 中，以便只有未屏蔽的索引填充了平均值。此外，我意识到我在之前的化身中不止一次迭代行业，所以我也解决了这个问题。但是， copyto() 调用产生了这个错误:

TypeError: Cannot cast array data from dtype('float64') to dtype('bool') according to the rule 'safe'

显然，我做错了什么；但是查看文档，我不明白它是什么。这看起来应该是从 mean(它是 np.float64 dtype)复制到 out(我之前没有定义过)，它应该使用 masked 作为 bool 数组来选择要复制的索引。任何人对这个问题有什么想法？

更新 3:

首先，感谢所有贡献者的所有反馈。

在深入研究这段代码之后，我想出了以下几点:

def GainPctInd(offset=0, nbars=2):
    class GainPctIndFact(CustomFactor):
        window_length = nbars + offset
        inputs = [USEquityPricing.close, ms.asset_classification.morningstar_industry_code]
        def compute(self, today, assets, out, close, industries):
            num_bars, num_assets = close.shape
            newest_bar_idx = (num_bars - 1) - offset
            oldest_bar_idx = newest_bar_idx - (nbars - 1)

            # Compute the gain percents for all stocks
            asset_gainpct = ((close[newest_bar_idx] - close[oldest_bar_idx]) / close[oldest_bar_idx]) * 100

            # For each industry, build a list of the per-stock gains over the given window 
            unique_ind = np.unique(industries[0,])
            for industry in unique_ind:
                ind_view = asset_gainpct[industries[0,] == industry]
                ind_mean = np.nanmean(ind_view)
                out[industries[0,] == industry] = ind_mean
    return GainPctIndFact()

出于某种原因，基于屏蔽 View 的计算没有产生正确的结果。此外，将这些结果放入 out 变量是行不通的。沿着这条线的某个地方，我偶然发现了一篇关于 numpy(默认情况下)在执行切片时如何创建数组 View 而不是副本的帖子，并且您可以根据 bool 条件执行稀疏切片。在这样的 View 上运行计算时，就计算而言，它看起来像一个完整的数组，但所有值实际上仍然在基本数组中。这有点像有一个指针数组，计算发生在指针指向的数据上。同样，您可以为稀疏 View 中的所有节点分配一个值，并让它更新所有节点的数据。这实际上大大简化了逻辑。

我仍然对任何人关于如何消除行业的最终循环并将该过程矢量化的任何想法感兴趣。我想知道 map/reduce 方法是否可行，但我对 numpy 仍然不够熟悉，无法弄清楚如何比这个 FOR 循环更有效地做到这一点。从好的方面来说，剩余的循环只有大约 140 次迭代要通过，而前两个循环每个都要通过 8000 次。除此之外，我现在避免构建 gain_by_industry 和 mean_cache dict 并避免随之而来的所有数据复制。因此，它不仅速度更快，而且内存效率也更高。

更新 4:

有人给了我一个更简洁的方法来完成这个，最终消除了额外的 FOR 循环。它基本上将循环隐藏在 Pandas DataFrame groupby 中，但它更简洁地描述了所需的步骤:

def GainPctInd2(offset=0, nbars=2):
    class GainPctIndFact2(CustomFactor):
        window_length = nbars + offset
        inputs = [USEquityPricing.close, ms.asset_classification.morningstar_industry_code]
        def compute(self, today, assets, out, close, industries):
            df = pd.DataFrame(index=assets, data={
                    "gain": ((close[-1 - offset] / close[(-1 - offset) - (nbars - 1)]) - 1) * 100,
                    "industry_codes": industries[-1]
                 })
            out[:] = df.groupby("industry_codes").transform(np.mean).values.flatten()
    return GainPctIndFact2()

根据我的基准测试，它根本不会提高效率，但验证正确性可能更容易。他们的例子的一个问题是它使用 np.mean而不是 np.nanmean , 和 np.nanmean如果您尝试使用 NaN 值，则会删除导致形状不匹配的 NaN 值。为了解决 NaN 问题，其他人提出了以下建议:

def GainPctInd2(offset=0, nbars=2):
    class GainPctIndFact2(CustomFactor):
        window_length = nbars + offset
        inputs = [USEquityPricing.close, ms.asset_classification.morningstar_industry_code]
        def compute(self, today, assets, out, close, industries):
            df = pd.DataFrame(index=assets, data={
                    "gain": ((close[-1 - offset] / close[(-1 - offset) - (nbars - 1)]) - 1) * 100,
                    "industry_codes": industries[-1]
                 })
            nans = isnan(df['industry_codes'])
            notnan = ~nans
            out[notnan] = df[df['industry_codes'].notnull()].groupby("industry_codes").transform(np.nanmean).values.flatten()
            out[nans] = nan
    return GainPctIndFact2()

Answer 1

有人给了我一个更简洁的方法来完成这个，最终消除了额外的 FOR 循环。它基本上将循环隐藏在 Pandas DataFrame groupby 中，但它更简洁地描述了所需的步骤:

def GainPctInd2(offset=0, nbars=2):
    class GainPctIndFact2(CustomFactor):
        window_length = nbars + offset
        inputs = [USEquityPricing.close, ms.asset_classification.morningstar_industry_code]
        def compute(self, today, assets, out, close, industries):
            df = pd.DataFrame(index=assets, data={
                    "gain": ((close[-1 - offset] / close[(-1 - offset) - (nbars - 1)]) - 1) * 100,
                    "industry_codes": industries[-1]
                 })
            out[:] = df.groupby("industry_codes").transform(np.mean).values.flatten()
    return GainPctIndFact2()

根据我的基准测试，它根本不会提高效率，但验证正确性可能更容易。他们的例子的一个问题是它使用 np.mean而不是 np.nanmean , 和 np.nanmean如果您尝试使用 NaN 值，则会删除导致形状不匹配的 NaN 值。为了解决 NaN 问题，其他人提出了以下建议:

def GainPctInd2(offset=0, nbars=2):
    class GainPctIndFact2(CustomFactor):
        window_length = nbars + offset
        inputs = [USEquityPricing.close, ms.asset_classification.morningstar_industry_code]
        def compute(self, today, assets, out, close, industries):
            df = pd.DataFrame(index=assets, data={
                    "gain": ((close[-1 - offset] / close[(-1 - offset) - (nbars - 1)]) - 1) * 100,
                    "industry_codes": industries[-1]
                 })
            nans = isnan(df['industry_codes'])
            notnan = ~nans
            out[notnan] = df[df['industry_codes'].notnull()].groupby("industry_codes").transform(np.nanmean).values.flatten()
            out[nans] = nan
    return GainPctIndFact2()

– 用户 36048