加速插值练习

2023-12-14

我正在对大约 120 万个观测值运行大约 45,000 个局部线性回归(本质上),所以我希望得到一些帮助来加快速度,因为我很不耐烦。

我基本上是在尝试为一堆公司构建逐年的工资合同——工资函数(给定公司、年份、职位的经验)。

这是我正在使用的数据集(基本结构):

> wages
         firm year position exp salary
      1: 0007 1996        4   1  20029
      2: 0007 1996        4   1  23502
      3: 0007 1996        4   1  22105
      4: 0007 1996        4   2  23124
      5: 0007 1996        4   2  22700
     ---                              
1175141:  994 2012        5   2  47098
1175142:  994 2012        5   2  45488
1175143:  994 2012        5   2  47098
1175144:  994 2012        5   3  45488
1175145:  994 2012        5   3  47098

我想为所有公司构建经验级别 0 到 40 的工资函数,如下:

> salary_scales
        firm year position exp   salary
     1: 0007 1996        4   0       NA
     2: 0007 1996        4   1 21878.67
     3: 0007 1996        4   2 23401.33
     4: 0007 1996        4   3 23705.00
     5: 0007 1996        4   4 24260.00
    ---                                
611019: 9911 2015        4  36       NA
611020: 9911 2015        4  37       NA
611021: 9911 2015        4  38       NA
611022: 9911 2015        4  39       NA
611023: 9911 2015        4  40       NA

为此,我一直在努力(根据@BondedDust的建议here)与COBS(COstrained B-Spline) 包,它允许我构建工资合同的单调性。

一些问题仍然存在;特别是,当我需要进行推断时(只要给定的公司没有任何非常年轻或非常老的员工),拟合度就有可能失去单调性或降至 0 以下。

为了解决这个问题,我一直在数据范围之外使用简单的线性外推法——将拟合曲线扩展到外部min_exp and max_exp这样它就会通过两个最低(或最高)的拟合点——虽然不完美,但看起来效果不错。

考虑到这一点,这就是我到目前为止的做法(请记住,我是一个data.table狂热者):

#get the range of experience for each firm
wages[,min_exp:=min(exp),by=.(year,firm,position)]
wages[,max_exp:=max(exp),by=.(year,firm,position)]
#Can't interpolate if there are only 2 or 3 unique experience cells represented
wages[,node_count:=length(unique(exp)),by=.(year,firm,position)]
#Nor if there are too few teachers
wages[,ind_count:=.N,by=.(year,firm,position)]
#Also troublesome when there is little variation in salaries like so:
wages[,sal_scale_flag:=mean(abs(salary-mean(salary)))<50,by=.(year,firm,position)]
wages[,sal_count_flag:=length(unique(salary))<5,by=.(year,firm,position)]

cobs_extrap<-function(exp,salary,min_exp,max_exp,
                      constraint="increase",print.mesg=F,nknots=8,
                      keep.data=F,maxiter=150){
  #these are passed as vectors
  min_exp<-min_exp[1]
  max_exp<-min(max_exp[1],40)
  #get in-sample fit
  in_sample<-predict(cobs(x=exp,y=salary,
                          constraint=constraint,
                          print.mesg=print.mesg,nknots=nknots,
                          keep.data=keep.data,maxiter=maxiter),
                     z=min_exp:max_exp)[,"fit"]

  #append by linear extension below min_exp
  c(if (min_exp==1) NULL else in_sample[1]-
      (min_exp:1)*(in_sample[2]-in_sample[1]),in_sample,
    #append by linear extension above max_exp
    if (max_exp==40) NULL else in_sample[length(in_sample)]+(1:(40-max_exp))*
      (in_sample[length(in_sample)]-in_sample[length(in_sample)-1]))
}

salary_scales<-
  wages[node_count>=7&ind_count>=10
               &sal_scale_flag==0&sal_count_flag==0,
               .(exp=0:40,
                 salary=cobs_extrap(exp,salary,min_exp,max_exp)),
               by=.(year,firm,position)]

注意到有什么特别可能会减慢我的代码速度的事情吗?还是我被迫要有耐心?

