[Rd] DSTEIN error (PR#7047)

[weigand.stephen at charter.net]

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)))