# R/GPois.R In Hanchao-Zhang/ProDenICA: Product Density Estimation for ICA using tilted Gaussian density estimates

#### Documented in GPois

```GPois <-
function(x, df = 6, B = 500, order = 1, widen = 1.2, density.return=FALSE, ...)
{
### x must be a unit vector with mean zero
### Note; if x has outliers, this will make the gaussian more robust
### order and widen are robustness parameters in computing the grid
### This function also computes derivitives
x <- drop(scale(x))
n <- length(x)
rangex <- range(x)
if(order == 1)
rx <- rangex
else {
rx <- sort(x)[c(order, n - order + 1)]
}
rx <- ylim.scale(rx, diff(rx) * widen)
xg <- seq(from = rx[1], to = rx[2], length = B)
gaps <- diff(rx)/(B - 1)
xcuts <- c(min(rangex[1], rx[1]) - gaps/2, xg[ - B] + gaps/2, max(
rangex[2], rx[2]) + gaps/2)
ys <- as.vector(table(cut(x, xcuts)))
gxg <- dnorm(xg)
bigdata <- list(xg=xg, gxg=gxg, ys=ys)
##  assign("bigdata",bigdata,frame=0)
## assign("df", df, frame = 0)
pois.fit <- gam(ys ~ s(xg, df) + offset(logb(gxg)), family = poisson,
data = bigdata, ...)
## Now to get the derivitives
Gs <- predict(pois.fit)-logb(gxg)
if(density.return){
## package up the function to return for plotting
density=list(x=xg,y=exp(Gs+logb(gxg)))
}
Gs <- Gs + logb(sum(gxg)/sum(fitted(pois.fit)))
resp <- Gs + residuals(pois.fit, type = "working")
weights <- pois.fit\$weights
df <- B - pois.fit\$df.residual
pois.refit <- smooth.spline(xg, resp, weights, df)
Gs <- predict(pois.refit, x, deriv = 0)\$y
gs <- predict(pois.refit, x, deriv = 1)\$y
gps <- predict(pois.refit, x, deriv = 2)\$y
rl=list(Gs = Gs, gs = gs, gps = gps)
if(density.return)rl\$density=density
rl
}
```
Hanchao-Zhang/ProDenICA documentation built on Jan. 17, 2022, 12:23 a.m.