这里可以使用一些较小的固定位置组合:

    firm year position exp salary count
 1: 0063 2010        5   2  37433    10
 2: 0063 2010        5   2  38749    10
 3: 0063 2010        5   4  38749    10
 4: 0063 2010        5   8  42700    10
 5: 0063 2010        5  11  47967    10
 6: 0063 2010        5  15  50637    10
 7: 0063 2010        5  19  51529    10
 8: 0063 2010        5  23  50637    10
 9: 0063 2010        5  33  52426    10
10: 0063 2010        5  37  52426    10
11: 9908 2006        4   1  26750    10
12: 9908 2006        4   6  36043    10
13: 9908 2006        4   7  20513    10
14: 9908 2006        4   8  45023    10
15: 9908 2006        4  13  33588    10
16: 9908 2006        4  15  46011    10
17: 9908 2006        4  15  37179    10
18: 9908 2006        4  22  43704    10
19: 9908 2006        4  28  56078    10
20: 9908 2006        4  29  44866    10

您的代码中有很多可以改进的地方,但让我们重点关注这里的主要瓶颈。当前的问题可以被视为尴尬地平行问题。这意味着您的数据可以分为多个较小的部分,每个较小的部分可以在单独的线程上单独计算,而无需任何额外的开销。

要了解当前问题的并行化可能性,您应该首先注意到您正在分别为每个单独的公司和/或年份执行完全相同的计算。例如,您可以将每年的计算拆分为较小的子任务,然后将这些子任务分配给不同的 CPU/GPU 核心。以这种方式可以获得显着的性能增益。 最后,当子任务处理完成后,您唯一需要做的就是合并结果。

但是,R 及其所有内部库作为单线程运行。您必须明确分割数据,然后将子任务分配给不同的核心。为了实现这一点,存在许多支持多线程的 R 包。我们将使用doparallel包在我们的示例中。

您没有提供足够大的显式数据集来有效测试性能,因此我们将首先创建一些随机数据:

set.seed(42)
wages<-data.table(firm=substr(10001:10010,2,5)[sample(10,size=1e6,replace=T)],
                  year=round(unif(1e6,1996,2015)),
                  position=round(runif(1e6,4,5)),
                  exp=round(runif(1e6,1,40)),
                  salary=round(exp(rnorm(1e6,mean=10.682,sd=.286))))
> wages
         firm year position exp salary
      1: 0001 1996        4  14  66136
      2: 0001 1996        4   3  42123
      3: 0001 1996        4   9  46528
      4: 0001 1996        4  11  35195
      5: 0001 1996        4   2  43926
     ---                              
 999996: 0010 2015        5  11  43140
 999997: 0010 2015        5  23  64025
 999998: 0010 2015        5  31  35266
 999999: 0010 2015        5  11  36267
1000000: 0010 2015        5   7  44315

现在,让我们运行代码的第一部分:

#get the range of experience for each firm
wages[,min_exp:=min(exp),by=.(year,firm,position)]
wages[,max_exp:=max(exp),by=.(year,firm,position)]
#Can't interpolate if there are only 2 or 3 unique experience cells represented
wages[,node_count:=length(unique(exp)),by=.(year,firm,position)]
#Nor if there are too few teachers
wages[,ind_count:=.N,by=.(year,firm,position)]
#Also troublesome when there is little variation in salaries like so:
wages[,sal_scale_flag:=mean(abs(salary-mean(salary)))<50,by=.(year,firm,position)]
wages[,sal_count_flag:=length(unique(salary))<5,by=.(year,firm,position)]
> wages
         firm year position exp salary min_exp max_exp node_count ind_count sal_scale_flag sal_count_flag
      1: 0001 1996        4  14  66136       1      40         40      1373          FALSE          FALSE
      2: 0001 1996        4   3  42123       1      40         40      1373          FALSE          FALSE
      3: 0001 1996        4   9  46528       1      40         40      1373          FALSE          FALSE
      4: 0001 1996        4  11  35195       1      40         40      1373          FALSE          FALSE
      5: 0001 1996        4   2  43926       1      40         40      1373          FALSE          FALSE
     ---                                                                                                 
 999996: 0010 2015        5  11  43140       1      40         40      1326          FALSE          FALSE
 999997: 0010 2015        5  23  64025       1      40         40      1326          FALSE          FALSE
 999998: 0010 2015        5  31  35266       1      40         40      1326          FALSE          FALSE
 999999: 0010 2015        5  11  36267       1      40         40      1326          FALSE          FALSE
