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.