tests/deriv-ex.R

In lpridge: Local Polynomial (Ridge) Regression

library(lpridge)

set.seed(17)
summary(x <- sort(round(250*runif(200))))
fn <- function(x) 10 - x/64 - (x/64 - 2)^2/2 + 5*exp(-(x-120)^2/200)
y <- fn(x) + round(sqrt(x+.1)/10 * rnorm(x), 2)
xx <- seq(0,250, by=1/4) # "fine grid"
plot(x,y)
lines(xx, fn(xx), col=adjustcolor("skyblue3", 0.5), lwd = 2.5)

.prt0 <- proc.time()
epa.15s	<- lpepa(x,y, bandw=15, x.out=xx, mnew=1)# "super-stable"
epa.20V   <- lpepa(x,y, bandw=20, x.out=xx,          var=TRUE)
cat('Time elapsed: ', (.pt <- proc.time()) - .prt0,'\n'); .prt0 <- .pt

with(epa.15s, lines(est ~ x.out, col="tomato"))
with(epa.20V, lines(est ~ x.out, col=3))

.d1 <- deriv(body(fn), "x", func=TRUE)
d1f <- function(x) drop(attr(.d1(x), "gradient"))
## manually:
df1 <- function(x) {c <- 1/64; c*(1-c*x) - (x-120)/20 *exp(-(x-120)^2/200)}
stopifnot(all.equal(d1f(xx), df1(xx)))

.prt0 <- proc.time()
epa.40Vd  <- lpepa(x,y,   bandw=40, x.out=xx, deriv=1, var=TRUE)
epa.40Vd2 <- lpepa(x,y,   bandw=40, x.out=xx, deriv=2, var=TRUE)
lpr.30Vd <-  lpridge(x,y, bandw=30, x.out=xx, deriv=1, var=TRUE)
cat('Time elapsed: ', (.pt <- proc.time()) - .prt0,'\n'); .prt0 <- .pt

curve(df1, 0, 250, n = 1001, lwd=2, col=2, main = "First derivative: true and estim.")
with(epa.40Vd, lines(est ~ x.out))
with(lpr.30Vd, lines(est ~ x.out, col=3))

parts <- c("est","est.var")
sub.s <- function(x, every.k= 10) x[seq(from=1, to=length(x), by=every.k)]
if(FALSE)
    dput(lapply(lapply(lpr.30Vd[parts], sub.s), signif, 7))