1000000: 0010 2015        5   7  44315       1      40         40      1326          FALSE          FALSE

我们现在将处理wages像您之前所做的那样以单线程方式。注意,我们首先保存原始数据,以便稍后对其进行多线程操作并比较结果:

start <- Sys.time()
salary_scales_1 <-
  wages[node_count>=7&ind_count>=10
        &sal_scale_flag==0&sal_count_flag==0,
        .(exp=0:40,salary=cobs_extrap(exp,salary,min_exp,max_exp)),
        by=.(firm,year,position)]
print(paste("No Parallelisation time: ",Sys.time()-start))
> print(paste("No Parallelisation time: ",Sys.time()-start))
[1] "No Parallelisation time:  1.13971961339315"
> salary_scales_1
       firm year position exp   salary
    1: 0001 1996        4   0 43670.14
    2: 0001 1996        4   1 43674.00
    3: 0001 1996        4   2 43677.76
    4: 0001 1996        4   3 43681.43
    5: 0001 1996        4   4 43684.99
   ---                                
16396: 0010 2015        5  36 44464.02
16397: 0010 2015        5  37 44468.60
16398: 0010 2015        5  38 44471.35
16399: 0010 2015        5  39 44472.27
16400: 0010 2015        5  40 43077.70

处理完所有内容大约花了1分8秒。请注意,我们的虚拟示例中只有 10 家不同的公司,这就是为什么处理时间与本地结果相比并不那么重要。

现在,让我们尝试以并行方式执行此任务。如前所述,对于我们的示例,我们希望拆分每年的数据并将较小的子部分分配给单独的核心。我们将使用doParallel用于此目的的包:

我们需要做的第一件事是创建一个具有特定数量核心的集群。在我们的示例中,我们将尝试使用所有可用的内核。接下来,我们必须注册集群并将一些变量导出到子节点的全局环境中。在这种情况下,子节点只需要访问wages。此外,还需要在节点上评估一些依赖库才能使其正常工作。在这种情况下,节点需要访问两个data.frame and cobs图书馆。代码如下所示:

library(doParallel)
start <- Sys.time()
cl <- makeCluster(detectCores()); 
registerDoParallel(cl); 
clusterExport(cl,c("wages"),envir=environment());
clusterEvalQ(cl,library("data.table"));
clusterEvalQ(cl,library("cobs"));
salary_scales_2 <- foreach(i = 1996:2015) %dopar%
  {
    subSet <- wages[.(i)] # binary subsetting
    subSet[node_count>=7&ind_count>=10
           &sal_scale_flag==0&sal_count_flag==0,
           .(exp=0:40,
             salary=cobs_extrap(exp,salary,min_exp,max_exp)),
           by=.(firm,year,position)]
  }
stopCluster(cl)
print(paste("With parallelisation time: ",Sys.time()-start))
> print(paste("With parallelisation time: ",Sys.time()-start))
[1] "With parallelisation time:  23.4177722930908"

我们现在有一个数据表列表salary_scales_2其中包含每个单独年份的子结果。注意处理时间的加速:这次只用了 23 秒,而不是原来的 1.1 分钟(改善 65%)。我们现在唯一需要做的就是合并结果。我们可以用do.call("rbind", salary_scales_2)为了将表的行合并在一起(这几乎不需要时间 - 一次运行 0.002 秒)。最后,我们还执行一个小检查来验证多线程结果确实与单线程运行的结果相同:

salary_scales_2<-do.call("rbind",salary_scales_2)
identical(salary_scales_1,salary_scales_2)
> identical(salary_scales_1,salary_scales_2)
[1] TRUE

回复评论这确实是一个非常有趣的例子,但我认为您可能会忽略这里更重要的问题。这data.table确实执行内存和结构相关的优化,以便您以更有效的方式查询和访问数据。然而,在这个例子中,不存在与内存或搜索相关的主要瓶颈,尤其是当您与实际的总数据处理时间进行比较时。cobs功能。例如,您更改的行subSet <- wages[year==uniqueYears[i],]当您计时时,每次调用只需要大约 0.04 秒。

