# R/x_probmv.r In probhat: Multivariate Generalized Kernel Smoothing and Related Statistical Methods

```#probhat: Multivariate Generalized Kernel Smoothing and Related Statistical Methods
#Copyright (C), Abby Spurdle, 2019 to 2021

#This program is distributed without any warranty.

#This program is free software.
#You can modify it and/or redistribute it, under the terms of:
#The GNU General Public License, version 2, or (at your option) any later version.

#Also, this license should be available at:

.binary.seq = function (m, n)
{	rep.vector = round (2 ^ (0:(m - 1) ) )
verts = matrix (0, nrow=n, ncol=m)
for (j in 1:m)
verts [,j] = rep (c (0, 1), times=rep.vector [j], each=rep.vector [1 + m - j])
verts
}

.binary.sign = function (binary)
{   sums = apply (binary, 1, sum)
sign = rep (1L, nrow (binary) )
if (ncol (binary) %% 2 == 0)
sign [(1 + sums) %% 2 == 0] = -1L
else
sign [sums %% 2 == 0] = -1L
sign
}

probmv = function (sf, a, b)
{	if (inherits (sf, "ccdfmv") )
{	if (all (sf %\$% ".low") )
.comb.prob (sf, a, b)
else
stop ("needs cdf based on lower tails")
}
else
stop ("needs ccdfmv object")
}

.comb.prob = function (F, a, b)
{   has.matrix.args = (is.matrix (a) && is.matrix (b) )
if (has.matrix.args)
{   m = ncol (a)
if (nrow (a) != nrow (b) )
stop ("nrow (a) needs to equal nrow (b)")
if (ncol (a) != ncol (b) )
stop ("ncol (a) needs to equal ncol (b)")
}
else
{   m = length (a)
if (length (a) != length (b) )
stop ("length (a) needs to equal length (b)")
}
nF = as.integer (round (2 ^ m) )
binary = .binary.seq (m, nF)
sign = .binary.sign (binary)
y = 0
for (i in 1:nF)
{   x = a
j = as.logical (binary [i,])
if (has.matrix.args)
x [,j] = b [,j]
else
x [j] = b [j]
y = y + sign [i] * F (x)
}
y
}

pwith = function (...)
UseMethod ("pwith")

pwith.cksuv = function (sf, xlim = c (a, b), ..., a=-Inf, b=Inf)
{	. = attributes (sf)
xlim = .val.XLIM.uv (xlim)
data = .\$data
.pwith.eval.uv (.\$bw, .\$kernel@F, data\$n, data\$x, xlim [1], xlim [2],
isw=.\$is.weighted, data\$w)
}

pwith.cksmv = function (sf, xlim = cbind (a, b), ..., a=-Inf, b=Inf)
{	. = attributes (sf)
xlim = .val.XLIM.mv (xlim, .\$m)
data = .\$data
.pwith.eval.mv (.\$bw, .\$kernel@F, data\$n, .\$m, data\$x, xlim [,1], xlim [,2], isw=.\$is.weighted, data\$w)
}

psv = function (sf)
{	x = attr (sf, "data")\$x
sf (x)
}

ph4.pcomp2 = function (fh, gh, sqrt.mse=TRUE, aggregate=TRUE, dfh = psv (fh), dgh = psv (gh) )
{	nf = length (dfh)
ng = length (dgh)
xf = attr (fh, "data")\$x
xg = attr (gh, "data")\$x
kf = sum ( (dfh - gh (xf) )^2) / nf #f data, evaluated with f and g
kg = sum ( (dgh - fh (xg) )^2) / ng #g data, evaluated with f and g
if (sqrt.mse)
{	kf = sqrt (kf)
kg = sqrt (kg)
}
if (aggregate)
(kf + kg) / 2
else
c (kf, kg)
}

pdist = function (sf, ..., sqrt.mse=TRUE)
{	if (! is.list (sf) )
stop ("list required")
nfh = length (sf)
fv0 = vector ("list", nfh)
dists = matrix (0, nfh, nfh)
for (i in seq_len (nfh) )
{	if (! is.pdf (sf [[i]]) )
stop ("sf should only contain (probhat) density functions")
fv0 [[i]] = psv (sf [[i]])
}
if (nfh > 1)
{	for (i in 1:(nfh - 1) )
{	for (j in (i + 1):nfh)
{	dists [i, j] = dists [j, i] =
ph4.pcomp2 (sf [[i]], sf [[j]], sqrt.mse, TRUE, fv0 [[i]], fv0 [[j]])
}
}
}
if (inherits (sf, "ph4.gset") )
rownames (dists) = colnames (dists) = attr (sf, "levnames")
dists
}

ph4.rdist = function (d, n)
{	if (missing (n) )
trim = FALSE
else
{	trim = TRUE
ntop = n
}

n = nrow (d)
m = ncol (d)
if (n != m || n < 2)
stop ("square matrix, with at least 2 rows required")
names = rownames (d)
if (is.null (names) )
names = paste0 ("v", 1:n)

n2 = (n - 1) * (n) / 2
if (trim && ntop > n2)
stop (sprintf ("n=%s (top-values), but only %s upper-right triangle values", ntop, n2) )
names2a = names2b = character (n2)
d2 = numeric (n2)

k = 1
for (i in 1:(n - 1) )
{	for (j in (i + 1):n)
{	names2a [k] = names [i]
names2b [k] = names [j]
d2 [k] = d [i, j]
k = k + 1
}
}

I = order (d2)
df = data.frame (a=names2a, b=names2b, d=d2)[I,]
df = data.frame (comb = paste (names2a, names2b, sep=":"), dist=d2)[I,]
rownames (df) = 1:n2
if (trim)
df [1:ntop,]
else
df
}
```

## Try the probhat package in your browser

Any scripts or data that you put into this service are public.

probhat documentation built on May 12, 2021, 5:08 p.m.