File: ucv.R

package info (click to toggle)
r-cran-mass 7.3-51.1-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster
  • size: 2,148 kB
  • sloc: ansic: 664; makefile: 2
file content (129 lines) | stat: -rw-r--r-- 3,897 bytes parent folder | download | duplicates (7)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
# file MASS/R/ucv.R
# copyright (C) 1994-9 W. N. Venables and B. D. Ripley
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 or 3 of the License
#  (at your option).
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  http://www.r-project.org/Licenses/
#

width.SJ <- function(x, nb=1000, lower=0.1*hmax, upper=hmax,
		     method = c("ste", "dpi"))
{
    fSD <- function(h, x, alph2, c1, n, d)
        (c1/SDh(x, alph2 * h^(5/7), n, d))^(1/5) - h
    SDh <- function(x, h, n, d)
        .C(VR_phi4_bin,
           as.integer(n),
           as.integer(length(x)),
           as.double(d),
           x,
           as.double(h),
           u = double(1))$u
    TDh <- function(x, h, n, d)
        .C(VR_phi6_bin,
           as.integer(n),
           as.integer(length(x)),
           as.double(d),
           x,
           as.double(h),
           u = double(1))$u

    method <- match.arg(method)
    n <- length(x)
    if(!n) stop("'x' has length zero")
    storage.mode(x) <- "double"
    Z <- .C(VR_den_bin,
            as.integer(n),
            as.integer(nb),
            d = double(1),
            x,
            cnt = integer(nb)
            )
    d <- Z$d; cnt <- as.integer(Z$cnt)
    hmax <- 1.144 * sqrt(var(x)) * n^(-1/5)
    scale <- min(sqrt(var(x)), IQR(x)/1.349)
    a <- 1.24 * scale * n^(-1/7)
    b <- 1.23 * scale * n^(-1/9)
    c1 <- 1/(2*sqrt(pi)*n)
    TD  <- -TDh(cnt, b, n, d)
    alph2 <- 1.357*(SDh(cnt, a, n, d)/TD)^(1/7)
    if(method == "dpi")
        res <- (c1/SDh(cnt,(2.394/(n * TD))^(1/7) , n, d))^(1/5)
    else {
        if (fSD(lower, cnt, alph2, c1, n, d) *
            fSD(upper, cnt, alph2, c1, n, d) > 0)
            stop("no solution in the specified range of bandwidths")
        res <- uniroot(fSD, c(lower, upper), tol=0.1*lower,
                       x=cnt, alph2=alph2, c1=c1, n=n, d=d)$root
    }
    4 * res
}


ucv <- function(x, nb=1000, lower=0.1*hmax, upper=hmax)
{
    fucv <- function(h, x, n, d)
        .C(VR_ucv_bin,
           as.integer(n),
           as.integer(length(x)),
           as.double(d),
           x,
           as.double(h),
           u = double(1))$u

    n <- length(x)
    if(!n) stop("'x' has length zero")
    hmax <- 1.144 * sqrt(var(x)) * n^(-1/5) * 4
    storage.mode(x) <- "double"
    Z <- .C(VR_den_bin,
            as.integer(n),
            as.integer(nb),
            d = double(1),
            x,
            cnt = integer(nb))
    d <- Z$d; cnt <- as.integer(Z$cnt)
    h <- optimize(fucv, c(lower, upper), tol=0.1*lower,
                  x=cnt, n=n, d=d)$minimum
    if(h < 1.1*lower | h > upper-0.1*lower)
        warning("minimum occurred at one end of the range")
    h
}

bcv <- function(x, nb=1000, lower=0.1*hmax, upper=hmax)
{
    fbcv <- function(h, x, n, d)
        .C(VR_bcv_bin,
           as.integer(n),
           as.integer(length(x)),
           as.double(d),
           x,
           as.double(h),
           u = double(1))$u

    n <- length(x)
    if(!n) stop("'x' has length zero")
    hmax <- 1.144 * sqrt(var(x)) * n^(-1/5) * 4
    storage.mode(x) <- "double"
    Z <- .C(VR_den_bin,
            as.integer(n),
            as.integer(nb),
            d = double(1),
            x,
            cnt = integer(nb)
            )
    d <- Z$d; cnt <- as.integer(Z$cnt)
    h<- optimize(fbcv, c(lower, upper), tol=0.1*lower,
                 x=cnt, n=n, d=d)$minimum
    if(h < 1.1*lower | h > upper-0.1*lower)
        warning("minimum occurred at one end of the range")
    h
}