Full_Name: Stephen Weigand
Version: 1.9.0
OS: Mac OS X 10.3.4
Submission from: (NULL) (68.115.89.235)
When running an iteratively reweighted least squares program R crashes and the
following is
written to the console.app (when using R GUI) or to stdout (when using R from
the command
line):
Parameter 5 to routine DSTEIN was incorrect
Mac OS BLAS parameter error in DSTEIN, parameter #0, (unavailable), is 0
In case it helps, here's the function and function call that causes the crash.
n <- 4
p <- 1
f <- 2; k <- 2
lms.bcn.univar <- function(y, B, p, f, k, lambda.0, mu.0, sigma.0){
n <- length(y)
D11 <- matrix(1, nrow = p, ncol = n)
D1 <- rbind(cbind(D11, matrix(0, nrow = p, ncol = f)),
cbind(matrix(0, nrow = f, ncol = n), diag(f)))
D22 <- rbind(rep(1:0,n),
rep(0:1,n))
x1 <- t(D22)[,1]
x2 <- t(D22)[,2]
D2 <- rbind(cbind(diag(n), matrix(0, nrow = n, ncol = 2*n)),
cbind(matrix(0, nrow = f, ncol = n), D22))
R <- matrix(0, nrow = n, ncol = k * n)
Ra <- diag(n + n*k)
A11 <- function(lambda, mu, sigma){
diag((1 + 2 * lambda^2 * sigma^2)/(mu^2 * sigma^2), nrow = n)
}
A22 <- function(lambda, sigma, n) {
A <- matrix(0, nrow = n*2, ncol = n*2)
block <- c(7 * sigma^2 / 4, lambda * sigma,
lambda * sigma, 2 / (sigma^2))
## code to get the indices of the elements of a block
## diagonal matrix when the matrix is treated as a vector.
m <- n*2
s <- integer(0)
for (i in seq(1, m-1, by = 2)) s <- c(s, rep(i:(i+1),2))
block.indices <- m * rep(0:(m-1), each = 2) + s
A[block.indices] <- block
return(A)
}
A12 <- function(lambda, mu, sigma, n) {
A <- matrix(0, nrow = n, ncol = n * 2)
block <- c(-1/(2 * mu), (2*lambda)/(mu * sigma))
block.indices <- n * 0:(2*n-1) + rep(1:n, each = 2)
A[block.indices] <- block
return(A)
}
I.exp <- function(A11,A12,A22) {
rbind(cbind(A11, A12),
cbind(t(A12),A22))
}
G <- function(D2, Ra, A11, A12, A22){
D2 %*% Ra %*% I.exp(A11,A12,A22) %*% t(Ra) %*% t(D2)
}
W <- function(G){
G11 <- G[1:n, 1:4]
G12 <- G[1:4, 5:6]
G22 <- G[5:6, 5:6]
}
W <- function(G,n,k){
G11 <- G[1:n, 1:n]
G12 <- G[1:n, (n+1):(n+k)]
G22 <- G[(n+1):(n+k), (n+1):(n+k)]
G11 - G12 %*% solve(G22) %*% t(G12)
}
zbc <- function(y,lambda, mu, sigma) {
((y/mu)^lambda - 1) / (lambda * sigma)
}
L1 <- function(y,lambda, mu, sigma) {
z <- zbc(y,lambda, mu, sigma)
return(z/(mu * sigma) + lambda * (z^2 - 1) / mu)
}
L2 <- function(y,lambda,mu,sigma){
z <- zbc(y,lambda, mu, sigma)
L.lambda <- z/lambda * (z - log(y/mu)/sigma) - log(y/mu) * (z^2 - 1)
L.sigma <- (z^2 - 1)/sigma
L <- c(L.lambda, L.sigma)
return(L[c(1, n + 1) + rep(0:(n-1), each = 2)])
}
u1 <- function(L1, L2, G, R, D22) {
G12 <- G[1:n, (n+1):(n+k)]
G22 <- G[(n+1):(n+k), (n+1):(n+k)]
return(L1 + (R - G12 %*% solve(G22) %*% D22) %*% L2)
}
u2 <- function(L2, A12, A22, R, D11, beta.new, beta.old){
L2 - (t(A12) + A22 %*% t(R)) %*% t(D11) %*% (beta.new - beta.old)
}
loglike <- function(y, l, m, s) {
sum( l * log(y/m) - log(s) - 0.5 * zbc(y,l,m,s)^2 )
}
loglikeOptim <- function(par,y){
lambda <- par[1]
mu <- par[2]
sigma <- par[3]
-1 * loglike(y, lambda, mu, sigma)
}
ll.0 <- loglike(y, lambda.0, mu.0, sigma.0)
require(MASS)
lambda.old <- lambda.0; mu.old <- mu.0; sigma.old <- sigma.0
parameters <- as.data.frame(matrix(NA, ncol = 4, nrow = B))
colnames(parameters) <- c("lambda", "mu", "sigma", "loglike")
for(i in 1:B) {
##cat(paste(i, ","))
a11 <- A11(lambda.old, mu.old, sigma.old)
a22 <- A22(lambda.old, sigma.old, n)
a12 <- A12(lambda.old, mu.old, sigma.old, n)
g <- G(D2, Ra, a11, a12, a22); w <- W(g, n ,2)
l1 <- L1(y, lambda.old, mu.old, sigma.old)
l2 <- L2(y, lambda.old, mu.old, sigma.old)
Y.mu <- solve(w) %*% u1(l1, l2, g, R, D22) + t(D11) %*% mu.old
mu.new <- coef(lm.gls(Y.mu ~ 1, W = w))
psi.old <- rbind(lambda.old, sigma.old)
Y.psi <- solve(a22) %*% u2(l2, a12, a22, R, D11, mu.new, mu.old) + t(D22)
%*% psi.old
psi.new <- coef(lm.gls(Y.psi ~ -1 + x1 + x2, W = a22))
lambda.new <- psi.new[1]; sigma.new <- psi.new[2]
parameters[i, 1] <- lambda.new
parameters[i, 2] <- mu.new
parameters[i, 3] <- sigma.new
parameters[i, 4] <- loglike(y, lambda.new, mu.new, sigma.new)
### update the old with the new
lambda.old <- lambda.new; mu.old <- mu.new; sigma.old <- sigma.new;
}
return(list(lambda.0 = lambda.0, mu.0 = mu.0, sigma.0 = sigma.0, loglike.0 =
ll.0,
parameters = parameters))
}
require(MASS)
set.seed(1676)
y <- exp(rnorm(20) + 2)
y <- round(y,2)
##print(lms.bcn.univar(y, B=15, p=1, f=2, k=2, 0.3, 8, 0.8))
y <- rgamma(100,3,2)
print(lms.bcn.univar(y, B=5, p=1, f=2, k=2, 0.3,
mu.0 = median(y),
sigma.0 = sd(y)/median(y)))
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-devel
