Commit ed0bbd34 by Eric Coissac

Make the resolution of the linear system using the QR decomposition like in the lm function

parent 31958291
......@@ -77,11 +77,25 @@ mcor = function(...) {
return(cov/vv)
}
pm.fit = function(x, y,
offset = offset,
pm.fit = function(covmat,y,xs,
singular.ok = singular.ok,
...) {
tol) {
xy.cov = covmat[xs,y]
xx.cov = covmat[xs,xs,drop=FALSE]
xx.qr = qr(xx.cov,tol=tol)
if (!singular.ok && xx.qr$rank < ncol(xx.cov))
stop("singular fit encountered")
coef = qr.coef(xx.qr,xy.cov)
return(list(coefficients=coef,
rank=xx.qr$rank + 1
)
)
}
.ctrace <- function(MAT) sum(MAT^2)
......@@ -166,55 +180,26 @@ pm = function (formula,data, subset, weights, na.action, method = "qr",
sw = sqrt(w)
vars.norm = lapply(vars.norm, function(v) v * sw)
subset.w=sw > 0
}
data.cov = mvar(vars.norm)
xy.cov = data.cov[seq_len(nvars)!=irep,irep,drop=FALSE]
xx.cov = data.cov[seq_len(nvars)!=irep,seq_len(nvars)!=irep]
if (length(xx.cov)>1) {
ixx.cov=solve(xx.cov)
vp = diag(ixx.cov)
vp.max = max(vp)
if (vp.max > 0)
rank = (vp > (vp.max * tol)) + 1
}
else {
ixx.cov=1/xx.cov
rank=2
vars.norm = vars.norm[subset.w,]
}
if (!singular.ok && rank < ncol(vars.norm))
stop("singular fit encountered")
data.cov = mvar(vars.norm)
scales=ixx.cov %*% xy.cov
dim(scales)=NULL
names(scales)=attr(mt, "term.labels")
z = pm.fit(data.cov,
y = irep, xs = seq_len(nvars)!=irep,
singular.ok = singular.ok,
tol = tol)
if (scale)
anova = scales^2
z$varpart = z$coefficients^2
else {
data.cor = mcor(vars.norm)
xy.cor = data.cor[seq_len(nvars)!=irep,irep,drop=FALSE]
xx.cor = data.cor[seq_len(nvars)!=irep,seq_len(nvars)!=irep]
if (length(xx.cor)>1) {
ixx.cor=solve(xx.cor)
vp = diag(ixx.cor)
vp.max = max(vp)
if (vp.max > 0)
rank = (vp > (vp.max * tol)) + 1
}
else {
ixx.cor=1/xx.cor
rank=2
}
anova=ixx.cor %*% xy.cor
dim(anova)=NULL
names(anova)=attr(mt, "term.labels")
anova = anova^2
anova = pm.fit(data.cor,
y = irep, xs = seq_len(nvars)!=irep,
singular.ok = singular.ok,
tol = tol)$coefficients
z$varpart = anova^2
}
Y = vars.norm[[irep]]
......@@ -237,6 +222,8 @@ pm = function (formula,data, subset, weights, na.action, method = "qr",
# Intercept estimation
scales = z$coefficients
b = Ymean - Reduce(function(a,b) a+b,
mapply(function(d,coef,rot) coef * (d %*% rot),
Xmeans, scales, A_xys,
......@@ -257,15 +244,12 @@ pm = function (formula,data, subset, weights, na.action, method = "qr",
yhat = sweep(yhat,2,b,"+")
z <- list(coefficients = scales,
varpart = anova,
residuals = Y - yhat,
fitted.values = yhat,
weights = w,
rank = rank,
df.residual = if (!is.null(w))
sum(w != 0)
else nrow(Y))
z$residuals = Y - yhat
z$fitted.values = yhat
z$weights = w
z$df.residual = if (!is.null(w))
sum(w != 0)
else nrow(Y)
if (ret.x)
z$x <- vars.norm[-irep]
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment