Commit eff483d7 authored by Martin Maechler's avatar Martin Maechler

Vectorized findSig() and findMu(); findSig.1() now works

parent d5fedbe9
......@@ -21,69 +21,46 @@ findSig <- function(Mu0, result, covar, equal.tol = 1e-6) { # equal.tol > 1e-16
}
## This does _not_ work (at least in some cases)
findSig.1 <- function(Mu0, result, covar, equal.tol){ # equal.tol > 1e-16
##' Now simplified and very close to findSig0() -- just returning weights additionally
findSig.1 <- function(Mu0, result, covar, equal.tol) { # equal.tol > 1e-16
stopifnot(length(Mu0) == 1L)
ms.w <- result$MS_weight
nd <- order(ms.w[, "Mu"])
n <- length(nd)
## revert order for all (w, mu, sig):
ms.w <- ms.w[nd, ]
w <- result$weights_set[, nd]
ind.uni <- c(TRUE, sapply(1:(n-1), function(i) !
isTRUE(all.equal(ms.w[i, ], ms.w[i+1, ], tol = equal.tol))))
p <- rep(1, n) # pointer, different weights of same ms have same pointer
k <- 1
for(i in 1:n){
if(ind.uni[i]){
p[i] <- k
k <- k + 1
}else{
p[i] <- p[i-1]
}
}
mu.w <- ms.w[ind.uni, "Mu"]
sig.w <- ms.w[ind.uni, "Sig"]
if(Mu0 < min(mu.w) || Mu0 > max(mu.w)) stop(sprintf("Mu0 must be in [%g, %g]",
min(mu.w), max(mu.w)))
mu.w <- ms.w[, "Mu"]
sig.w <- ms.w[, "Sig"]
if(Mu0 < min(mu.w) || Mu0 > max(mu.w))
stop(sprintf("Mu0 must be in [%g, %g]", min(mu.w), max(mu.w)))
i <- findInterval(Mu0, mu.w) # [..)...[..]
if(i == n) { # Mu0 is max(mu.w)
if(i == n || # Mu0 is max(mu.w)
isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))) {
list(Sig = sig.w[i], weight = w[, i])
} else{
m1 <- which(p == i)
m2 <- which(p == i+1)
} else {
a <- (Mu0 - mu.w[i+1])/(mu.w[i] - mu.w[i+1])
# solve for a : Mu0 = a* Mu1 + (1-a)* Mu2 :
wm <- matrix(0, n, m1*m2)
for(s in m1) {
wm[,(s-1)*m2 + 1:m2] <- a*w[, s] + (1-a)*w[, m2]
}
list(Sig = sqrt(colSums(wm *(covar %*% wm))), #return Sig0
weight = wm) # might have replicate columns
wm <- a*w[, i] + (1-a)*w[, i+1]
list(Sig = sqrt(colSums(wm *(covar %*% wm) )),
weight = wm)
}
}
##' This is much cleaner and logical than findSig.1 above ==> use it everywhere
##' This is much cleaner and logical than the original findSig() ==> use a version of it
##' A version of findSig() only used inside findMu():
findSig0 <- function(Mu0, result, covar, equal.tol = 1e-8){ # equal.tol > 1e-16
stopifnot(length(Mu0) == 1L)
ms.w <- result$MS_weight
nd <- order(ms.w[, "Mu"])
n <- length(nd)
## revert order for all (w, mu, sig):
ms.w <- ms.w[nd, ]
w <- result$weights_set[, nd]
mu.w <- ms.w[ , "Mu"]
sig.w <- ms.w[ , "Sig"]
mu.w <- ms.w[, "Mu"]
sig.w <- ms.w[, "Sig"]
i <- findInterval(Mu0, mu.w) # [..)...[..]
if(i == n || # Mu0 is max(mu.w)
isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))) {
if(i == n || isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))) {
sig.w[i]
} else {
a <- (Mu0 - mu.w[i+1])/(mu.w[i] - mu.w[i+1])
......
......@@ -43,7 +43,7 @@ cov21 <- muS.sp500$covar[iS, iS]
CLsp.21 <- CLA(muS.sp500$mu[iS], covar=cov21, lB=0, uB = 1/2)
CLsp.21 # 14 turning points
summary(tpM <- CLsp.21$MS_weights[,"Mu"])
str(m0 <- c(min(tpM),seq(0.00205, 0.00525, by = 0.00005))
str(m0 <- c(min(tpM),seq(0.00205, 0.00525, by = 0.00005), max(tpM)))
sig. <- findSig(m0, result=CLsp.21, covar=cov21)
str(sig.)
stopifnot(dim(sig.$weight) == c(21, length(m0)))
......
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