如果您在运行时使用分析器,那么您会注意到它不是data.table或其任何需要并行化的操作或分组,它是cobs占用几乎所有处理时间的函数(并且该函数甚至不使用data.table作为输入)。在这个例子中我们试图做的是重新分配我们的总工作负载cobs函数到不同的内核以实现加速。我们的意图是not来分割data.table操作,因为它们根本不昂贵。然而,我们确实必须拆分 data.table,因为我们需要将数据拆分为单独的数据。cobs函数运行。事实上,我们甚至利用了这样一个事实:data.table拆分和合并表时在各方面都很高效。这根本没有花费额外的时间。

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

加速插值练习 的相关文章

  • 按不规则时间间隔对数据进行分组求和(R语言)

    我正在看这里的 stackoverflow 帖子 R 计算一组内的观察次数 https stackoverflow com questions 65366412 r count number of observations within a
  • R 编程常用工具

    如果已经以不同的方式问过这个问题 我深表歉意 但我找不到任何达到我想要的东西 我真的是从其他软件包 SPSS 开始接触 R 的 当我了解真正可以做什么时 我意识到我还需要其他 工具 这让我想到了我的问题 您有哪些用于开发 R 代码的设置 我
  • 在 igraph 中为社区分配颜色

    我在 igraph 中使用 fastgreedy community 检测算法在 R 中生成社区 代码返回 12 个社区 但是在绘图时很难识别它们 因为它返回的图的颜色数量有限 我怎样才能用十二种不同的颜色绘制这个图表 l2 lt layo
  • C# 数据表来保存表格(无限嵌套)

    我相对较新C 但来自C C 背景 我需要一个类似于的数据类型 类 DataTable 但允许存储的列保存 简单 类型 int float boolean string 以及相同类型的数据 以便一个列可以保存另一个表 该表也具有存储表等的列
  • 更快的 %in% 运算符

    The 快速匹配 https cran r project org web packages fastmatch index html包实现了更快的版本match对于重复匹配 例如在循环中 set seed 1 library fastma
  • 排序因素与水平

    有人能解释一下 R 中 ordered 参数的用途吗 R says ordered逻辑标志来确定级别是否应被视为有序 按给定的顺序 所以如果我有一个名为名称的因素并设置ordered TRUE names lt factor c fred
  • WPF ListBox - 如何从数据表中放入值?

    I have ListBox并希望将值放入此列表框中DataTable listBoxVisibleFields DataContext SelectedFields Where SelectedFields is a DataTable充
  • 数据表中每一行的工具提示

    这个问题尖叫着是重复的JSF 2 0 Primefaces 2 x 数据表行的工具提示 https stackoverflow com questions 9980155 jsf 2 0 primefaces 2 x tooltip for
  • 如何使用 R 将每个文件的数据添加为附加行,从而将不同的 .csv 文件合并为一个完整的文件?

    我有几个不同的文件夹 它们都包含一个 csv 文件 所有这些 csv 文件都有一个单独的列 其中包含实验的一种条件的数据 我想以将每个文件的数据添加为新列的方式合并这些 csv 文件 目前 它看起来像这样 C1 csv 102 106 15
  • dplyr 返回每个组的全局平均值,而不是每个组的平均值

    有人可以解释一下我在这里做错了什么 library dplyr temp lt data frame a c 1 2 3 1 2 3 1 2 3 b c 1 2 3 1 2 3 1 2 3 temp gt group by temp 1 g
  • 将第 N 行上的 NA 行插入 data.frames 列表,其中 N 来自列表

    经过几个小时后 我发现自己无法解决以下问题 我有一个数据框列表 我想分别向每个 DF 插入 而不是替换 一行或多行 NA 始终至少一行 要插入的 NA 数量存储在单独的列表中 为了说明这一点 我有以下两个列表 list of datafra
  • 跟踪循环迭代

    抛硬币 成功 你赢100 否则你输50 你会一直玩 直到你口袋里有钱a 的价值如何a在任何迭代中都被存储 a lt 100 while a gt 0 if rbinom 1 1 0 5 1 a lt a 100 else a lt a 50
  • rpart“as.character(x) 中的错误:无法强制类型 'builtin' 为类型 'character' 的向量”消息是什么意思?

    我一直在用头撞rpart几天了 尝试为我拥有的这个数据集制作分类树 我认为现在是时候询问生命线了 我确信这是我没有看到的愚蠢的事情 但这里是我一直在做什么 EuropeWater lt read csv file paste Users a
  • R中的字典数据结构

    在 R 中 我有 例如 gt foo lt list a 1 b 2 c 3 如果我输入foo I get a 1 1 b 1 2 c 1 3 我怎样才能看透foo仅获取 键 列表 在这种情况下 a b c R 列表可以具有命名元素 因此可
  • 将字符串列拆分为多个虚拟变量

    作为 R 中 data table 包的相对缺乏经验的用户 我一直在尝试将一个文本列处理为大量指示符列 虚拟变量 每列中的 1 表示特定的子字符串是在字符串列中找到 例如我想处理这个 ID String 1 a b 2 b c 3 c 进入
  • 为什么数据帧上的 is.vector 不返回 TRUE?

    tl dr R 中的向量到底是什么 长版 R 中很多东西都是向量 例如 数字是长度为 1 的数值向量 is vector 1 1 TRUE 列表也是一个向量 is vector list 1 1 TRUE 好的 所以列表是一个向量 显然 数
  • 使用 R 下载压缩数据文件、提取和导入数据

    EZGraphs 在 Twitter 上写道 很多在线 csv 都被压缩了 有没有办法下载 解压缩存档并使用 R 将数据加载到 data frame Rstats 我今天也尝试这样做 但最终只是手动下载 zip 文件 我尝试过类似的东西 f
  • Quantmod 的简单功能不再起作用

    我明天要交论文 我收到了一条关于 quantmod 的非常奇怪的错误消息 这是我在过去几周使用这个包时从未遇到过的 我无法导入特定于道琼斯指数 DJI 的数据 我收到以下错误消息 getSymbols DJI src yahoo from
  • R 中两个时间戳之间的左连接

    我的目标是执行左连接intervals哪里的bike id比赛和created at时间戳在records在 之间start and end in the intervals table gt class records 1 data ta
  • 在R中循环子文件夹

    我正在 R 环境中包含多个子文件夹的文件夹中工作 我想要循环遍历多个子文件夹 然后在每个子文件夹中调用 R 脚本来执行 我想出了下面的代码 但我的代码似乎添加了 到子文件夹列表 我收到错误 文件中的错误 文件名 r 编码 编码 无效的 描述

