代码之家 › 专栏 › 技术社区 › bdecaf

如何在分布不均的数据集中拟合平滑滞后?

smoothing spline r

1

bdecaf · 技术社区 · 6 年前

How can I fit a smooth hysteresis in R? . 直接应用 smooth.spline 虽然它被证明是解决这类问题的一个有用的通用思想,并且在我的模拟玩具数据集上工作得很好,但它却很不符合我的实际情况。

here .

下面的图像是由最后的代码创建的。

Thysteresis还没有完全关闭-但这不是一个问题。应用 平滑.spline 产生的方式太乱输出(红色)。使用参数 spar 在里面 平滑.spline

然而真正让我不安的是左下角点周围的“结”——即使是平滑曲线(蓝色)也是如此。当它从左下角移动到右上角时,它还有一个人工“循环”,尽管附近没有数据点。

library(dplyr)
library(ggplot2)
library(tidyr)

exdata <- readRDS("data/example.Rdata")

exdata %<>% 
  mutate(y_smooth = smofun(y),
         x_smooth = smofun(x),
         y_smooth2 = smofun(y, spar = 0.7),
         x_smooth2= smofun(x, spar = 0.7)
         )

ggplot(exdata, mapping = aes(x=x, y=y)) +
  geom_point(alpha = 0.3) +
  geom_path(mapping = aes(x=x_smooth, y=y_smooth), colour = "red") +
  geom_path(mapping = aes(x=x_smooth2, y=y_smooth2), colour = "blue") +
  theme_light()

1 回复 | 直到 6 年前

1

Zheyuan Li 6 年前

正如我在上一篇文章中所说: How can I fit a smooth hysteresis in R? ,我们要分别平滑每个坐标。

## you have 713 data
plot(1:713, exdata$x); points(1:713, exdata$y)

smooth.spline 每个人。它在你最初的玩具数据集上运行得很好,但在你上面绘制的真实数据集上运行得不太好。直接使用 R smooth.spline(): smoothing spline is not smooth but overfitting my data

平滑.spline 不允许我们手动放置结。你可以告诉它你想要多少节(通过使用参数) nknots )它通过分位数在内部放置节点。事实上,我们不必依赖 平滑.spline . 这只是一个回归问题。我们可以用任何一种回归方法,只要你满意。在下文中,我使用由提供的回归B样条曲线 splines 包(一个带有R,所以不需要安装)和适合使用普通最小二乘法(通过 lm ).

library(splines)

test <- function (b1, b2, n1, n2, n3, degree = 1, exdata) {

  x <- exdata$x; y <- exdata$y; t <- 1:length(x)

  ## place knots:
  ## n1 knots on [1, b1)
  ## n2 knots on [b1, b2]
  ## n3 knots on (b3, 713]
  t1 <- seq(1, b1, length = n1 + 1)
  t2 <- seq(b1, b2, length = n2)
  t3 <- seq(b2, 713, length = n3 + 1)
  kt <- c(t1[1:n1], t2, t3[-1])

  kt_inner <- kt[2:(length(kt) - 1)]

  fitx <- lm( x ~ bs(t, degree = degree, knots = kt_inner) )
  fity <- lm( y ~ bs(t, degree = degree, knots = kt_inner) )

  xh <- fitted(fitx); yh <- fitted(fity)

  par(mfrow = c(1, 2))

  plot(t, x, col = 8); points(t, y, col = 8)
  abline(v = kt, lty = 2, col = 8)
  lines(t, xh, lwd = 2); lines(t, yh, lwd = 2)

  plot(x, y, col = 8)
  lines(xh, yh, lwd = 2)

  z <- list(x = xh, y = yh)
  }

test函数允许您指定两个断点 b1 b2 [1, b1) n1 结,开 [b1, b2] n2 结,等等 (b3, 713] 你想要什么 n3 结。设置 degree 允许我们尝试线性样条、二次样条和三次样条。测试函数将生成拟合的样条曲线,并以不可见的方式返回平滑的样条曲线 x , y 价值观。

