利用numpy广播/矢量化技术从其他阵列构建新的阵列

Question

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

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

当调用compute函数时，assets是资产名称的一维数组，close是一个多维numpy数组，其中assets中列出的每种资产都有窗口长度收盘价（使用相同的索引号），而industries是与一维数组中每个资产相关联的行业代码列表。在这个向量的计算中

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

结果表明，资产收益是每只股票所有计算收益的一维数组。我不清楚的一点是如何使用numpy完成计算，而不需要手动循环数组。在

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

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

在我看来，在numpy中使用一些经过高度优化的数组遍历应该是可能的，但是我似乎无法理解。在今天之前，我从未使用过numpy，而且我对Python相当陌生，所以这可能没什么用。在

---

更新：

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

^{pr2}$

它给了我以下错误：

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

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

上面的错误来自妈妈，蒙面的你在哪（）打电话。行业数组是20x8412，所以我假设资产收益是列在（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'

很明显，我做错了什么；但是翻看文件，我看不出是什么。这看起来应该是从平均值复制的（即np.浮动64dtype）到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_-cachedict，并避免所有随之发生的数据复制。因此，它不仅速度更快，而且内存效率也要高得多。在

---

更新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问题，其他人建议：

^{pr2}$

–用户36048