SErob:

Usage Arguments Examples

Usage

1
SErob(Z, mu, Sigma, theta, d, r, tau)

Arguments

Z
mu
Sigma
theta
d
r
tau

Examples

 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
##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (Z, mu, Sigma, theta, d, r, tau) 
{
    n = nrow(Z)
    p = ncol(Z)
    ps = p * (p + 1)/2
    q = 6
    Dup = Dp(p)
    DupPlus = solve(t(Dup) %*% Dup) %*% t(Dup)
    InvSigma = solve(Sigma)
    sigma = vech(Sigma)
    W = 0.5 * t(Dup) %*% (InvSigma %x% InvSigma) %*% Dup
    Zr = matrix(NA, n, p)
    A = matrix(0, p + q, p + q)
    B = matrix(0, p + q, p + q)
    sumg = rep(0, p + q)
    for (i in 1:n) {
        zi = Z[i, ]
        zi0 = zi - mu
        di = d[i]
        if (di <= r) {
            u1i = 1
            u2i = 1/tau
            du1i = 0
            du2i = 0
        }
        else {
            u1i = r/di
            u2i = u1i^2/tau
            du1i = -r/di^2
            du2i = -2 * r^2/tau/di^3
        }
        Zr[i, ] = sqrt(u2i) * t(zi0)
        g1i = u1i * zi0
        vTi = vech(zi0 %*% t(zi0))
        g2i = u2i * vTi - sigma
        gi = rbind(g1i, g2i)
        sumg = gi + sumg
        B = B + gi %*% t(gi)
        ddmu = -1/di * t(zi0) %*% InvSigma
        ddsigma = -t(vTi) %*% W/di
        dg1imu = -u1i * diag(p) + du1i * zi0 %*% ddmu
        dg1isigma = du1i * zi0 %*% ddsigma
        dg2imu = -u2i * DupPlus %*% (zi0 %x% diag(p) + diag(p) %x% 
            zi0) + du2i * vTi %*% ddmu
        dg2isigma = du2i * vTi %*% ddsigma - diag(q)
        dgi = rbind(cbind(dg1imu, dg1isigma), cbind(dg2imu, dg2isigma))
        A = A + dgi
    }
    A = -1 * A/n
    B = B/n
    invA = solve(A)
    OmegaSW = invA %*% B %*% t(invA)
    OmegaSW = OmegaSW[(p + 1):(p + q), (p + 1):(p + q)]
    SEsw = getSE(theta, OmegaSW, n)
    SEinf = SEML(Zr, theta)$inf
    Results = list(inf = SEinf, sand = SEsw, Zr = Zr)
    return(Results)
  }

musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.