经过几次尝试,我发现以下可能是最好的:

z <- test(b1 = 90, b2 = 130, n1 = 4, n2 = 5, n3 = 10, degree = 3, exdata)

这里没有处罚。但是,通过使用 mgcv::gam 代替 流光溢彩 + splines::bs mgcv 允许我们指定结(请参见: mgcv: How to set number and / or locations of knots for splines ).

如果我们不想麻烦地指定结,我们可以使用“自适应样条线”平滑在 ,它允许平滑参数局部变化(它本质上产生一组平滑参数)。这确实克服了 平滑.spline 其中只有一个控制全局惩罚的平滑参数。

请注意,在样条线方法中,选择断点隐藏在节点的选择中,因为断点在节点集中。

也没有基于样条线的方法。包裹 segmented 分段线性回归是通过自动选择变化点/断点来实现的(但是你必须通过使用参数提前告诉它有多少断点) psi

不同的方法有很大的比较和讨论空间。作为将来可能重新调查的备份,我在这里附上了您的数据集(以防dropbox链接中断)。

exdata<- 124.6366, 123.3217, 127.4646, 128.0695, 110.6105, 105.3012, 98.3239, 89.5727, 101.3341, 103.7006, 100.9025, 101.2322, 107.0066, 99.4793, 101.446, 110.8227, 115.8369, 113.8645, 119.1815, 115.382, 106.6097, 99.0487, 99.2332, 97.9979, 98.5008, 99.4881, 96.2035, 99.3862, 108.1677, 107.2952, 107.5039, 109.8523, 107.0338, 98.331, 106.3463, 96.0376, 95.0978, 92.9253, 94.8009, 96.7117, 101.732, 107.1491, 107.0559, 107.6543, 108.8192, 101.6353, 94.3223, 94.3889, 93.4375, 99.0691, 101.1942, 102.7105, 105.8936, 104.8027, 93.6381, 99.6715, 98.8952, 100.2534, 104.2578, 105.9263, 106.0875, 118.8101, 121.278, 535.1731, 276.7803, 330.6044, 439.8171, 541.2278, 648.0471, 794.4545, 865.3579, 909.2367, 978.2347, 1067.6434, 1088.8355, 1153.1032, 1301.0525, 1300.5524, 1319.1291, 1332.0036, 1339.5618, 1342.3142, 1351.3678, 1355.4577, 1364.5274, 1386.3259, 1388.0149, 1393.6825, 1407.7905, 1410.5657, 1389.9132, 1386.5312, 1376.0732, 1371.785, 1375.1907, 1387.2856, 1386.5125, 1387.6686, 1379.9445, 1378.8248, 1346.6911, 1352.2614, 1347.0287, 1349.8785, 1349.8588, 1345.208, 1335.5, 1343.194, 1331.6572, 1323.6614, 1314.4852, 1322.7351, 1320.4661, 1329.9164, 1326.9096, 1336.5571, 1333.9171, 1330.2649, 1293.1904, 1292.3817, 1281.3893, 1285.0823, 1292.0862, 1294.0099, 1286.4846, 1297.9813, 1289.5944, 1286.8668, 1282.2538, 1275.3989, 1273.0393, 1275.9338, 1272.8439, 1270.2929, 1263.9435, 1253.0364, 1247.754, 1250.2937, 1259.3518, 1273.8355, 1286.2065, 1293.7979, 1286.6026, 1290.399, 1285.4735, 1273.8435, 1261.1018, 1264.7923, 1250.2116, 1245.0338, 1247.965, 1253.4547, 1248.2883, 1250.0202, 1254.8524, 1257.7671, 1263.4184, 1261.9403, 1260.0229, 1258.0052, 1232.7507, 1240.5921, 1235.1443, 1230.6716, 1254.277, 1244.1764, 1236.4669, 1253.6932, 1263.0974, 1248.3615, 1255.2735, 1257.1036, 1251.1783, 1247.3818, 1249.0185, 1238.181, 1232.3126, 1224.0232, 1215.3999, 1207.859, 1207.6245, 1204.4897, 1212.6945, 1206.7439, 1211.7789, 1211.6164, 1221.6207, 1226.8173, 1224.421, 1223.8749, 1217.5781, 1204.7366, 1192.8738, 1196.071, 1195.3558, 1198.7697, 1196.0119, 1191.8406, 1181.6458, 1184.4447, 1180.3189, 1187.6651, 1193.5248, 1189.1359, 1185.1723, 1172.775, 1170.4184, 1175.5543, 1128.5597, 1131.7037, 1126.5585, 1128.371, 1141.5219, 1149.2273, 1136.943, 1143.498, 1148.9494, 1147.0205, 1151.1578, 1165.1582, 1157.0784, 1155.4953, 1149.6255, 1142.4007, 1128.9415, 1129.4853, 1136.5361, 1136.5032, 1138.8941, 1144.9796, 1146.6936, 1141.7335, 1146.273, 1139.3171, 1144.506, 1142.9864, 1144.5049, 1147.8466, 1152.4488, 1136.4797, 1134.9687, 1132.3612, 1131.8932, 1129.8746, 1135.4942, 1137.3314, 1135.515, 1123.3012, 1123.8599, 1127.5974, 1121.888, 1121.8971, 1122.7738, 1121.3582, 1101.4597, 1101.6568, 1106.9945, 1111.0663, 1109.69, 1123.6524, 1117.9505, 1113.164, 1102.1349, 1106.6481, 1114.2183, 1117.4973, 1115.0706, 1105.0425, 1098.5148, 1095.9327, 1091.0444, 1094.8351, 1097.9752, 1102.8526, 1101.8913, 1104.7084, 1105.5634, 1094.34, 1098.789, 1101.8951, 1093.274, 1082.6922, 1088.4664, 1088.8069, 1081.3171, 1082.6755, 1081.0162, 1074.4439, 1072.8442, 1082.235, 1090.5553, 1100.5808, 1095.3863, 1098.7245, 1097.1755, 1092.7324, 1088.4808, 1088.091, 1082.7414, 1084.135, 1087.6769, 1080.6349, 1079.5889, 1073.2832, 1072.3528, 1063.9176, 1074.6075, 1068.1074, 1071.7748, 1073.2347, 1081.5981, 1071.6388, 1072.2588, 1064.7244, 1055.2577, 1044.7628, 1040.6255, 1047.1508, 1051.5531, 1050.5133, 1063.7744, 1070.2546, 1040.4452, 1051.2625, 1064.2338, 1067.7722, 1059.7189, 1065.0454, 1062.3471, 1059.7083, 1055.7133, 1057.4425, 1057.8175, 1052.1135, 1059.3531, 1061.0893, 1071.2829, 1074.8224, 1080.928, 1068.002, 1071.509, 1070.4365, 1071.788, 1066.6481, 1068.2753, 1071.1143, 1049.311, 1056.0742, 1042.2977, 1041.6723, 1043.9296, 1042.7619, 1032.842, 1042.7006, 1046.3895, 1051.801, 1054.1023, 1059.6545, 1051.8906, 1052.3833, 1055.4235, 1051.0044, 1051.7467, 1051.2107, 1035.2799, 1034.5771, 1032.9374, 1037.5497, 1041.925, 1043.5744, 1040.8266, 1041.4861, 1036.9859, 1024.7875, 1020.2694, 1017.749, 1014.9565, 1011.8741, 1015.188, 1018.6312, 1009.4944, 1010.4238, 1017.4137, 1019.2012, 1021.8013, 1028.7364, 1023.8553, 1014.1613, 1023.4766, 1015.2452, 1015.9395, 1018.8826, 1016.6993, 1013.6742, 1019.1466, 1021.0248, 1026.475, 1027.4478, 1034.2363, 1032.2496, 1027.5767, 1020.5724, 1019.5044, 1004.9575, 997.7392, 1001.0576, 1019.6292, 1018.1884, 1021.1877, 1023.7156, 1025.1062, 1017.3725, 1013.8538, 1009.8373, 1015.0838, 1013.4391, 1010.1612, 1010.2299, 1023.8341, 1011.8823, 1000.5451, 999.6638, 998.7858, 987.4057, 995.7576, 987.0883, 982.2467, 982.1057, 976.0313, 974.0374, 977.333, 985.7118, 985.4352, 989.2399, 994.7814, 994.221, 993.302, 992.1595, 974.4546, 974.6581, 978.4625, 988.9793, 987.0331, 984.4553, 983.4771, 974.0751, 973.6206, 965.2969, 967.768, 971.067, 971.0889, 965.449, 971.4746, 968.8909, 973.6257, 982.4815, 990.0573, 996.3394, 992.512, 983.0908, 974.1946, 965.6558, 959.3316, 956.7886, 943.7902, 948.5879, 963.673, 960.4892, 965.7126, 967.3483, 970.1678, 958.4293, 961.2468, 955.0185, 951.1211, 946.7029, 942.4834, 947.1571, 947.0276, 961.0235, 954.7518, 949.5334, 944.0529, 935.3442, 924.5955, 927.5307, 929.7883, 919.1908, 919.4065, 918.371, 921.3777, 929.4295, 942.7737, 949.4334, 954.6624, 962.7175, 956.2436, 950.7582, 958.6938, 956.3965, 944.6777, 956.2773, 954.0037, 932.8004, 928.7932, 940.8626, 926.6375, 937.3218, 958.6777, 964.0768, 952.0505, 950.5115, 948.4479, 939.3607, 935.7749, 770.943, 757.8363, 745.7519, 763.9224, 784.2032, 788.0772, 801.267, 815.823, 805.0596, 785.733, 791.6119, 783.7698, 777.2689, 776.5628, 783.5845, 799.8222, 792.5408, 800.3527, 795.215, 800.7316, 793.3939, 802.6706, 811.1229, 811.8536, 816.8235, 824.7065, 818.9618, 816.8607, 808.3944, 795.1955, 776.3799, 768.6471, 761.8164, 765.1176, 766.0153, 779.6374, 786.0757, 773.0423, 765.5792, 764.6384, 782.6235, 810.0114, 818.5287, 818.1307, 806.6282, 792.5812, 775.8072, 764.1966, 772.322, 768.9524, 772.9565, 761.8561, 757.6279, 747.5509, 751.0812, 746.2558, 759.4676, 776.8661, 781.1691, 784.5929, 798.0746, 794.1206, 780.7399, 778.5812, 792.9143, 813.2369, 814.7838, 819.026, 820.9111, 885.9339, 969.7707, 1003.4679, 1016.2571, 1030.1364, 972.6434, 912.3593, 905.3346, 930.5082, 928.4695, 956.5865, 989.6098, 1031.7731, 1069.3736, 1130.5232, 1170.21, 1199.5514, 1221.1681, 1237.044, 1256.4062, 1270.5876, 1286.3091, 1295.0264, 1313.8072, 1329.798, 1335.5542, 1336.4502, 1343.9396, 1339.8556, 1351.9698, 1362.3803, 1374.9839, 1362.5866, 1365.8521, 1366.4755, 1365.9438, 1361.402, 1368.8723, 1368.0455, 1370.7146, 1364.1305, 1377.6213, 1377.9086, 1368.8143, 1367.461, 1370.4544, 1371.3131, 1377.6303, 1377.8511, 1381.5674, 1379.2977, 1381.5947, 1361.9709, 1370.2396, 1374.2347, 1366.8974, 1355.6494, 1364.5881, 1349.1647, 1338.2199, 1345.7384, 1361.5699, 1367.6073, 1380.1442, 1395.4343, 1385.6615, 1372.8999, 1370.573, 1370.7435, 1360.4678, 1364.7567, 1359.9372, 1349.8359, 1351.2706, 1365.1294, 1364.0234, 1392.4452, 1401.6108, 1400.6695, 1386.2857, 1388.6885, 1377.9086, 1372.2517, 1371.9156, 1379.2685, 1385.6579, 1373.927, 1371.754, 1368.5747, 1367.0114, 1362.5276, 1363.9995, 1369.0521, 1368.8681, 1384.4765, 1378.8405, 1378.6151, 1382.5593, 1387.9971, 1371.732, 1376.1441, 1380.364, 1378.7928, 1368.9345, 1368.0102, 1376.1059, 1383.2482, 1378.0472, 1370.8553, 1370.4587, 1375.6658, 1366.2872, 1379.8434, 1390.603, 1395.1468, 1390.2211, 1384.3907, 1374.8753, 1377.2276, 1377.0982, 1371.5311, 1378.3185, 1385.8293, 1381.8495, 1367.9317, 1380.1028, 1379.8653, 1367.4705, 1351.2862, 1363.7038, 1349.9331, 1360.0009, 1361.8892, 1377.0465, 1382.0664, 1380.6035, 1380.4988, 1380.8638, 1372.7751, 1358.9864, 1364.3868, 1357.5578, 1357.5977, 1355.0698, 1358.8151, 1353.26, 1348.2656, 1347.9021, 1343.5709, 1343.5168, 1348.2359, 1354.7574, 1355.4504, 1367.8658, 1370.2762, 1373.4483, 1370.7956, 1363.3388, 1355.879, 1357.7108, 1360.4041, 1360.3489, 1371.9408, 1380.8992, 1377.2768, 1368.4755, 1363.3053, 1366.4712, 1363.8744, 1359.4646, 1360.1674, 1359.3972, 1352.4229, 1349.1245, 1359.1523, 1364.564, 1369.0246, 1374.0508, 1367.073, 1365.7079, 1357.6626, 1343.3174, 1335.4061, 1338.2838, 1334.7667, 1333.3732, 1328.9323, 1330.4992, 1335.2204, 1335.2879, 1339.3959, 1343.8426, 1345.2876, 1354.1573, 1349.1856, 1353.8924, 1359.5196, 1357.9431, 1346.0194, 1357.7789, 1357.2725, 1356.8047, 1353.7432, 1357.8778, 1357.279, 1354.2659, 1357.1141, 1360.9292, 1367.7369, 1372.6273, 1371.5874, 1367.3136, 1369.3025, 1360.0364, 1354.0252, 1353.3615, 1357.2325, 1345.4091, 1340.5862, 1342.8216, 1341.6592, 1334.9786, 1343.5994, 1349.0226, 1339.9034, 1339.8304, 1346.0953, 1339.9869, 1335.315, 1337.3079, 1336.5715, 1340.4444, 1354.2357, 1370.1224, 1378.7307, 1379.9508, 1365.8797, 1356.5842, 1350.4367, 1335.5538, 1344.1722, 1343.9125, 1351.5426, 1348.199, 1360.899, 1359.9313, 1366.1409, 1371.2208, 1379.0304, 1382.5708, 1384.7204, 1383.8992, 1379.1008, 1374.7729, 1359.9191, 1354.1172, 1351.8529, 1349.0833, 1341.5454, 1351.5255, 1350.431, 1356.8146, 1358.7656, 1365.3269, 1364.9431, 1366.631, 1360.4874, 1355.2152, 1358.7005, 1368.2299, 1363.7082, 1370.059, 1376.8738, 1366.2659, 1349.0049, 1366.3064, 1365.0829, 1376.7078, 1388.0081, 1404.1403, 1390.7425, 1395.0293, 1378.9801, 1367.3256, 1362.0991, 1364.6094, 1354.9031, 1353.7332, 1355.2581, 1342.6576, 1348.3014, 1356.1723, 1371.4283, 1367.0735, 1373.578, 1382.0668, 1372.1796, 1378.207, 1378.8784, 1388.2788, 1384.4918, 1366.59, 1362.3637, 1377.1559, 1385.4569, 1380.4629, 1388.1607, 1392.9759, 1383.4105, 1379.1023, 1386.516, 1399.0417, 1401.0841, 1394.361, 1386.4877, 1380.187, 1365.4649, 1353.3016, 1355.7279, 1357.9931, 1350.9256, 1355.7186, 1369.6754, 1373.8288, 1382.3468, 1395.4189, 1390.5029, 1386.8633, 1388.8058, 1381.5241, 1393.0759, 1375.6446, 1365.7619, 1358.0941, 1367.3018, 1350.1741, 1371.0745, 1385.838, 1385.603, 1369.2039, 1376.0055, 1384.2491, 1384.5682, 1396.8139, 1418.3658, 1416.778, 1410.158, 1408.9264, 1403.2051, 1391.3539, 1396.5059, 1388.0211, 1374.9082, 1369.2135, 1378.1975, 1364.6066, 1359.5322, 1354.2922, 1368.703, 1368.5208, 1373.6731, 1383.5117, 1400.3879, 1393.4152, 1378.1168, 1364.1499, 1359.5792, 1353.0026, 1357.5476, 1359.2456, 1369.5849, 1373.1943, 1372.4231, 1353.5337, 1340.596, 1335.708, 1328.0028, 1330.2829, 1336.832, 1343.4294, 1366.4506, 1375.6982, 1380.9705, 1384.9823, 1382.3623, 1369.2053, 1365.4555, 1360.5891, 1347.5065, 1354.8904, 1348.4929, 1348.3129, 1356.0756, 1350.554, 1339.8713, 1346.8866, 1355.5239, 1381.4118, 1382.4912, 1378.934, 1380.8874, 1368.8742, 1370.4616, 1365.5845, 1369.9524, 1349.8104, 1348.22, 1344.7895, 1350.4568, 1344.6752, 1353.1054, 1354.3282, 1352.2624, 1346.3093, 1347.3563, 1346.6762, 1336.3459, 1340.3673, 1348.9338, 1342.7555, 1344.3949, 1342.5281, 1333.1926, 1320.0272, 1323.2783, 1333.2061, 1338.7266, 1339.7386, 1347.0615, 1343.7917, 1332.6601, 1336.0047, 1335.4876, 1324.5035, 1321.8359, 1321.1697, 1319.0138, 1315.2431, 1312.5121, 1343.4203, 1338.8717, 1331.6775, 1322.9244, 1318.0362, 1310.816, 1321.7368, 1314.1894, 1314.4201, 1322.6703, 1330.1945, 1335.717, 1352.1948, 1357.4488, 1354.8395, 1350.7415, 1332.5945, 1330.2049, 1325.509, 1314.4949, 1305.2541, 1307.5292, 1314.601, 1305.8255, 1308.7152, 1308.4823, 1311.1685, 1318.4754, 1334.7864, 1338.8564, 1293.19, 1293.2667, 1292.7449, 1294.4715, 1306.9834, 1306.4165, 1294.5289, 1295.0713, 1290.9031, 1281.81, 1291.6441, 1302.0391, 1298.9444, 1306.18, 1313.127, 1308.9509, 1308.0279, 1307.2502, 1309.6825, 1317.9792, 1319.1761, 1314.9458, 1309.848, 1298.2314, 1277.2732, 1270.0605, 1273.5297, 1267.1429, 1266.8957, 1275.6693, 1285.4705, 1274.4266, 1278.2602, 1283.2233, 1292.7671, 1286.6295, 1281.6284, 1286.9236, 1292.2375, 1289.8843, 1282.707, 1299.9073, 1279.6421, 1284.4193, 1286.7726, 1291.7005, 1294.8435, 1296.9812,