使用 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 是资产名称的一维数组，close 是一个多维 numpy 数组，其中有 window_length assets 中列出的每个资产的收盘价（使用相同的指数编号），industries 是与 1 中每个资产相关联的行业代码列表-d 数组。我知道 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,))

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

上面的错误来自 ma.masked_where() 调用。 industries 数组是 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 em>) 同时应用我以前的掩码，以便只有未掩码的索引填充平均值。此外，我意识到我在之前的版本中不止一次地迭代了行业，所以我也修复了这个问题。但是，copyto() 调用产生了这个错误：

TypeError: 无法将数组数据从 dtype('float64') 转换为 dtype('bool') 根据规则 'safe'

显然，我做错了什么；但是查看文档，我看不到它是什么。这看起来应该是从 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 次。除此之外，我现在正在避免构造 gains_by_industry 和 mean_cache 字典，并避免随之而来的所有数据复制。因此，它不仅速度更快，而且内存效率也更高。

---

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

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

--user36048