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,
