[Rd] DSTEIN error (PR#7047)

From: <weigand.stephen_at_charter.net>
Date: Sat 03 Jul 2004 - 13:09:38 EST


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

______________________________________________
R-devel@stat.math.ethz.ch mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-devel Received on Sat Jul 03 13:13:51 2004

This archive was generated by hypermail 2.1.8 : Fri 18 Mar 2005 - 08:59:07 EST