通过 R 函数识别拟合平滑样条的所有局部极值 'smooth.spline'
Identify all local extrema of a fitted smoothing spline via R function 'smooth.spline'
我有一个二维数据集。
我使用 R 的 smooth.spline
函数按照本文中的示例平滑我的点图:
https://stat.ethz.ch/R-manual/R-devel/library/stats/html/predict.smooth.spline.html
这样我就得到了类似于这张图绿线的样条图
我想知道 X
值,其中平滑样条的一阶导数等于零(以确定确切的最小值或最大值)。
我的问题是我输入 predict()
函数的初始数据集(或我可以自动生成的数据集)不包含与平滑样条对应的精确 X
值极值。
如何找到这样的 X
值?
这是上面绿色样条线的一阶导数的图片
但极值的确切 X 坐标仍然不准确。
我生成图片的近似 R 脚本如下所示
sp1 <- smooth.spline(df)
pred.prime <- predict(sp1, deriv=1)
pred.second <- predict(sp1, deriv=2)
d1 <- data.frame(pred.prime)
d2 <- data.frame(pred.second)
dfMinimums <- d1[abs(d1$y) < 1e-4, c('x','y')]
我认为这里有两个问题。
- 您使用的是原始 x 值,它们之间的距离太远并且
- 由于 x 的间距较宽,您认为导数 "close enough" 为零的阈值过高。
这基本上是您的代码,但具有更多的 x 值并且需要更小的导数。由于您没有提供任何数据,我对其进行了粗略的近似,足以说明问题。
## Coarse approximation of your data
x = runif(300, 0,45000)
y = sin(x/5000) + sin(x/950)/4 + rnorm(300, 0,0.05)
df = data.frame(x,y)
sp1 <- smooth.spline(df)
样条代码
Sx = seq(0,45000,10)
pred.spline <- predict(sp1, Sx)
d0 <- data.frame(pred.spline)
pred.prime <- predict(sp1, Sx, deriv=1)
d1 <- data.frame(pred.prime)
Mins = which(abs(d1$y) < mean(abs(d1$y))/150)
plot(df, pch=20, col="navy")
lines(sp1, col="darkgreen")
points(d0[Mins,], pch=20, col="red")
极值看起来不错。
plot(d1, type="l")
points(d1[Mins,], pch=20, col="red")
识别出的点看起来像导数的零点。
我发现@G5W 实现中的一个警告是,它有时 returns 多个 记录围绕极值而不是单个极值收盘。在图表上看不到,因为它们实际上都落在一个点上。
来自 的以下代码片段用一阶导数的最小值过滤掉单个极值点:
library(tidyverse)
df2 <- df %>%
group_by(round(y, 4)) %>%
filter(abs(d1) == min(abs(d1))) %>%
ungroup() %>%
select(-5)
你可以使用我的R包SplinesUtils
:https://github.com/ZheyuanLi/SplinesUtils,可以通过
安装
devtools::install_github("ZheyuanLi/SplinesUtils")
要使用的函数是SmoothSplinesAsPiecePoly
和solve
。我将只使用文档下的示例。
library(SplinesUtils)
## a toy dataset
set.seed(0)
x <- 1:100 + runif(100, -0.1, 0.1)
y <- poly(x, 9) %*% rnorm(9)
y <- y + rnorm(length(y), 0, 0.2 * sd(y))
## fit a smoothing spline
sm <- smooth.spline(x, y)
## coerce "smooth.spline" object to "PiecePoly" object
oo <- SmoothSplineAsPiecePoly(sm)
## plot the spline
plot(oo)
## find all stationary / saddle points
xs <- solve(oo, deriv = 1)
#[1] 3.791103 15.957159 21.918534 23.034192 25.958486 39.799999 58.627431
#[8] 74.583000 87.049227 96.544430
## predict the "PiecePoly" at stationary / saddle points
ys <- predict(oo, xs)
#[1] -0.92224176 0.38751847 0.09951236 0.10764884 0.05960727 0.52068566
#[7] -0.51029209 0.15989592 -0.36464409 0.63471723
points(xs, ys, pch = 19)
我有一个二维数据集。
我使用 R 的 smooth.spline
函数按照本文中的示例平滑我的点图:
https://stat.ethz.ch/R-manual/R-devel/library/stats/html/predict.smooth.spline.html
这样我就得到了类似于这张图绿线的样条图
我想知道 X
值,其中平滑样条的一阶导数等于零(以确定确切的最小值或最大值)。
我的问题是我输入 predict()
函数的初始数据集(或我可以自动生成的数据集)不包含与平滑样条对应的精确 X
值极值。
如何找到这样的 X
值?
这是上面绿色样条线的一阶导数的图片
但极值的确切 X 坐标仍然不准确。
我生成图片的近似 R 脚本如下所示
sp1 <- smooth.spline(df)
pred.prime <- predict(sp1, deriv=1)
pred.second <- predict(sp1, deriv=2)
d1 <- data.frame(pred.prime)
d2 <- data.frame(pred.second)
dfMinimums <- d1[abs(d1$y) < 1e-4, c('x','y')]
我认为这里有两个问题。
- 您使用的是原始 x 值,它们之间的距离太远并且
- 由于 x 的间距较宽,您认为导数 "close enough" 为零的阈值过高。
这基本上是您的代码,但具有更多的 x 值并且需要更小的导数。由于您没有提供任何数据,我对其进行了粗略的近似,足以说明问题。
## Coarse approximation of your data
x = runif(300, 0,45000)
y = sin(x/5000) + sin(x/950)/4 + rnorm(300, 0,0.05)
df = data.frame(x,y)
sp1 <- smooth.spline(df)
样条代码
Sx = seq(0,45000,10)
pred.spline <- predict(sp1, Sx)
d0 <- data.frame(pred.spline)
pred.prime <- predict(sp1, Sx, deriv=1)
d1 <- data.frame(pred.prime)
Mins = which(abs(d1$y) < mean(abs(d1$y))/150)
plot(df, pch=20, col="navy")
lines(sp1, col="darkgreen")
points(d0[Mins,], pch=20, col="red")
极值看起来不错。
plot(d1, type="l")
points(d1[Mins,], pch=20, col="red")
识别出的点看起来像导数的零点。
我发现@G5W 实现中的一个警告是,它有时 returns 多个 记录围绕极值而不是单个极值收盘。在图表上看不到,因为它们实际上都落在一个点上。
来自
library(tidyverse)
df2 <- df %>%
group_by(round(y, 4)) %>%
filter(abs(d1) == min(abs(d1))) %>%
ungroup() %>%
select(-5)
你可以使用我的R包SplinesUtils
:https://github.com/ZheyuanLi/SplinesUtils,可以通过
devtools::install_github("ZheyuanLi/SplinesUtils")
要使用的函数是SmoothSplinesAsPiecePoly
和solve
。我将只使用文档下的示例。
library(SplinesUtils)
## a toy dataset
set.seed(0)
x <- 1:100 + runif(100, -0.1, 0.1)
y <- poly(x, 9) %*% rnorm(9)
y <- y + rnorm(length(y), 0, 0.2 * sd(y))
## fit a smoothing spline
sm <- smooth.spline(x, y)
## coerce "smooth.spline" object to "PiecePoly" object
oo <- SmoothSplineAsPiecePoly(sm)
## plot the spline
plot(oo)
## find all stationary / saddle points
xs <- solve(oo, deriv = 1)
#[1] 3.791103 15.957159 21.918534 23.034192 25.958486 39.799999 58.627431
#[8] 74.583000 87.049227 96.544430
## predict the "PiecePoly" at stationary / saddle points
ys <- predict(oo, xs)
#[1] -0.92224176 0.38751847 0.09951236 0.10764884 0.05960727 0.52068566
#[7] -0.51029209 0.15989592 -0.36464409 0.63471723
points(xs, ys, pch = 19)