diff -Nru sm-2.2-5.7/data/geyser.tab sm-2.2-5.7.1/data/geyser.tab --- sm-2.2-5.7/data/geyser.tab 2021-01-01 17:32:00.000000000 +0000 +++ sm-2.2-5.7.1/data/geyser.tab 1970-01-01 00:00:00.000000000 +0000 @@ -1,300 +0,0 @@ -waiting duration - 80.0000000 4.0166667 - 71.0000000 2.1500000 - 57.0000000 4.0000000 - 80.0000000 4.0000000 - 75.0000000 4.0000000 - 77.0000000 2.0000000 - 60.0000000 4.3833333 - 86.0000000 4.2833333 - 77.0000000 2.0333333 - 56.0000000 4.8333333 - 81.0000000 1.8333333 - 50.0000000 5.4500000 - 89.0000000 1.6166667 - 54.0000000 4.8666667 - 90.0000000 4.3833333 - 73.0000000 1.7666667 - 60.0000000 4.6666667 - 83.0000000 2.0000000 - 65.0000000 4.7333333 - 82.0000000 4.2166667 - 84.0000000 1.9000000 - 54.0000000 4.9666667 - 85.0000000 2.0000000 - 58.0000000 4.0000000 - 79.0000000 2.0000000 - 57.0000000 4.0000000 - 88.0000000 2.8333333 - 68.0000000 4.5000000 - 76.0000000 4.0666667 - 78.0000000 3.7166667 - 74.0000000 3.5166667 - 85.0000000 4.4666667 - 75.0000000 2.2166667 - 65.0000000 4.8833333 - 76.0000000 2.6000000 - 58.0000000 4.1500000 - 91.0000000 2.2000000 - 50.0000000 4.7666667 - 87.0000000 1.8333333 - 48.0000000 4.6000000 - 93.0000000 2.2666667 - 54.0000000 4.1333333 - 86.0000000 2.0000000 - 53.0000000 4.0000000 - 78.0000000 2.0000000 - 52.0000000 4.0000000 - 83.0000000 1.8833333 - 60.0000000 4.2666667 - 87.0000000 2.0833333 - 49.0000000 4.4666667 - 80.0000000 2.5000000 - 60.0000000 4.0000000 - 92.0000000 1.7666667 - 43.0000000 4.3333333 - 89.0000000 2.1833333 - 60.0000000 4.4833333 - 84.0000000 3.8833333 - 69.0000000 3.3333333 - 74.0000000 3.7333333 - 71.0000000 4.0000000 -108.0000000 1.9500000 - 50.0000000 5.2666667 - 77.0000000 2.0000000 - 57.0000000 4.0000000 - 80.0000000 2.0000000 - 61.0000000 4.0000000 - 82.0000000 2.0000000 - 48.0000000 4.0000000 - 81.0000000 3.5333333 - 73.0000000 2.1666667 - 62.0000000 4.5000000 - 79.0000000 2.0166667 - 54.0000000 4.1500000 - 80.0000000 4.2000000 - 73.0000000 4.3333333 - 81.0000000 1.9333333 - 62.0000000 4.6500000 - 81.0000000 3.8166667 - 71.0000000 4.0333333 - 79.0000000 4.1666667 - 81.0000000 4.6666667 - 74.0000000 1.8166667 - 59.0000000 4.0000000 - 81.0000000 3.0000000 - 66.0000000 4.0000000 - 87.0000000 2.0000000 - 53.0000000 4.4500000 - 80.0000000 2.0500000 - 50.0000000 4.2500000 - 87.0000000 1.9166667 - 51.0000000 4.6666667 - 82.0000000 1.7333333 - 58.0000000 4.3833333 - 81.0000000 1.7666667 - 49.0000000 4.6000000 - 92.0000000 1.8666667 - 50.0000000 4.4500000 - 88.0000000 1.6333333 - 62.0000000 5.0333333 - 93.0000000 1.8166667 - 56.0000000 5.1000000 - 89.0000000 1.6333333 - 51.0000000 4.2833333 - 79.0000000 2.0000000 - 58.0000000 4.0000000 - 82.0000000 2.0000000 - 52.0000000 4.5333333 - 88.0000000 2.0000000 - 52.0000000 4.0000000 - 78.0000000 2.9333333 - 69.0000000 4.7333333 - 75.0000000 3.9000000 - 77.0000000 1.9500000 - 53.0000000 4.1166667 - 80.0000000 1.8000000 - 55.0000000 4.6666667 - 87.0000000 1.8333333 - 53.0000000 4.7000000 - 85.0000000 2.1166667 - 61.0000000 4.7833333 - 93.0000000 1.8166667 - 54.0000000 4.1000000 - 76.0000000 4.6500000 - 80.0000000 4.0000000 - 81.0000000 2.0000000 - 59.0000000 4.0000000 - 86.0000000 4.0000000 - 78.0000000 4.2166667 - 71.0000000 4.1333333 - 77.0000000 3.9333333 - 76.0000000 3.7500000 - 94.0000000 4.4166667 - 75.0000000 2.4666667 - 50.0000000 4.1666667 - 83.0000000 3.8000000 - 82.0000000 4.3166667 - 72.0000000 3.8666667 - 77.0000000 4.6833333 - 75.0000000 1.7000000 - 65.0000000 4.9666667 - 79.0000000 4.2666667 - 72.0000000 4.5833333 - 78.0000000 4.0000000 - 77.0000000 4.0000000 - 79.0000000 4.0000000 - 75.0000000 4.0000000 - 78.0000000 1.9833333 - 64.0000000 4.6000000 - 80.0000000 0.8333333 - 49.0000000 4.9166667 - 88.0000000 1.7333333 - 54.0000000 4.5833333 - 85.0000000 1.7000000 - 51.0000000 4.7500000 - 96.0000000 1.8333333 - 50.0000000 4.5000000 - 80.0000000 1.8666667 - 78.0000000 4.4500000 - 81.0000000 4.4500000 - 72.0000000 4.0000000 - 75.0000000 4.8000000 - 78.0000000 4.0000000 - 87.0000000 4.0000000 - 69.0000000 2.0000000 - 55.0000000 4.0000000 - 83.0000000 1.9333333 - 49.0000000 4.5833333 - 82.0000000 2.0000000 - 57.0000000 3.7000000 - 84.0000000 2.8666667 - 57.0000000 4.8333333 - 84.0000000 3.4500000 - 73.0000000 4.3833333 - 78.0000000 1.8000000 - 57.0000000 4.4000000 - 79.0000000 2.4833333 - 57.0000000 4.5166667 - 90.0000000 2.1000000 - 62.0000000 4.3500000 - 87.0000000 4.3666667 - 78.0000000 1.7833333 - 52.0000000 4.9166667 - 98.0000000 1.8166667 - 48.0000000 4.0000000 - 78.0000000 4.0000000 - 79.0000000 4.0000000 - 65.0000000 3.8666667 - 84.0000000 1.8500000 - 50.0000000 4.7000000 - 83.0000000 2.0166667 - 60.0000000 4.4666667 - 80.0000000 1.8666667 - 50.0000000 4.1666667 - 88.0000000 1.9000000 - 50.0000000 4.2500000 - 84.0000000 3.2500000 - 74.0000000 4.2166667 - 76.0000000 1.8833333 - 65.0000000 4.9833333 - 89.0000000 1.8500000 - 49.0000000 4.0000000 - 88.0000000 1.9666667 - 51.0000000 4.7666667 - 78.0000000 4.0000000 - 85.0000000 2.0000000 - 65.0000000 4.0000000 - 75.0000000 4.0000000 - 77.0000000 2.3833333 - 69.0000000 4.4166667 - 92.0000000 4.2166667 - 68.0000000 4.3666667 - 87.0000000 2.0000000 - 61.0000000 4.4500000 - 81.0000000 1.7500000 - 55.0000000 4.5000000 - 93.0000000 1.6166667 - 53.0000000 4.7000000 - 84.0000000 2.5666667 - 70.0000000 3.7000000 - 73.0000000 4.2333333 - 93.0000000 1.9333333 - 50.0000000 4.3500000 - 87.0000000 4.0000000 - 77.0000000 4.0000000 - 74.0000000 4.0000000 - 72.0000000 4.2166667 - 82.0000000 4.0000000 - 74.0000000 4.1333333 - 80.0000000 1.8833333 - 49.0000000 4.4666667 - 91.0000000 1.9500000 - 53.0000000 4.2166667 - 86.0000000 1.7166667 - 49.0000000 4.4500000 - 79.0000000 4.2500000 - 89.0000000 3.9666667 - 87.0000000 4.3833333 - 76.0000000 1.9666667 - 59.0000000 4.4500000 - 80.0000000 4.2666667 - 89.0000000 1.9166667 - 45.0000000 4.4166667 - 93.0000000 3.0000000 - 72.0000000 4.0000000 - 71.0000000 2.0000000 - 54.0000000 4.0000000 - 79.0000000 3.2833333 - 74.0000000 1.8333333 - 65.0000000 4.6166667 - 78.0000000 1.8333333 - 57.0000000 4.6166667 - 87.0000000 4.6000000 - 72.0000000 4.2500000 - 84.0000000 1.9333333 - 47.0000000 4.9833333 - 84.0000000 1.9666667 - 57.0000000 4.3000000 - 87.0000000 4.2000000 - 68.0000000 4.5333333 - 86.0000000 4.4000000 - 75.0000000 4.6166667 - 73.0000000 2.0000000 - 53.0000000 4.0000000 - 82.0000000 4.0000000 - 93.0000000 3.9166667 - 77.0000000 2.0000000 - 54.0000000 4.5000000 - 96.0000000 1.8000000 - 48.0000000 4.0000000 - 89.0000000 2.7500000 - 63.0000000 4.7333333 - 84.0000000 3.9666667 - 76.0000000 1.9500000 - 62.0000000 4.9666667 - 83.0000000 1.8500000 - 50.0000000 4.8000000 - 85.0000000 4.0000000 - 78.0000000 4.0000000 - 78.0000000 4.0000000 - 81.0000000 4.0000000 - 78.0000000 4.0000000 - 76.0000000 4.0000000 - 74.0000000 4.0000000 - 81.0000000 2.0000000 - 66.0000000 4.0000000 - 84.0000000 1.9333333 - 48.0000000 4.3333333 - 93.0000000 1.6666667 - 47.0000000 4.7666667 - 87.0000000 1.9500000 - 51.0000000 4.6833333 - 78.0000000 1.9333333 - 54.0000000 4.4166667 - 87.0000000 2.1333333 - 52.0000000 4.0833333 - 85.0000000 2.0666667 - 58.0000000 4.0000000 - 88.0000000 4.0000000 - 79.0000000 2.0000000 Binary files /tmp/tmpvb7kb_s4/FeQaMLFFqF/sm-2.2-5.7/data/geyser.tab.gz and /tmp/tmpvb7kb_s4/M4JMIwRDmR/sm-2.2-5.7.1/data/geyser.tab.gz differ diff -Nru sm-2.2-5.7/debian/changelog sm-2.2-5.7.1/debian/changelog --- sm-2.2-5.7/debian/changelog 2021-09-20 02:48:51.000000000 +0000 +++ sm-2.2-5.7.1/debian/changelog 2022-07-08 00:24:22.000000000 +0000 @@ -1,3 +1,12 @@ +sm (2.2-5.7.1-1) unstable; urgency=medium + + * New upstream release + + * debian/control: Set Build-Depends: to current R version + * debian/control: Set Standards-Version: to current version + + -- Dirk Eddelbuettel Thu, 07 Jul 2022 19:24:22 -0500 + sm (2.2-5.7-1) unstable; urgency=medium * New upstream release diff -Nru sm-2.2-5.7/debian/control sm-2.2-5.7.1/debian/control --- sm-2.2-5.7/debian/control 2021-09-20 02:48:06.000000000 +0000 +++ sm-2.2-5.7.1/debian/control 2022-07-08 00:24:22.000000000 +0000 @@ -2,8 +2,8 @@ Section: gnu-r Priority: optional Maintainer: Dirk Eddelbuettel -Build-Depends: debhelper-compat (= 11), dh-r, r-base-dev (>= 4.1.1) -Standards-Version: 4.6.0 +Build-Depends: debhelper-compat (= 11), dh-r, r-base-dev (>= 4.2.1) +Standards-Version: 4.6.1 Vcs-Browser: https://salsa.debian.org/edd/r-cran-sm Vcs-Git: https://salsa.debian.org/edd/r-cran-sm.git Homepage: https://cran.r-project.org/package=sm diff -Nru sm-2.2-5.7/DESCRIPTION sm-2.2-5.7.1/DESCRIPTION --- sm-2.2-5.7/DESCRIPTION 2021-09-13 12:10:02.000000000 +0000 +++ sm-2.2-5.7.1/DESCRIPTION 2022-07-04 11:29:42.000000000 +0000 @@ -2,7 +2,7 @@ Type: Package Title: Smoothing Methods for Nonparametric Regression and Density Estimation -Version: 2.2-5.7 +Version: 2.2-5.7.1 Date: 2021-09-13 Author: Adrian Bowman and Adelchi Azzalini. Ported to R by B. D. Ripley up to version 2.0, @@ -10,13 +10,13 @@ version 2.2 by Adrian Bowman. Maintainer: Adrian Bowman Depends: R (>= 3.1.0) -Suggests: rgl, misc3d, akima, gam, tkrplot, rpanel (>= 1.1-4), tcltk +Suggests: rgl, misc3d, interp, gam, tkrplot, rpanel (>= 1.1-4), tcltk Description: This is software linked to the book 'Applied Smoothing Techniques for Data Analysis - The Kernel Approach with S-Plus Illustrations' Oxford University Press. License: GPL (>= 2) LazyData: TRUE NeedsCompilation: yes -Packaged: 2021-09-13 11:48:34 UTC; adrianbowman +Packaged: 2022-07-04 08:59:51 UTC; ripley Repository: CRAN -Date/Publication: 2021-09-13 12:10:02 UTC +Date/Publication: 2022-07-04 11:29:42 UTC diff -Nru sm-2.2-5.7/man/sm.ancova.Rd sm-2.2-5.7.1/man/sm.ancova.Rd --- sm-2.2-5.7/man/sm.ancova.Rd 2021-09-09 20:52:48.000000000 +0000 +++ sm-2.2-5.7.1/man/sm.ancova.Rd 2022-07-04 08:53:37.000000000 +0000 @@ -75,8 +75,8 @@ a plot on the current graphical device is produced, unless \code{display="none"} } \details{ -see Sections 6.4 and 6.5 of the book by Bowman \& Azzalini, and -the papers by Young \& Bowman listed below. +see Sections 6.4 and 6.5 of the book by Bowman & Azzalini, and +the papers by Young & Bowman listed below. This function is a developed version of code originally written by Stuart Young. } \references{ diff -Nru sm-2.2-5.7/man/sm.discontinuity.Rd sm-2.2-5.7.1/man/sm.discontinuity.Rd --- sm-2.2-5.7/man/sm.discontinuity.Rd 2021-09-09 21:02:01.000000000 +0000 +++ sm-2.2-5.7.1/man/sm.discontinuity.Rd 2022-07-04 08:54:04.000000000 +0000 @@ -96,7 +96,7 @@ Bowman, A.W., Pope, A. and Ismail, B. (2006). Detecting discontinuities in nonparametric regression curves and surfaces. - \emph{Statistics \& Computing}, 16, 377--390. + \emph{Statistics & Computing}, 16, 377--390. } \seealso{ \code{\link{sm.regression}}, \code{\link{sm.options}} diff -Nru sm-2.2-5.7/man/sm.sigma2.compare.Rd sm-2.2-5.7.1/man/sm.sigma2.compare.Rd --- sm-2.2-5.7/man/sm.sigma2.compare.Rd 2021-01-01 17:32:00.000000000 +0000 +++ sm-2.2-5.7.1/man/sm.sigma2.compare.Rd 2022-07-04 08:54:26.000000000 +0000 @@ -24,10 +24,10 @@ \details{see the reference below.} -\references{Bock, M., Bowman, A.W.\ \& Ismail, B. (2007). +\references{Bock, M., Bowman, A.W.\ & Ismail, B. (2007). Estimation and inference for error variance in bivariate nonparametric regression. - \emph{Statistics \& Computing}, to appear.} + \emph{Statistics & Computing}, to appear.} \seealso{\code{\link{sm.sigma}}} diff -Nru sm-2.2-5.7/man/sm.sigma.Rd sm-2.2-5.7.1/man/sm.sigma.Rd --- sm-2.2-5.7/man/sm.sigma.Rd 2021-01-01 17:32:00.000000000 +0000 +++ sm-2.2-5.7.1/man/sm.sigma.Rd 2022-07-04 08:54:19.000000000 +0000 @@ -35,10 +35,10 @@ a p-value for the test of constant variance.} \section{Side Effects}{none.} \details{see the reference below.} -\references{Bock, M., Bowman, A.W.\ \& Ismail, B. (2007). +\references{Bock, M., Bowman, A.W.\ & Ismail, B. (2007). Estimation and inference for error variance in bivariate nonparametric regression. - \emph{Statistics \& Computing}, to appear.} + \emph{Statistics & Computing}, to appear.} \seealso{\code{\link{sm.sigma2.compare}}} \examples{ \dontrun{ diff -Nru sm-2.2-5.7/MD5 sm-2.2-5.7.1/MD5 --- sm-2.2-5.7/MD5 2021-09-13 12:10:02.000000000 +0000 +++ sm-2.2-5.7.1/MD5 2022-07-04 11:29:42.000000000 +0000 @@ -1,5 +1,5 @@ 7555e62b2f78e21c4fe2421de9c904e9 *ChangeLog -c323a85b1fdc8cf5e1e1b82cabd2014c *DESCRIPTION +545dda96cdfe5ca4960c83d3329de1f2 *DESCRIPTION 36128133045efeb0fabc4abf4b84d9ac *NAMESPACE 5b09398cb5221935c79a47bde6d8d962 *NEWS.md b1c0be3283f34d3b8df4af7bb2fddd63 *R/ancova.r @@ -12,11 +12,11 @@ 467f10cc62682cd119a2568954f4c1d0 *R/ps-normal.r 063cd395a89edfde65dfb7a24b0ad26e *R/regression.r 5f430339b0ab9ff420bc5b973412ad3f *R/rpanel.r -55746c572357257a5a3345b24ad406fb *R/sm.r +7451131c0b88fe2619e2e37c11b66766 *R/sm.r 3fc95761d3a67b6d8bdbbcc777bf8508 *R/sphere.r f044ed7d790dd89af28c8ef9abf83457 *R/survival.r 5a4838807f8064f5d8826e3a2696b017 *R/utilities.r -a40a42f7cd644a2683dcc482f84dc393 *R/variogram.r +4a507919071cf82f62ccecd54232267f *R/variogram.r e9ff1bce747db2a48fe7b64a346fda87 *R/zzz.R 2adb64c8e856f6367c806993979144ac *data/aircraft.rda ca1386a986a51ab9cb8abb4363fb2673 *data/airpc.rda @@ -29,7 +29,7 @@ 56e642fe3f1704f579daf5ded4b47118 *data/dogs.rda 72dbe71a11191af4d4319d07cb90eeb0 *data/follicle.rda 0a4394287ff1a533185d030297544ce2 *data/geys3d.rda -0b247c4faca947d14f00e69aa78e5039 *data/geyser.tab +98b422bbac933b73998d775982c594a2 *data/geyser.tab.gz 15de7a46319ec1554e24287173a47cfb *data/lcancer.rda 9936554a9782c5c507917d76fc427d2c *data/mackerel.rda 2b4573cf5a89636bb31400f291ebc507 *data/magrem.rda @@ -214,13 +214,13 @@ 4a4994a4365bc098acefcb9daef8616d *man/sig.trace.Rd 479235458a5a73dd5b9c8e51d8b511d7 *man/sm-internal.Rd 5ae1a7393157440efd1ded51149f20c4 *man/sm.Rd -03dd91ba0b6f1f3f58970025dfe72d41 *man/sm.ancova.Rd +267442310c996547666dfd25ff83d741 *man/sm.ancova.Rd ea20163062c9e452bea4b3b7934b51b7 *man/sm.autoregression.Rd 6ebe6595a79c5bb77f612c518e4461a6 *man/sm.binomial.Rd 0697eb7b8e86266b5ad7b34497424660 *man/sm.binomial.bootstrap.Rd 54e29fd661171dd70a6db6904a720391 *man/sm.density.Rd 9440022c6f7985ca5989bcba3a454710 *man/sm.density.compare.Rd -50e314f14e9496ae6ad6c1778d18014a *man/sm.discontinuity.Rd +6c6919b95354b9f38058b493f6c3872f *man/sm.discontinuity.Rd 250add6d8ea3fc7eee0bd9a65b2a2887 *man/sm.monotonicity.Rd b23c695cd0faec1e12ccc2226354c8ca *man/sm.options.Rd 20b869b7b4759c869984359ef251175a *man/sm.pca.Rd @@ -230,8 +230,8 @@ bae9697d0ce71d91680037ec1d518a67 *man/sm.regression.autocor.Rd b60a62560fd2dddc4a842edee340e5b5 *man/sm.rm.Rd 0186bfb9ccbb37f33357f8087fd524d8 *man/sm.script.Rd -addd459f820b0b0ad4862836c26f1f9a *man/sm.sigma.Rd -b3dff66a7f4eed91ffc6d77318084934 *man/sm.sigma2.compare.Rd +5311e1d06f5ea5e1b50eb9837527de2a *man/sm.sigma.Rd +9ac1d4401071327d65f8e00a2e79df5f *man/sm.sigma2.compare.Rd 9a4f5e1491f10e5adf959ff7874d0481 *man/sm.sphere.Rd 11c40a2aab77a6e8fc84a2207609a691 *man/sm.surface3d.Rd d3fcd9d0f8c41c4df35ea713df6d1f7e *man/sm.survival.Rd diff -Nru sm-2.2-5.7/R/sm.r sm-2.2-5.7.1/R/sm.r --- sm-2.2-5.7/R/sm.r 2021-09-09 17:13:19.000000000 +0000 +++ sm-2.2-5.7.1/R/sm.r 2022-07-04 08:59:44.000000000 +0000 @@ -11,15 +11,15 @@ if (missing(model)) model <- "none" return(sm.regression(x, y, h, model = model, weights = weights, ...)) } - else if (class(x) != "formula") { + else if (!inherits(x, "formula")) { x.name <- deparse(substitute(x)) if (weights.missing) weights <- NA if (missing(model)) model <- "none" return(sm.density(x, h, model = model, weights = weights, xlab = x.name, ...)) } - + opt <- sm.options(list(...)) - + replace.na(opt, display, "lines") # replace.na(opt, reference, "none") opt$reference <- "none" @@ -43,7 +43,7 @@ formula.linear <- as.formula(formula.linear) names(vars.inf) <- rownames(involved) bricks.type <- sapply(vars.inf[-response.ind], mode) - ind <- (bricks.type == "numeric") & + ind <- (bricks.type == "numeric") & sapply(vars.inf[-response.ind], is.factor) bricks.type[ind] <- "factor" Xlinear <- vars.inf[-response.ind][bricks.type != "list"] @@ -65,12 +65,12 @@ fixed <- list() fac <- list() xmissing <- FALSE - + if (any(apply(involved, 2, sum) > 3)) stop("four-way interactions not yet implemented.") for (i in 1:length(terms.smooth)) { - + inv <- which(involved[ , terms.smooth[i]] == 1) ilinear <- which(bricks.type[names(inv)] == "numeric") ifactor <- which(bricks.type[names(inv)] == "factor") @@ -85,7 +85,7 @@ } else fac[[i]] <- NA - + nvars <- length(inv) X[[i]] <- matrix(nrow = length(y), ncol = 0) xlabels[[i]] <- vector("character") @@ -104,11 +104,11 @@ newvar <- eval.parent(parse(text = vars.inf[[j]]$variables)) if (is.matrix(newvar)) { nms <- colnames(newvar) - if (any(is.null(colnames(newvar)))) - nms <- paste(vars.inf[[j]]$variables, + if (any(is.null(colnames(newvar)))) + nms <- paste(vars.inf[[j]]$variables, "[", 1:ncol(newvar), "]", sep = "") } - else + else nms <- vars.inf[[j]]$variables xlab[[i]] <- c(xlab[[i]], nms) newvar <- as.matrix(newvar) @@ -119,7 +119,7 @@ if (length(prd) != ndims.new) stop("period does not match the columns of x.") period[[i]] <- c(period[[i]], prd) - if (any(!is.na(prd))) { + if (any(!is.na(prd))) { for (k in 1:ndims.new) if (!is.na(prd[k])) newvar[ , k] <- newvar[ , k] %% prd[k] } @@ -151,7 +151,7 @@ nseg[[i]] <- rep(switch(sum(ndims[[i]]), 100, 17, 7), sum(ndims[[i]])) if (any(is.na(X[[i]]))) xmissing <- TRUE } - + # Remove observations which have missing data. ind.missing <- lapply(X, function(x) apply(x, 1, function(z) any(is.na(z)))) ind.missing <- cbind(is.na(y), matrix(unlist(ind.missing), ncol = length(X))) @@ -161,15 +161,15 @@ for (i in 1:length(X)) X[[i]] <- as.matrix(X[[i]][!ind.missing, ]) cat("warning: missing data removed.\n") } - + if (opt$verbose > 1) tim <- proc.time() - + P <- list(length = length(terms.smooth)) B <- model.matrix(formula.linear, parent.frame()) m <- ncol(B) - + for (i in 1:length(terms.smooth)) { - mat <- ps.matrices(X[[i]], xrange[[i]], ndims = ndims[[i]], + mat <- ps.matrices(X[[i]], xrange[[i]], ndims = ndims[[i]], nseg = nseg[[i]], period = period[[i]]) if (all(is.na(fac[[i]]))) { B <- cbind(B, mat$B) @@ -193,7 +193,7 @@ } xrange[[i]] <- mat$xrange } - + if (opt$verbose > 1) { cat("Timings:\nconstructing matrices", (proc.time() - tim)[1], "seconds\n") tim <- proc.time() @@ -202,7 +202,7 @@ b.ind <- list(length = length(m)) for (i in 1:length(terms.smooth)) b.ind[[i]] <- (cumsum(m)[i] + 1):cumsum(m)[i + 1] - + if (weights.missing) { # btb <- t(B) %*% B btb <- crossprod(B) @@ -230,7 +230,7 @@ B1 <- solve(btb + lambda * P) sum(diag(btb %*% B1)) } - + for (i in 1:length(terms.smooth)) { if (any(is.na(df[[i]]))) df[[i]] <- switch(sum(ndims[[i]]), 6, 12, 18) # code doesn't currently handle more than one df for terms with more than one variable. @@ -240,7 +240,7 @@ if (any(is.na(lambda[[i]]))) lambda[[i]] <- lambda.select(btb[b.ind[[i]], b.ind[[i]]], bty[b.ind[[i]]], P[[i]], df[[i]]) } - + if (opt$verbose > 1) { cat("selecting smoothing parameters", (proc.time() - tim)[1], "seconds\n") tim <- proc.time() @@ -250,7 +250,7 @@ Pall <- matrix(0, nrow = ncol(B), ncol = ncol(B)) for (i in 1:length(terms.smooth)) Pall[b.ind[[i]], b.ind[[i]]] <- lambda[[i]] * P[[i]] - B1 <- solve(btb + Pall) + B1 <- solve(btb + Pall) # btb <- t(B[,-1]) %*% diag(rep(1, length(y))) %*% B[,-1] # return(list(btb = btb)) # B1 <- solve(btb[-1, -1] + Pall[-1, -1] - lambda[[i]] * mat$cmat) @@ -258,23 +258,23 @@ # alpha <- as.vector(B1 %*% bty[-1]) # mu <- c(B[, -1] %*% alpha) # return(list(alpha = alpha, B = B[ , -1], btb = btb[-1, -1])) - + if (opt$verbose > 1) { cat("fitting", (proc.time() - tim)[1], "seconds\n") tim <- proc.time() } - + # Force the estimate to pass through fixed points (1 covariate only) if (length(terms.smooth) == 1 & ndims[[1]] == 1 & all(!is.na(fixed[[1]]))) { fxd <- fixed[[1]] - if (any(fxd[,1] < xrange[[1]][1]) | + if (any(fxd[,1] < xrange[[1]][1]) | any(fxd[,1] > xrange[[1]][2])) stop("fixed points must be inside the range of the data.") A <- cbind(1, ps.matrices(as.matrix(fxd[ , 1]), xrange[[1]], ndims[[1]], nseg[[1]])$B) - alpha <- alpha + B1 %*% t(A) %*% solve(A %*% B1 %*% t(A)) %*% + alpha <- alpha + B1 %*% t(A) %*% solve(A %*% B1 %*% t(A)) %*% (fxd[ , 2] - A %*% alpha) - } + } mu <- c(B %*% alpha) df.model <- sum(btb * t(B1)) @@ -284,7 +284,7 @@ rss <- sum((y - mu)^2) tss <- sum((y - mean(y))^2) R.squared <- 100 * (tss - rss) / tss - + if (opt$verbose > 1) { cat("summaries", (proc.time() - tim)[1], "seconds\n") tim <- proc.time() @@ -293,22 +293,22 @@ # If there is only one term, include the mean # if (nterms == 1) b.ind[[1]] <- c(1, b.ind[[1]]) - result <- list(fitted = mu, alpha = alpha, m = m, B = B, + result <- list(fitted = mu, alpha = alpha, m = m, B = B, bty = bty, btb = btb, B1 = B1, Pall = Pall, xlabels = xlabels, linear.matrix = model.matrix(formula.linear, parent.frame()), terms.linear = terms.linear, terms.smooth = terms.smooth, xlab = xlab, ylab = ylab, term.labels = term.labels, - lambda = lambda, ndims = ndims, + lambda = lambda, ndims = ndims, y = y, X = X, fac = fac, Xlinear = Xlinear, bricks.type = bricks.type, sigma = sigma, cov.alpha = cov.alpha, b.ind = b.ind, - df = df, df.model = df.model, df.error = df.error, + df = df, df.model = df.model, df.error = df.error, rss = rss, R.squared = R.squared, xrange = xrange, nseg = nseg, bdeg = bdeg, period = period, pam.formula = pam.formula, involved = involved, nterms = nterms) if (!weights.missing) result$weights <- weights class(result) <- "pam" - + if (nterms == 1 & ndims[[1]] <= 2) { if (opt$panel) { replace.na(opt, df.max, switch(ndims[[1]], 20, 50, 100)) @@ -327,7 +327,7 @@ llambda.df <- c(llambda.df, lambda.df(exp(min(llambda)), btb, Pall)) } df.fun <- approxfun(llambda.df, llambda) - + sm.pam.draw <- function(pam.panel) { plot(pam.panel$model, options = pam.panel$opt) title(pam.panel$df) @@ -361,7 +361,7 @@ else if (opt$display != "none") plot(result, ...) } - + invisible(result) } @@ -421,25 +421,25 @@ inrange <- rep(TRUE, nnew) for (i in 1:model$nterms) for (j in 1:ncol(X[[i]])) - inrange <- inrange & X[[i]][ , j] >= model$xrange[[i]][j, 1] & - X[[i]][ , j] <= model$xrange[[i]][j, 2] + inrange <- inrange & X[[i]][ , j] >= model$xrange[[i]][j, 1] & + X[[i]][ , j] <= model$xrange[[i]][j, 2] outrange <- which(!inrange) if (length(outrange) > 0 & verbose > 0) warning("some evaluation points are out of range and have been omitted.") nnew <- length(which(inrange)) B <- rep(1, nnew) for (i in 1:model$nterms) { - mat <- ps.matrices(as.matrix(X[[i]][inrange, ]), xrange = model$xrange[[i]], + mat <- ps.matrices(as.matrix(X[[i]][inrange, ]), xrange = model$xrange[[i]], ndims = model$ndims[[i]], nseg = model$nseg[[i]], period = model$period[[i]]) B <- cbind(B, mat$B) } - + fv <- rep(NA, nnew) fv[inrange] <- c(B %*% model$alpha) if (model$nterms == 1 & model$ndims[[1]] == 1 & deriv > 0) { - mat <- ps.matrices(as.matrix(X[[1]][inrange, ]), xrange = model$xrange[[1]], - ndims = model$ndims[[1]], nseg = model$nseg[[1]], bdeg = model$bdeg - deriv, + mat <- ps.matrices(as.matrix(X[[1]][inrange, ]), xrange = model$xrange[[1]], + ndims = model$ndims[[1]], nseg = model$nseg[[1]], bdeg = model$bdeg - deriv, period = model$period[[1]]) alpha1 <- diff(model$alpha[-1], differences = deriv) h <- model$xrange[[1]][,2] - model$xrange[[1]][,1] @@ -447,10 +447,10 @@ } results <- list(fit = fv, inrange = inrange, B = B) - + if (se.fit) results$se.fit = sqrt(diag(B %*% model$cov.alpha %*% t(B))) - + return(invisible(results)) } @@ -477,11 +477,11 @@ #---------------------------------------------------------------------------- sm.mask <- function(x, eval.points, mask.method = "hull") { - + ngrid <- nrow(eval.points) - grid.points <- cbind(rep(eval.points[, 1], ngrid), + grid.points <- cbind(rep(eval.points[, 1], ngrid), rep(eval.points[, 2], rep(ngrid, ngrid))) - + if (mask.method == "hull") { hull.points <- as.matrix(x[order(x[, 1], x[, 2]), ]) dh <- diff(hull.points) @@ -516,7 +516,7 @@ } else mask <- matrix(1, ncol = ngrid, nrow = ngrid) - + return(invisible(mask)) } @@ -540,19 +540,19 @@ # Compute a set of basis functions and a penalty matrix associated with x. # An intercept term and the main effect of any interaction terms are removed. - + ndimx <- ncol(x) if (ndimx > 3) stop("terms with more than three dimensions cannot be used.") n <- nrow(x) - + if (missing(nseg)) nseg <- rep(switch(ndimx, 100, 17, 7), ndimx) - + # Compute B-spline basis - + b <- list(length = ndimx) m <- vector(length = ndimx) for (i in 1:ndimx) { - b[[i]] <- bbase(x[,i], xl = xrange[i , 1], xr = xrange[i, 2], nseg = nseg[i], + b[[i]] <- bbase(x[,i], xl = xrange[i , 1], xr = xrange[i, 2], nseg = nseg[i], deg = bdeg) m[i] <- ncol(b[[i]]) } @@ -562,9 +562,9 @@ B <- t(apply(cbind(b[[1]], b[[2]]), 1, function(x) c(x[1:m[1]] %x% x[-(1:m[1])]))) if (ndimx == 3) - B <- t(apply(cbind(B, b[[3]]), 1, + B <- t(apply(cbind(B, b[[3]]), 1, function(x) c(x[1:(m[1]*m[2])] %x% x[-(1:(m[1]*m[2]))]))) - + # Construct smoothness penalty matrices P <- list() for (i in 1:ndimx) { @@ -613,7 +613,7 @@ pmat <- pmat + diag(m[1]) %x% matrix(1, nrow = m[2], ncol = m[2]) %x% diag(m[3]) pmat <- pmat + diag(m[1]) %x% diag(m[2]) %x% matrix(1, nrow = m[3], ncol = m[3]) } - + result <- list(B = B, P = pmat, xrange = xrange, nseg = nseg, bdeg = bdeg, pord = pord) if (length(ndims) == 1) result$cmat <- cmat invisible(result) diff -Nru sm-2.2-5.7/R/variogram.r sm-2.2-5.7.1/R/variogram.r --- sm-2.2-5.7/R/variogram.r 2021-09-13 09:58:24.000000000 +0000 +++ sm-2.2-5.7.1/R/variogram.r 2022-07-01 07:58:02.000000000 +0000 @@ -1,16 +1,16 @@ "sm.variogram" <- function(x, y, h, - df.se = "automatic", max.dist = NA, n.zero.dist = 1, original.scale = TRUE, + df.se = "automatic", max.dist = NA, n.zero.dist = 1, original.scale = TRUE, varmat = FALSE, ...) { type <- "binned" bin.type <- "log" type.se <- "smooth-monotonic-original" - + if ("geodata" %in% class(x)) { y <- x$data x <- x$coords } - + opt <- sm.options(list(...)) data <- sm.check.data(x = x, y = y, ...) x <- data$x @@ -19,7 +19,7 @@ ndim <- data$ndim opt <- data$options # rawdata <- list(x = x, y = y, nbins = opt$nbins, nobs = n, ndim = ndim) - + model <- opt$model if (!(model %in% c("none", "independent", "isotropic", "stationary"))) { cat("The 'model' argument is not recognised - reverting to 'none'.\n") @@ -80,23 +80,23 @@ i1 <- (as.vector(imat))[ind > 0] i2 <- (as.vector(t(imat)))[ind > 0] ipair <- cbind(i1, i2) - + # if (!is.na(max.dist)) { # ind <- (hall <= max.dist) # hall <- hall[ind] # dall <- dall[ind] # ipair <- ipair[ind, ] # } - + results <- list(distance = hall, sqrtdiff = dall, ipair = ipair) - - + + #------------------------------------------------------------ # Bin the differences #------------------------------------------------------------ if (!(model %in% c("isotropic", "stationary"))) { - + if (bin.type == "regular") { bins <- binning(hall, dall, nbins = opt$nbins) hh <- bins$x @@ -125,7 +125,7 @@ dd <- tapply(dall, ibin, mean) # hh <- breaks[-1] - diff(breaks) / 2 hh <- tapply(results$distance, ibin, mean) - wts <- table(ibin) + wts <- table(ibin) } else if (bin.type == "unique") { hh <- sort(unique(signif(hall, 6))) @@ -138,7 +138,7 @@ ind <- (hall < 2 * .Machine$double.eps) nzero <- length(which(ind)) if (nzero >= max(n.zero.dist, 1)) breaks <- c(breaks, 2 * .Machine$double.eps) - breaks <- c(breaks, exp(min(log(hall[!ind])) + + breaks <- c(breaks, exp(min(log(hall[!ind])) + (1:opt$nbins) * diff(range(log(hall[!ind]))) / opt$nbins)) nbrks <- length(breaks) breaks[nbrks] <- breaks[nbrks] + 1 @@ -159,7 +159,7 @@ results$ibin <- ibin results$breaks <- breaks results$nbins <- length(hh) - + nbins <- length(results$sqrtdiff.mean) if (!is.numeric(df.se)) df.se <- round(0.8 * nbins) @@ -173,7 +173,7 @@ gamma.hat <- dd } else { - if (missing(h)) h <- h.select(hh, dd, weights = wts, + if (missing(h)) h <- h.select(hh, dd, weights = wts, nbins = 0, df = df.se, method = opt$method) replace.na(opt, eval.points, seq(min(hall), max(hall), length = opt$ngrid)) ev <- opt$eval.points @@ -196,14 +196,14 @@ if (type.se == "cressie") gamma.hat.V <- 0.5 * dd^4 / (0.457 + 0.494 / wts) if (type.se == "smooth-monotonic-original") { - sm.model <- ps.normal(hh, 0.5 * dd0, df = df.se, weights = wts, + sm.model <- ps.normal(hh, 0.5 * dd0, df = df.se, weights = wts, eval.points = hh, increasing = TRUE, display = "none") gamma.hat.V <- sm.model$estimate # results$gamma.hat.V <- gamma.hat.V } if (type.se == "smooth-original") { - sm.model <- ps.normal(hh, 0.5 * dd0, df = df.se, weights = wts, + sm.model <- ps.normal(hh, 0.5 * dd0, df = df.se, weights = wts, eval.points = hh, increasing = FALSE, display = "none") gamma.hat.V <- sm.model$estimate @@ -212,12 +212,12 @@ if (type.se %in% c("smooth", "smooth-monotonic", "smooth-w", "smooth-monotonic-w")) { # ind <- (hh <= max.dist) # gamma.hat.V <- rep(0, length(hh)) - # gamma.hat.V[ind] <- sm.regression(hh, dd, weights = wts, eval.points = hh[ind], + # gamma.hat.V[ind] <- sm.regression(hh, dd, weights = wts, eval.points = hh[ind], # display = "none")$estimate # gamma.hat.V[!ind] <- gamma.hat.V[hh == max(hh[ind])] inc <- (type.se %in% c("smooth-monotonic", "smooth-monotonic-w")) wp <- (type.se %in% c("smooth-w", "smooth-monotonic-w")) - sm.model <- ps.normal(hh, dd, df = df.se, weights = wts, + sm.model <- ps.normal(hh, dd, df = df.se, weights = wts, eval.points = hh, increasing = inc, kappa = 1e8, weights.penalty = wp, display = "none") gamma.hat.V <- sm.model$estimate @@ -246,12 +246,12 @@ # messages = FALSE) # gamma.hat.V <- variofit(vgm.emp, ini = c(var(y), 0.2), fix.nugget = TRUE, fix.kappa = TRUE, # max.dist = max.dist, messages = FALSE) - # results$cov.pars <- gamma.hat.V$cov.pars + # results$cov.pars <- gamma.hat.V$cov.pars # results$kappa <- gamma.hat.V$kappa # # plot(vgm.emp) # # lines(gamma.hat) - # gamma.hat.V <- gamma.hat.V$cov.pars[1] - - # cov.spatial(results$distance.mean, cov.pars = gamma.hat.V$cov.pars, + # gamma.hat.V <- gamma.hat.V$cov.pars[1] - + # cov.spatial(results$distance.mean, cov.pars = gamma.hat.V$cov.pars, # kappa = gamma.hat.V$kappa) # } # Smoothing - small df and tilted to small distances @@ -260,9 +260,9 @@ # gamma.hat <- sm.regression(hh, dd, df = 4, weights = wts, h.weights = hw, # eval.points = hh, display = "none")$estimate # gamma.hat <- (gamma.hat / 0.977741)^4 - + gamma.hat.V <- pmax(gamma.hat.V, 0) - + if (type == "binned") { se <- rep(0, nbins) # LAST MODIFICATION @@ -292,8 +292,8 @@ # for (i in 1:nbins) { # # if (opt$verbose > 0) cat(i, "") # for (j in i:nbins){ - # # ib <- sort(unique(results$ibin))[i] - # # jb <- sort(unique(results$ibin))[j] + # # ib <- sort(unique(results$ibin))[i] + # # jb <- sort(unique(results$ibin))[j] # V[i, j] <- cov.bin.fun(i, j, results, gamma.hat.V) # if (j > i) V[j, i] <- V[i, j] # } @@ -324,7 +324,7 @@ #------------------------------------------------------------ if (model == "independent") { - + vv <- 0.1724 cv <- 0.03144 Sigma <- table(c(igp, igp), c(i1, i2)) @@ -345,12 +345,12 @@ A <- t(diag(nb) - A) %*% diag(wts) %*% (diag(nb) - A) A <- A - (1 + tobs) * t(diag(nb) - W) %*% diag(wts) %*% (diag(nb) - W) pval <- p.quad.moment(A, Sigma, 0, 0) - if (opt$verbose > 0) + if (opt$verbose > 0) cat("Test of spatial independence: p = ",round(pval, 3), "\n") results$h <- h results$p <- pval } - + } #------------------------------------------------------------ @@ -358,9 +358,9 @@ #------------------------------------------------------------ if ((opt$display != "none") & !(model %in% c("isotropic", "stationary"))) { - + fn <- if (original.scale) function(x) (x / 0.977741)^4 else I - + if (opt$display %in% c("bins", "means")) { xx <- hh yy <- fn(dd) @@ -369,9 +369,9 @@ xx <- hall yy <- if (original.scale) dall^4 else dall } - + if (!opt$add) { - + replace.na(opt, xlab, "Distance") replace.na(opt, ylab, if (original.scale) " Squared difference" else "Square-root abs. difference") replace.na(opt, xlim, range(xx)) @@ -418,7 +418,7 @@ # lines(ev, gamma.hat - 2 * se, lty = 2, col = opt$col) # } # } - + } #------------------------------------------------------------ @@ -432,7 +432,7 @@ amat <- atan2(t(x2mat) - x2mat, t(x1mat) - x1mat) amat[amat < 0] <- amat[amat < 0] + pi angles <- (as.vector(amat))[ind > 0] - + centres1 <- seq(0, pi, length = opt$nbins + 1) centres1 <- centres1[-1] - diff(centres1) / 2 ind <- (hall < 2 * .Machine$double.eps) @@ -450,7 +450,7 @@ if (length(ind.pt) > 1) ind.pt <- ind.pt[gpt[, 2] == min(gpt[, 2])] ind.pt } - ibin <- apply(cbind(angles, hall), 1, identify.grid, centres) + ibin <- apply(cbind(angles, hall), 1, identify.grid, centres) ibin <- match(ibin, unique(ibin)) dd <- as.vector(tapply(dall, ibin, mean)) dd0 <- as.vector(tapply(dall0, ibin, mean)) @@ -461,10 +461,10 @@ # sm.regression(cbind(hh, ang), dd, df = 12, weights = wts, nbins = 0, display = "rgl", # col.points = "red", alpha = 0.7, size = 2, period = c(NA, pi)) - + h <- h.select(cbind(hh, ang), dd, weights = wts, df = opt$df, nbins = 0, period = c(NA, pi)) if (!is.numeric(df.se)) df.se <- round(0.8 * opt$nbins) - + results$distance.mean <- hh results$sqrtdiff.mean <- dd results$angles <- angles @@ -472,36 +472,36 @@ results$weights <- wts results$ibin <- ibin results$h <- h - + if (is.na(max.dist)) max.dist <- max(hh) + 1 ind <- (hh <= max.dist) gamma.hat.V <- rep(0, length(hh)) - + if (type.se == "binned") gg <- dd[ind] else if (type.se == "smooth") - gg <- sm.regression(hh, dd, weights = wts, eval.points = hh[ind], + gg <- sm.regression(hh, dd, weights = wts, eval.points = hh[ind], display = "none", nbins = 0)$estimate else if (type.se == "smooth-monotonic") gg <- ps.normal(hh, dd, df = df.se, weights = wts, eval.points = hh[ind], increasing = TRUE, weights.penalty = FALSE, display = "none")$estimate else if (type.se == "smooth-monotonic-original") - gg <- ps.normal(hh, 0.5 * dd0, df = df.se, weights = wts, + gg <- ps.normal(hh, 0.5 * dd0, df = df.se, weights = wts, negative = FALSE, increasing = TRUE, eval.points = hh[ind], display = "none")$estimate - + gamma.hat.V[ind] <- pmax(gg, 0) gamma.hat.V[!ind] <- gamma.hat.V[hh == max(hh[ind])] # results$gamma.hat.V <- gamma.hat.V if (type.se != "smooth-monotonic-original") gamma.hat.V <- (gamma.hat.V / 0.977741)^4 - + # plot(hh, dd) # plot(hh, 0.5 * dd0) # ps.normal(hh, 0.5 * dd0, df = df.se, weights = wts, negative = FALSE, increasing = TRUE) # print(cbind(hh, gamma.hat.V)) # stop() - + V <- matrix(0, nrow = nbins, ncol = nbins) @@ -512,8 +512,8 @@ # for (i in 1:nbins) { # if (opt$verbose > 0) cat(i, "") # for (j in i:nbins){ -# # ib <- sort(unique(results$ibin))[i] -# # jb <- sort(unique(results$ibin))[j] +# # ib <- sort(unique(results$ibin))[i] +# # jb <- sort(unique(results$ibin))[j] # V[i, j] <- cov.bin.fun(i, j, results, gamma.hat.V) # if (j > i) V[j, i] <- V[i, j] # } @@ -532,11 +532,11 @@ V <- matrix(data=output$res, nrow=length(gamma.hat.V), ncol=length(gamma.hat.V), byrow= TRUE) # if (opt$verbose > 0) cat("\n") - + opt$period <- c(NA, pi) - + if (opt$test | (opt$display != "none")) { - mdl1 <- sm(dd ~ s(cbind(hh, ang), df = opt$df, period = opt$period), + mdl1 <- sm(dd ~ s(cbind(hh, ang), df = opt$df, period = opt$period), weights = wts, display = "none") df0 <- ceiling(sqrt(opt$df)) mdl0 <- sm(dd ~ s(hh, df = df0), weights = wts, display = "none") @@ -544,16 +544,16 @@ # mdl0 <- sm(dd ~ s(hh, df = opt$df / 2), weights = wts, display = "none") # mdl0 <- sm(dd ~ s(hh, df = opt$df / 3), weights = wts, display = "none") # mdl0 <- sm(dd ~ s(hh, df = 4), weights = wts, display = "none") - # mdl0 <- sm(dd ~ s(hh, lambda = mdl$lambda[[1]][1] * (mdl$nseg[[1]][1] + 3)), + # mdl0 <- sm(dd ~ s(hh, lambda = mdl$lambda[[1]][1] * (mdl$nseg[[1]][1] + 3)), # weights = wts, display = "none") } - + if (opt$test) { - + # save(model, V, wts, dd, hh, ang, h, opt, file = "temp.dmp") - + S1 <- mdl1$B %*% mdl1$B1 %*% t(mdl1$B * wts) - S0 <- mdl0$B %*% mdl0$B1 %*% t(mdl0$B * wts) + S0 <- mdl0$B %*% mdl0$B1 %*% t(mdl0$B * wts) # est1 <- S0 %*% dd # S1 <- S0 - S1 @@ -562,8 +562,8 @@ # ds <- S1 %*% dd # tobs <- c(t(ds) %*% V1 %*% ds) # pval1 <- p.quad.moment.adjusted(t(S1) %*% V1 %*% S1, V, tobs) -# - # S0 <- sm.weight(hh, hh, h[1], weights = wts) +# + # S0 <- sm.weight(hh, hh, h[1], weights = wts) # S1 <- sm.weight2(cbind(hh, ang), cbind(hh, ang), h, weights = wts, options = opt) # cat("traces:", round(sum(diag(S0))), round(sum(diag(S1))), "\n") # est0 <- S0 %*% dd @@ -587,10 +587,10 @@ ds <- S1 %*% dd tobs <- c(t(ds) %*% V1 %*% ds) pval <- p.quad.moment.adjusted(t(S1) %*% V1 %*% S1, V, tobs) - + # cat(round(pval, 3), round(pval1, 3), "\n") - -# mdl <- sm(dd ~ s(hh, df = 4) * s(ang, df = 4, period = pi), + +# mdl <- sm(dd ~ s(hh, df = 4) * s(ang, df = 4, period = pi), # weights = wts, display = "none") # # save(model, V, wts, file = "temp.dmp") # ind <- c(mdl$b.ind[[2]], mdl$b.ind[[3]]) @@ -599,7 +599,7 @@ # I.i[ind] <- 1 # A <- mdl$B1 %*% t(mdl$B * wts) # pval <- p.quad.moment.adjusted(t(A) %*% diag(I.i) %*% A, V, tobs) - + if (opt$verbose > 0) cat("Test of isotropy: p = ", round(pval, 3), "\n") results$h <- h results$p <- pval @@ -615,33 +615,33 @@ replace.na(op, zlab, "Square-root difference") replace.na(op, ylab, "Angle") op$display <- "none" - + u <- list(length = 2) for (j in 1:2) u[[j]] <- seq(mdl1$xrange[[1]][j, 1], mdl1$xrange[[1]][j, 2], length = opt$ngrid) U <- as.matrix(expand.grid(u)) mask <- sm.mask(cbind(hh, ang), cbind(u[[1]], u[[2]]), mask.method = opt$mask.method) - + B1 <- ps.matrices(U, mdl1$xrange[[1]], 2, nseg = mdl1$nseg[[1]], period = mdl1$period[[1]])$B B1 <- cbind(1, B1) S1 <- B1 %*% mdl1$B1 %*% t(mdl1$B * wts) est1 <- matrix(S1 %*% dd, nrow = opt$ngrid) * mask - + B0 <- ps.matrices(as.matrix(U[ , 1]), mdl0$xrange[[1]], 1, nseg = mdl0$nseg[[1]])$B B0 <- cbind(1, B0) S0 <- B0 %*% mdl0$B1 %*% t(mdl0$B * wts) est0 <- matrix(S0 %*% dd, nrow = opt$ngrid) * mask - + stde <- sqrt(diag((S1 - S0) %*% V %*% t(S1 - S0))) sdiff <- (est1 - est0) / matrix(stde, ncol = opt$ngrid) ev <- u gamma.hat <- est1 - + results$eval.points <- u results$estimate <- est1 results$sdiff <- sdiff - # surf <- sm.regression.2d(cbind(hh, ang), dd, h, model = "isotropic", weights = wts, + # surf <- sm.regression.2d(cbind(hh, ang), dd, h, model = "isotropic", weights = wts, # rawdata = list(x = cbind(hh, ang), y = dd), options = op) # ngrid <- nrow(surf$eval.points) # ev <- rep(surf$eval.points[, 1], ngrid) @@ -655,7 +655,7 @@ # alpha = 0.5, alpha.mesh = 0) } else if (opt$display == "image") { - if (!requireNamespace("akima", quietly = TRUE)) stop("this option requires the akima package.") + if (!requireNamespace("interp", quietly = TRUE)) stop("this option requires the interp package.") a <- U a <- rbind(a, cbind(a[ , 1], a[ , 2] + pi)) a1 <- a[ , 1] * cos(a[ , 2]) @@ -663,8 +663,8 @@ b <- rep(c(est1), 2) sdiff <- rep(c(sdiff), 2) ind <- !is.na(b) & !duplicated(cbind(a1, a2)) - inte <- akima::interp(a1[ind], a2[ind], b[ind]) - ints <- akima::interp(a1[ind], a2[ind], sdiff[ind]) + inte <- interp::interp(a1[ind], a2[ind], b[ind]) + ints <- interp::interp(a1[ind], a2[ind], sdiff[ind]) cts <- contourLines(ints) lvls <- rep(0, length(cts)) for (i in 1:length(cts)) lvls[i] <- cts[[i]]$level @@ -676,7 +676,7 @@ axis(2) # op <- opt # op$display <- "none" - # surf <- sm.regression.2d(cbind(hh, ang), dd, h, model = "isotropic", weights = wts, + # surf <- sm.regression.2d(cbind(hh, ang), dd, h, model = "isotropic", weights = wts, # rawdata = list(x = cbind(hh, ang), y = dd), options = op) if (length(lvls) > 0) contour(ints, levels = lvls, add = TRUE, col = "red") # a <- as.matrix(expand.grid(surf$eval.points[ , 1], surf$eval.points[ , 2])) @@ -693,7 +693,7 @@ # axis(2) # op <- opt # op$display <- "none" -# surf <- sm.regression.2d(cbind(hh, ang), dd, h, model = "isotropic", weights = wts, +# surf <- sm.regression.2d(cbind(hh, ang), dd, h, model = "isotropic", weights = wts, # rawdata = list(x = cbind(hh, ang), y = dd), options = op) # sdiff <- rep(c(surf$sdiff), 2) # cts <- contourLines(interp(a1[ind], a2[ind], sdiff[ind])) @@ -704,7 +704,7 @@ } ) } - + } #------------------------------------------------------------ @@ -739,7 +739,7 @@ if (length(ind.pt) > 1) ind.pt <- ind.pt[gpt[, 3] == min(gpt[, 3])] ind.pt } - ibin <- apply(cbind(av1, av2, hall), 1, identify.grid, centres) + ibin <- apply(cbind(av1, av2, hall), 1, identify.grid, centres) ibin <- match(ibin, unique(ibin)) dd <- as.vector(tapply(dall, ibin, mean)) dd0 <- as.vector(tapply(dall0, ibin, mean)) @@ -748,7 +748,7 @@ a2 <- as.vector(tapply(av2, ibin, mean)) wts <- as.vector(table(ibin)) nbins <- length(unique(ibin)) - + # h <- h.select(cbind(hh, ang), dd, weights = wts, df = opt$df, nbins = 0, period = c(NA, pi)) if (!is.numeric(df.se)) df.se <- round(0.8 * opt$nbins) @@ -760,32 +760,32 @@ # results$a2 <- a2 results$weights <- wts results$ibin <- ibin - + xx <- cbind(x = a1, y = a2, Distance = hh) - + if (is.na(max.dist)) max.dist <- max(hh) + 1 ind <- (hh <= max.dist) gamma.hat.V <- rep(0, length(hh)) - + if (type.se == "binned") gg <- dd[ind] else if (type.se == "smooth") - gg <- sm.regression(hh, dd, weights = wts, eval.points = hh[ind], + gg <- sm.regression(hh, dd, weights = wts, eval.points = hh[ind], display = "none", nbins = 0)$estimate else if (type.se == "smooth-monotonic") gg <- ps.normal(hh, dd, df = df.se, weights = wts, eval.points = hh, increasing = TRUE, weights.penalty = FALSE, display = "none")$estimate else if (type.se == "smooth-monotonic-original") - gg <- ps.normal(hh, 0.5 * dd0, df = df.se, weights = wts, + gg <- ps.normal(hh, 0.5 * dd0, df = df.se, weights = wts, negative = FALSE, increasing = TRUE, eval.points = hh[ind], display = "none")$estimate - + gamma.hat.V[ind] <- pmax(gg, 0) gamma.hat.V[!ind] <- gamma.hat.V[hh == max(hh[ind])] # results$gamma.hat.V <- gamma.hat.V if (type.se != "smooth-monotonic-original") gamma.hat.V <- (gamma.hat.V / 0.977741)^4 - + V <- matrix(0, nrow = nbins, ncol = nbins) # if (opt$verbose > 0) cat(nbins, ": ") @@ -817,13 +817,13 @@ # if (opt$verbose > 0) cat("\n") model1 <- sm(dd ~ s(xx, df = opt$df), weights = wts, display = "none") - + if (opt$test) { df0 <- ceiling(opt$df^(1/3)) # df0 <- ceiling(opt$df / 3) # df0 <- ceiling(opt$df / 2) model0 <- sm(dd ~ s(hh, df = df0), weights = wts, display = "none") - # mdl0 <- sm(dd ~ s(hh, lambda = mdl$lambda[[1]][1] * (mdl$nseg[[1]][1] + 3)^2), + # mdl0 <- sm(dd ~ s(hh, lambda = mdl$lambda[[1]][1] * (mdl$nseg[[1]][1] + 3)^2), # weights = wts, display = "none") # model0 <- sm(dd ~ s(hh, df = 3), weights = wts, display = "none") S0 <- model0$B %*% model0$B1 %*% t(model0$B * wts) @@ -832,7 +832,7 @@ S0 <- t(S0 - S1) %*% V1 %*% (S0 - S1) tobs <- c(dd %*% S0 %*% dd) pval <- p.quad.moment.adjusted(S0, V, tobs) - if (opt$verbose > 0) + if (opt$verbose > 0) cat("Test of stationarity: p = ", round(pval, 3), "\n") results$p <- pval # results$df0 <- df0 @@ -841,18 +841,18 @@ u <- list(length = 3) for (j in 1:3) u[[j]] <- seq(model1$xrange[[1]][j, 1], model1$xrange[[1]][j, 2], length = opt$ngrid) - U <- as.matrix(expand.grid(u)) + U <- as.matrix(expand.grid(u)) B1 <- ps.matrices(U, model1$xrange[[1]], 3, nseg = model1$nseg[[1]], period = model1$period[[1]])$B B1 <- cbind(1, B1) S1 <- B1 %*% model1$B1 %*% t(model1$B * wts) est1 <- S1 %*% dd - + B0 <- ps.matrices(as.matrix(U[ , 3]), model0$xrange[[1]], 1, nseg = model0$nseg[[1]])$B B0 <- cbind(1, B0) S0 <- B0 %*% model0$B1 %*% t(model0$B * wts) est0 <- S0 %*% dd - + stde <- diag((S1 - S0) %*% V %*% t(S1 - S0)) stde[stde < 0] <- NA stde <- sqrt(stde) @@ -865,14 +865,14 @@ # Add V into the estimated surface. # Return the estimate from the fitted surface. # plot(model1) - + # replace.na(opt, xlab, "Distance") # replace.na(opt, zlab, "Square-root difference") # replace.na(opt, ylab, "Angle") gamma.hat <- model1$fitted ev <- xx # results$model <- model1 - + names(u) <- c("x1", "x2", "distance") results$eval.points <- u results$estimate <- array(est1, dim = rep(opt$ngrid, 3)) @@ -889,12 +889,12 @@ # results$estimate <- gamma.hat # results$gamma.hat.V <- gamma.hat.V results$df <- opt$df - + if (opt$se & model == "none") results$se <- se if ((model == "independent") & opt$band) results$se.band <- se.band if (varmat | model %in% c("isotropic", "stationary")) results$V <- V invisible(results) - + }