lpr.30.T <- list(
    est = c(0.02229422, 0.02222213, 0.02252096, 0.02342268,
            0.02441539, 0.02143815, 0.01499614, 0.007546031, 0.003249612,
            0.001089905, 0.0001838037, -0.0001373865, -0.0007027075, 0.001103937,
            0.00475638, 0.0065422, 0.008629815, 0.01287751, 0.01544041, 0.0175178,
            0.0190212, 0.01683765, 0.01229409, 0.007850588, 0.005160776,
            -0.0002930656, -0.004053807, -0.008135056, -0.01169963, -0.00867799,
            -0.003745031, 0.0006846396, 0.007440006, 0.01597122, 0.02660977,
            0.04034368, 0.05737227, 0.07538167, 0.09257077, 0.107556, 0.103767,
            0.09828084, 0.09084512, 0.07977893, 0.06323039, 0.05239796, 0.03386389,
            0.01385851, -0.01033999, -0.03081899, -0.0513893, -0.06863713,
            -0.08249654, -0.09545513, -0.1078559, -0.1203647, -0.1247685,
            -0.1146073, -0.1030085, -0.09261767, -0.0841685, -0.07661196,
            -0.07033584, -0.06093106, -0.0526194, -0.04492832, -0.03801281,
            -0.03376433, -0.0290251, -0.02622674, -0.02297068, -0.01834848,
            -0.01458384, -0.01296454, -0.01243847, -0.01467232, -0.01861357,
            -0.02130569, -0.02286563, -0.02498196, -0.02740975, -0.02819129,
            -0.02755911, -0.02649983, -0.02946556, -0.03115722, -0.02943064,
            -0.03391023, -0.0394636, -0.04549067, -0.05155681, -0.05801209,
            -0.06759929, -0.08252602, -0.1035041, -0.1188784, -0.1267757,
            -0.1267948, -0.1221933, -0.1226147, -0.1249049),
    est.var = c(0.001274448,
                0.001289882, 0.001278549, 0.001128746, 0.0008462655, 0.0005655128,
                0.000374932, 0.0002471065, 0.0001692552, 0.0001230342, 9.697478e-05,
                8.4373e-05, 8.859465e-05, 8.807464e-05, 8.70996e-05, 9.675577e-05,
                0.0001097343, 0.00011335, 0.0001190441, 0.0001365236, 0.0001703944,
                0.0001848792, 0.0001441367, 0.0001259812, 0.000117191, 0.0001040061,
                9.135686e-05, 9.000982e-05, 9.822586e-05, 9.39134e-05, 9.38495e-05,
                9.756342e-05, 9.374591e-05, 9.367529e-05, 0.0001056997, 0.0001237378,
                0.0001434425, 0.0001625226, 0.0001900015, 0.0001955096, 0.0001636594,
                0.0001430319, 0.0001173526, 0.0001066161, 9.859755e-05, 0.0001044603,
                9.404138e-05, 9.419802e-05, 0.0001002991, 0.0001048351, 9.715271e-05,
                9.338953e-05, 9.671686e-05, 9.790364e-05, 0.0001028107, 0.0001203195,
                0.000142119, 0.000129301, 0.0001129553, 0.0001048131, 9.196553e-05,
                8.226805e-05, 7.754489e-05, 7.527381e-05, 7.276773e-05, 7.525802e-05,
                8.658272e-05, 9.462503e-05, 9.939682e-05, 8.960028e-05, 8.113668e-05,
                7.477402e-05, 7.719303e-05, 7.779136e-05, 7.973575e-05, 8.443209e-05,
                9.440942e-05, 9.316252e-05, 9.004642e-05, 8.918413e-05, 0.00010154,
                0.0001283644, 0.0001627717, 0.0001599295, 0.0001271949, 0.0001121226,
                0.0001438388, 0.000117949, 9.915013e-05, 9.992847e-05, 0.0001223734,
                0.0001696298, 0.0002550821, 0.0004011194, 0.0006473775, 0.0009639706,
                0.00134593, 0.00156556, 0.001647994, 0.001672598, 0.002068033))


if(FALSE)
    dput(lapply(lapply(epa.40Vd2[parts], sub.s), signif, 7))

