使用 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()

当计算函数被调用，资产是资产名称的1-d阵列，靠近是多维numpy的阵列，其中有window_length靠近价格中列出的每个资产的资产（使用相同的索引号），和工业是与一维数组中的每个资产相关联的行业代码列表。我知道 numpy 在这一行中矢量化了 Gainpct 的计算：

asset_gainpct = (close[-1] - close[offset]) / close[offset] * 100

结果是asset_gainpct是每只股票的所有计算收益的一维数组。我不清楚的部分是我将如何使用 numpy 来完成计算，而无需我手动循环遍历数组。

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

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

在我看来，这应该可以在 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：条件和输入之间的形状不一致（得到 (20, 8412) 和 (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() 调用产生了这个错误：

类型错误：无法根据规则“安全”将数组数据从 dtype('float64') 转换为 dtype('bool')

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

---

更新 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()

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

我仍然对任何人关于如何消除行业的最终循环并将该过程矢量化的任何想法感兴趣。我想知道 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 问题，其他人提出了以下建议：

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 中的循环，但它更简洁地描述了所需的步骤：

\n\n

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

\n\n

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

\n\n

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

\n\n

\xe2\x80\x93 用户36048

\n