1  | 
Z | 
|
mu | 
|
Sigma | 
|
theta | 
|
d | 
|
r | 
|
tau | 
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)
  }
 | 
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.