计算从原点到R中每个时间点的梯形积分
Calculating trapezoidal integration from origin to each time point in R
我有一个数据框x:
head(x)
# time Qfr
#1 1 0.004751271
#2 2 0.005405618
#3 3 0.005785781
#4 4 0.006028213
#5 5 0.006179973
#6 6 0.006263814
我正在尝试计算从 time = 0 到每个时间点的数值积分，即积分：
\integral_{u=0}^t Qfr du
我的数据集看起来像
plot(x, type = "p", cex = 0.2)
到目前为止，我只能使用包 pracma:
计算总积分
require(pracma)
trapz(x$time, x$Qfr)
# [1] 0.1536843
如何对从原点到该行给定时间的积分进行编码？
非常感谢任何帮助！
x <-
structure(list(time = 1:100, Qfr = c(0.00475127142639315, 0.00540561802578535,
0.00578578141896237, 0.00602821304872631, 0.00617997318815436,
0.00626381438010966, 0.0062930341038365, 0.00627650284793016,
0.00622076547748955, 0.00613104312485634, 0.00601175416200995,
0.00586680072681021, 0.00569973138194467, 0.00551383427584607,
0.00531218958660475, 0.00509769744944577, 0.00487309097312275,
0.00464094029551979, 0.0044036514994002, 0.00416346290979426,
0.00392244046575488, 0.00368247330791138, 0.00344527034180023,
0.00321235826358148, 0.00298508133306843, 0.00276460302703881,
0.00255190958997126, 0.0023478154110241, 0.00215297008955578,
0.00196786700285879, 0.00179285315617775, 0.00162814007427384,
0.00147381548391774, 0.00132985553610085, 0.00119613732394456,
0.00107245146585054, 0.000958514542040229, 0.000853981195025623,
0.000758455729566888, 0.000671503074231956, 0.000592658993812166,
0.000521439468716574, 0.000457349183339767, 0.000399889089676022,
0.000348563034666554, 0.00030288345957139, 0.000262376196826176,
0.000226584404259194, 0.000195071688184451, 0.000167424475817082,
0.000143253703811788, 0.000122195893694217, 0.000103913686771705,
8.80959110345906e-05, 7.44572508696844e-05, 6.27375873853488e-05,
5.27010730706422e-05, 4.4134999643845e-05, 3.68485125344791e-05,
3.06712197100413e-05, 2.54517366998868e-05, 2.1056203851463e-05,
1.73668062169167e-05, 1.42803211212964e-05, 1.17067134905708e-05,
9.56779447711296e-06, 7.79595484853561e-06, 6.33298101979921e-06,
5.1289585088234e-06, 4.14126496950611e-06, 3.33365277945052e-06,
2.67541940141655e-06, 2.14066236056745e-06, 1.70761464370064e-06,
1.35805559020556e-06, 1.07679186575234e-06, 8.51202848467075e-07,
6.7084467545985e-07, 5.27107259938671e-07, 4.12918764032332e-07,
3.22492271674051e-07, 2.51109725048683e-07, 1.94938546369442e-07,
1.50876746740867e-07, 1.16422711484835e-07, 8.95662353428025e-08,
6.86977528360132e-08, 5.25330624541746e-08, 4.00511738714265e-08,
3.0443212344273e-08, 2.30705923615172e-08, 1.74309232147824e-08,
1.31303328068787e-08, 9.86109389078818e-09, 7.38361045310242e-09,
5.51197296231586e-09, 4.102421614833e-09, 3.044168569724e-09,
2.25212541292965e-09, 1.66116273443706e-09)), .Names = c("time",
"Qfr"), class = "data.frame", row.names = c(NA, -100L))
您需要遍历索引，在本例中为 i，并计算梯形直到第 i 个数据点。
sapply(1:100, function(i) trapz(x$time[1:i],x$Qfr[1:i]))
因为另一个答案告诉你如何使用 pracma::trapz 来达到你的目的，我不能那样做。我本来打算那样写一个答案，但由于我花了很多时间编辑你的问题，@shayaa 排在第一位。幸运的是，我有一个更好的主意。
梯形数值积分并不复杂。您已经在常规网格 1、2、3、... 100 上有 time，bin 大小为 1，以及网格上的已知函数值 Qfr。每个 bin 上的数值积分就是梯形的面积。所以，你可以计算：
## integration on each bin cell
cell <- with(x, (Qfr[1:99] + Qfr[2:100]) / 2)
## Note that precisely I should write
## cell <- with(x, (Qfr[1:99] + Qfr[2:100]) / 2 * diff(time))
## But as I said, you have equally spaced bin points with bin size 1
## `diff(time)` is always 1, hence left out
## You need to bear this in mind, once you work on more general cases.
那么，你想要的累计积分值就是：
cumsum(c(0, cell))
这个方法超级快！假设你有 N 个数据点，它的计算成本为 O(N)，但它是完全矢量化的。使用 sapply 的另一个答案未矢量化，将花费您 O(N^2) 计算。
我有一个数据框x:
head(x)
# time Qfr
#1 1 0.004751271
#2 2 0.005405618
#3 3 0.005785781
#4 4 0.006028213
#5 5 0.006179973
#6 6 0.006263814
我正在尝试计算从 time = 0 到每个时间点的数值积分，即积分：
\integral_{u=0}^t Qfr du
我的数据集看起来像
plot(x, type = "p", cex = 0.2)
到目前为止，我只能使用包 pracma:
require(pracma)
trapz(x$time, x$Qfr)
# [1] 0.1536843
如何对从原点到该行给定时间的积分进行编码？
非常感谢任何帮助！
x <-
structure(list(time = 1:100, Qfr = c(0.00475127142639315, 0.00540561802578535,
0.00578578141896237, 0.00602821304872631, 0.00617997318815436,
0.00626381438010966, 0.0062930341038365, 0.00627650284793016,
0.00622076547748955, 0.00613104312485634, 0.00601175416200995,
0.00586680072681021, 0.00569973138194467, 0.00551383427584607,
0.00531218958660475, 0.00509769744944577, 0.00487309097312275,
0.00464094029551979, 0.0044036514994002, 0.00416346290979426,
0.00392244046575488, 0.00368247330791138, 0.00344527034180023,
0.00321235826358148, 0.00298508133306843, 0.00276460302703881,
0.00255190958997126, 0.0023478154110241, 0.00215297008955578,
0.00196786700285879, 0.00179285315617775, 0.00162814007427384,
0.00147381548391774, 0.00132985553610085, 0.00119613732394456,
0.00107245146585054, 0.000958514542040229, 0.000853981195025623,
0.000758455729566888, 0.000671503074231956, 0.000592658993812166,
0.000521439468716574, 0.000457349183339767, 0.000399889089676022,
0.000348563034666554, 0.00030288345957139, 0.000262376196826176,
0.000226584404259194, 0.000195071688184451, 0.000167424475817082,
0.000143253703811788, 0.000122195893694217, 0.000103913686771705,
8.80959110345906e-05, 7.44572508696844e-05, 6.27375873853488e-05,
5.27010730706422e-05, 4.4134999643845e-05, 3.68485125344791e-05,
3.06712197100413e-05, 2.54517366998868e-05, 2.1056203851463e-05,
1.73668062169167e-05, 1.42803211212964e-05, 1.17067134905708e-05,
9.56779447711296e-06, 7.79595484853561e-06, 6.33298101979921e-06,
5.1289585088234e-06, 4.14126496950611e-06, 3.33365277945052e-06,
2.67541940141655e-06, 2.14066236056745e-06, 1.70761464370064e-06,
1.35805559020556e-06, 1.07679186575234e-06, 8.51202848467075e-07,
6.7084467545985e-07, 5.27107259938671e-07, 4.12918764032332e-07,
3.22492271674051e-07, 2.51109725048683e-07, 1.94938546369442e-07,
1.50876746740867e-07, 1.16422711484835e-07, 8.95662353428025e-08,
6.86977528360132e-08, 5.25330624541746e-08, 4.00511738714265e-08,
3.0443212344273e-08, 2.30705923615172e-08, 1.74309232147824e-08,
1.31303328068787e-08, 9.86109389078818e-09, 7.38361045310242e-09,
5.51197296231586e-09, 4.102421614833e-09, 3.044168569724e-09,
2.25212541292965e-09, 1.66116273443706e-09)), .Names = c("time",
"Qfr"), class = "data.frame", row.names = c(NA, -100L))
您需要遍历索引，在本例中为 i，并计算梯形直到第 i 个数据点。
sapply(1:100, function(i) trapz(x$time[1:i],x$Qfr[1:i]))
因为另一个答案告诉你如何使用 pracma::trapz 来达到你的目的，我不能那样做。我本来打算那样写一个答案，但由于我花了很多时间编辑你的问题，@shayaa 排在第一位。幸运的是，我有一个更好的主意。
梯形数值积分并不复杂。您已经在常规网格 1、2、3、... 100 上有 time，bin 大小为 1，以及网格上的已知函数值 Qfr。每个 bin 上的数值积分就是梯形的面积。所以，你可以计算：
## integration on each bin cell
cell <- with(x, (Qfr[1:99] + Qfr[2:100]) / 2)
## Note that precisely I should write
## cell <- with(x, (Qfr[1:99] + Qfr[2:100]) / 2 * diff(time))
## But as I said, you have equally spaced bin points with bin size 1
## `diff(time)` is always 1, hence left out
## You need to bear this in mind, once you work on more general cases.
那么，你想要的累计积分值就是：
cumsum(c(0, cell))
这个方法超级快！假设你有 N 个数据点，它的计算成本为 O(N)，但它是完全矢量化的。使用 sapply 的另一个答案未矢量化，将花费您 O(N^2) 计算。