epa.40d2.T <- list(
    est = c(0.005973069, 0.00525159, 0.004265993, 5.543538e-05,
            -0.003643127, -0.004129886, -0.003571208, -0.002730489, -0.001991292,
            -0.001416728, -0.000786919, -0.0001962287, 0.0002237254, 0.0005535253,
            0.000638903, 0.000656796, 0.0006444588, 0.0003512309, 2.432636e-05,
            -0.0001374263, -0.0002169834, -0.000492684, -0.0006789096, -0.0008669604,
            -0.001060009, -0.0009372228, -0.0006881592, -0.0003489099, -3.460457e-05,
            0.0006400346, 0.001478762, 0.00234706, 0.003252709, 0.00419436,
            0.005028455, 0.005317064, 0.004104386, 0.003320867, 0.002594748,
            0.001372943, -9.628925e-05, -0.001334469, -0.00237591, -0.003761749,
            -0.005209561, -0.005916782, -0.006284806, -0.006347289, -0.006156965,
            -0.006377393, -0.005951912, -0.005491826, -0.004690328, -0.003953238,
            -0.002743524, -0.00162772, -0.0006884371, 0.0004363056, 0.001531293,
            0.002810299, 0.003789386, 0.003766839, 0.003539994, 0.003116755,
            0.0026953, 0.002353148, 0.002187459, 0.002055807, 0.001821656,
            0.001512261, 0.001216675, 0.0007704807, 0.0002910974, 0.0001802986,
            0.0001580854, -2.192169e-05, -0.0002479372, -0.0004419121, -0.0004892296,
            -0.0004730236, -0.000783885, -0.0008715044, -0.0008964609, -0.0008418996,
            -0.0008380501, -0.0009902708, -0.001255034, -0.001705962, -0.002377678,
            -0.003131305, -0.003915043, -0.004967261, -0.006274099, -0.007962944,
            -0.01046467, -0.01199683, -0.01286726, -0.01096378, -0.003806234,
            -0.0004184816, 0.005782601),
    est.var = c(0.0003689218, 0.0001914982,
                0.00012878, 6.976101e-05, 3.604109e-05, 2.104798e-05, 1.25106e-05,
                7.56681e-06, 4.684791e-06, 2.89014e-06, 1.84113e-06, 1.213595e-06,
                8.457589e-07, 6.487753e-07, 5.480269e-07, 5.501742e-07, 6.657365e-07,
                6.639384e-07, 4.103037e-07, 4.157615e-07, 4.437498e-07, 4.112107e-07,
                3.861566e-07, 4.242876e-07, 5.589571e-07, 5.894606e-07, 5.856352e-07,
                5.464175e-07, 5.274938e-07, 4.729866e-07, 4.722855e-07, 5.39823e-07,
                5.876371e-07, 5.931668e-07, 6.082821e-07, 6.101139e-07, 4.279594e-07,
                3.878046e-07, 4.104848e-07, 4.387643e-07, 4.491503e-07, 4.480333e-07,
                5.764631e-07, 6.153242e-07, 6.464596e-07, 6.331725e-07, 5.106045e-07,
                4.363449e-07, 3.999944e-07, 4.706624e-07, 4.28506e-07, 4.267681e-07,
                4.530885e-07, 4.858428e-07, 4.359352e-07, 3.847176e-07, 3.667357e-07,
                3.556107e-07, 3.635952e-07, 4.139824e-07, 5.483039e-07, 5.35252e-07,
                5.254244e-07, 5.197202e-07, 4.766399e-07, 3.929885e-07, 3.366216e-07,
                3.482119e-07, 3.416057e-07, 3.269638e-07, 3.293217e-07, 3.51608e-07,
                4.155127e-07, 3.886885e-07, 3.720152e-07, 3.752141e-07, 4.468934e-07,
                5.795008e-07, 7.095282e-07, 6.032648e-07, 4.541418e-07, 3.979016e-07,
                3.911268e-07, 3.673768e-07, 3.861204e-07, 4.654381e-07, 6.571598e-07,
                1.004496e-06, 1.572744e-06, 2.429328e-06, 3.880271e-06, 6.531213e-06,
                1.073932e-05, 1.69096e-05, 2.770833e-05, 4.302006e-05, 7.009248e-05,
                0.0001093674, 0.0001932428, 0.0002821257, 0.0005402126))

if(FALSE)
    dput(lapply(lapply(epa.40Vd[parts], sub.s), signif, 7))