随机推荐

  • 安装 debug_inspector (0.0.2) 时出错,Bundler 无法继续

    将我的 Rails 应用程序从 3 2 迁移到 4 0 0 并将 ruby 2 2 2 迁移到 2 2 4 在进行捆绑安装时遇到 安装 debug inspector 0 0 2 时发生错误 并且 Bundler 无法继续 问题 使用 ub
  • 自定义ListCellRenderer不会改变背景颜色

    我有这门课 SuppressWarnings serial private class DataCellRenderer extends JLabel implements ListCellRenderer public DataCellR
  • 为什么 String.intern() 在 Oracle JDK 1.7 中的行为不同?

    这是一个java片段 public class TestIntern public static void main String argvs String s1 new StringBuilder ja append va toStrin
  • 使用 Delphi 2007 将 Base64 字符串作为二进制文件保存到磁盘

    我有一个 Base64 二进制字符串 它是由第 3 方供应商发送给我们的 XML 文档的一部分 我希望能够将其保存回其原始文件格式 jpg 使用此问题中接受的答案 使用 php 将 Base64 字符串作为二进制文件保存到磁盘 我可以轻松地
  • 是否可以被动安装 .EXE 但仍使用 Powershell 显示 GUI?

    正如标题所说 是否可以使用 Powershell 被动 静默安装 EXE 但仍然显示安装程序 GUI 我希望下一个自动 单击 但仍然希望 GUI 显示为进度指示器 UPDATE 有一个用于 Windows Installer 的 Power
  • jQuery 通过 XPath 选择元素

    我有一个 XPath 选择器 如何使用 jQuery 获取与该选择器匹配的元素 我见过https developer mozilla org en Introduction to using XPath in JavaScript但它没有使
  • VS Code Markdown - 有没有一种聪明的方法来使链接智能感知?

    我开始使用 VisualStudio Code 和 Markdown 作为笔记工具 类似于泽特卡斯滕 作为此笔记系统的一部分 笔记应始终链接到其他笔记 文件 作为正在恢复的开发人员 我开始写 决策技术 当我输入 然后按 ctrl space
  • 将 2d 数组转换为 1d PHP

    我有一个数组 下面的示例 其中包含 3 个内容非常相似的数组 我想知道是否有一种简单的方法可以将所有内容放入一个数组中 而不是 3 个单独的数组中 我尝试的所有操作似乎都只是覆盖数据 而我只留下了第三个数组中的信息 array 0 gt a
  • 每个产品返回一个 SQL 行,其中包含价格和最新日期

    我不是一个 SQL 专家 我过去使用过它 很少遇到 google 无法解决的问题 但是这次我需要询问社区 我有一个数据库 其中有一个名为 交易 的表 它的数据如下 ProdNo Price TransactionDate Purchased
  • 如何更改 Activity 中的片段文本视图文本

    我找不到如何更改片段的textview从一个Activity 我有 4 个文件 MainActivity java activity main xml FragmentClass java frag class xml frag class
  • 开发证书推送通知

    如何向在开发配置文件下运行的构建发送推送通知 谁能帮帮我吗 您应该使用您的开发版本证书在 Apple Push Notification Center 上进行身份验证 要测试推送通知 您可以使用 MacOs 应用程序 推我宝贝 完整的使用教
  • pandas 日期转字符串

    我有一个约会时间pandas Series 一栏称为 日期 我想在像字符串一样的循环中获取 i 元素 s apply lambda x x strftime Y m d or astype str tail 1 reset index da
  • Java 数据库连接

    我有一个java项目 其中有许多连接到数据库的文件 谁能告诉我是否可以使用 java 类文件连接到数据库 这样我就不会为每个文件创建数据库连接 请教我如何操作 感谢您的帮助 D 这是我使用的代码 但它不起作用 dbConnect java
  • 在 Django 中,如何向用户发送登录后消息?

    我有一个 Django 网站 当用户登录时 我会向他们显示一条成功消息 我使用如下信号来执行此操作 def post login actions sender user request kwargs messages success req
  • 使用 OpenSL 进行通话录音

    自从 Galaxy S5 中的棒棒糖更新后 我尝试修复应用程序中的通话录音 作为基础 我使用这里的谷歌示例项目 Sample 这是代码的主要部分 AudioRecorder AudioRecorder SampleFormat sample
  • 使用D3多次更新表不起作用

    我有一个页面需要加载新数据并使用 d3 将其更新为现有表元素 我正在写以下内容但没有用 completeTableData head1 va1 head2 val2 haed1 val4 head2 val5 var tr d3 selec
  • Android 上的新电子邮件通知

    有没有办法让应用程序在收到新电子邮件时收到通知 我的目标是 如果收到具有特定主题的电子邮件 将其视为 过滤器 或 规则 则发出警报 播放声音 振动等 但我不想检查电子邮件我自己做服务器 我想我正在寻找类似 android telephony
  • 如何在android中从图库中选择多张图像?

    我正在制作一个项目 我想从图库中选择多张照片并将其保存在图像视图数组中 我可以导入单个图像并保存在 imageview 谁能告诉我如何导入多个图像并保存在数组或不同的图像视图中 MainActivity extends Activity i
  • 点击式 tkinter 窗口

    该函数复制自使用 TkInter 绑定设置不可交互 点击 覆盖 不仅窗口无法点击 PNG 也不透明 PNG在这里 https drive google com file d 1tlLl2hjPq38mc c PpMhkKDlP1HqvDY5
  • 加速插值练习

    我正在对大约 120 万个观测值运行大约 45 000 个局部线性回归 本质上 所以我希望得到一些帮助来加快速度 因为我很不耐烦 我基本上是在尝试为一堆公司构建逐年的工资合同 工资函数 给定公司 年份 职位的经验 这是我正在使用的数据集 基