epa.40.T <- list(
    est = c(0.04011727, 0.04654933, 0.0440402, 0.03461681, 0.02169244,
            0.01382214, 0.008939529, 0.006007981, 0.004316259, 0.003164288,
            0.003315283, 0.004228259, 0.005430772, 0.006873652, 0.007711286,
            0.008365776, 0.008919946, 0.008747495, 0.006356444, 0.0051292,
            0.004157813, 0.003257316, 0.003234971, 0.003143546, 0.002770634,
            0.002868093, 0.002691931, 0.002539738, 0.003842211, 0.006443181,
            0.01012377, 0.01443586, 0.02120269, 0.03131376, 0.04306628, 0.05337821,
            0.05567451, 0.05655216, 0.05845041, 0.06050494, 0.06024862, 0.05814509,
            0.05712463, 0.0506945, 0.04081804, 0.02885176, 0.01497677, 0.001749387,
            -0.009533884, -0.01921083, -0.0308403, -0.04205738, -0.05342632,
            -0.06363766, -0.07176413, -0.07835563, -0.08485778, -0.08959265,
            -0.09252134, -0.09298555, -0.08921572, -0.08009209, -0.07050485,
            -0.06215323, -0.0546681, -0.0487828, -0.04355158, -0.03819087,
            -0.0337525, -0.03083543, -0.02874027, -0.02767962, -0.02600663,
            -0.02462118, -0.02326894, -0.02216083, -0.02096434, -0.02040704,
            -0.0204265, -0.02058528, -0.02244195, -0.02486518, -0.02796743,
            -0.03073023, -0.0344798, -0.03773674, -0.04092551, -0.04379927,
            -0.04590977, -0.04882758, -0.05365846, -0.06037438, -0.06743245,
            -0.0748644, -0.08354326, -0.09749254, -0.1190881, -0.1504599,
            -0.1911819, -0.2135679, -0.2285541),
    est.var = c(0.004105202,
                0.002381893, 0.00167919, 0.001156419, 0.0007539224, 0.0005081698,
                0.0003436571, 0.0002337784, 0.0001611879, 0.0001122577, 8.070577e-05,
                6.11617e-05, 4.99236e-05, 4.489221e-05, 4.389396e-05, 4.723976e-05,
                5.534353e-05, 6.230872e-05, 5.480998e-05, 5.630862e-05, 5.690259e-05,
                5.469279e-05, 5.206698e-05, 5.231974e-05, 5.616829e-05, 5.437301e-05,
                5.071612e-05, 4.697193e-05, 4.425646e-05, 4.043741e-05, 3.928156e-05,
                4.143786e-05, 4.421631e-05, 4.775833e-05, 5.289455e-05, 5.878149e-05,
                5.199678e-05, 4.888727e-05, 4.994033e-05, 5.149355e-05, 5.304381e-05,
                5.257264e-05, 5.641748e-05, 5.249267e-05, 4.828498e-05, 4.439159e-05,
                3.929277e-05, 3.698952e-05, 3.693465e-05, 4.0169e-05, 3.904751e-05,
                3.98237e-05, 4.179261e-05, 4.379268e-05, 4.333012e-05, 4.323019e-05,
                4.46398e-05, 4.566203e-05, 4.734766e-05, 5.028663e-05, 5.154964e-05,
                4.678282e-05, 4.339737e-05, 4.226519e-05, 4.142497e-05, 3.91447e-05,
                3.66799e-05, 3.57632e-05, 3.461693e-05, 3.332138e-05, 3.305036e-05,
                3.378918e-05, 3.622651e-05, 3.58515e-05, 3.53884e-05, 3.605553e-05,
                3.964627e-05, 4.52568e-05, 5.235753e-05, 5.54541e-05, 5.017117e-05,
                4.323935e-05, 4.004631e-05, 3.838974e-05, 3.839393e-05, 4.07111e-05,
                4.621292e-05, 5.521834e-05, 6.754224e-05, 8.556588e-05, 0.0001136882,
                0.0001639589, 0.0002433641, 0.0003603769, 0.0005448669, 0.0008158418,
                0.001290934, 0.002045996, 0.003437501, 0.00526647, 0.009107401
                ))

stopifnot(TRUE,
          all.equal(lapply(lpr.30Vd[parts], sub.s), lpr.30.T,  tol = 1e-6)# Lnx.64: {1.23,1.395}e-7
        , all.equal(lapply(epa.40Vd [parts], sub.s), epa.40.T,   tol=6e-7)#  "      {7.42.7.97}e-8
        , all.equal(lapply(epa.40Vd2[parts], sub.s), epa.40d2.T, tol=7e-7)#  "      {8.70,9.26}e-8
)