On Sat, Mar 28, 2009 at 2:33 AM, Jim Silverton <jim.silver...@gmail.com> wrote:
> Can someone explain why I am getting the following error: in the r code
> below?
Are you asking what the error message means or how you are getting a
computationally singular system? (My guess is that currvar is very
small but you would need to check with tracebacks, etc. If this is
inside a long-running function call then you should check
?dump.frames
to determine how to set up post-mortem debugging.)
However, you would do yourself a favor if you checked on the numerical
properties of the way you are performing your calculation. See, for
example, the "Comparisons" vignette in the Matrix package. The
timings described there may not be that important to you (it appears
that you are dealing with a 2 by 2 system - if so, why do you
hard-wire the dimension into this statement when you could use
diag(nrow(XX1)) instead?) but enhancing the numerical properties of
the calculation may be.
Naively translating a formula into an expression like the one below is
definitely suboptimal. I can't see what is going to be done with this
inverse because that expression doesn't occur in the code fragment you
show but just because we see a formula like
(X'X)^{-1} X'y
doesn't mean this is a good way to calculate the result. The first
rule of numerical linear algebra is that you never (well, hardly ever)
calculate the inverse of a matrix. To calculate the inverse of a
matrix you must use some kind of an decomposition like the QR or the
LU or the Cholesky. If you then look at your formula in terms of the
decomposition you usually find that you don't need to evaluate the
inverse at all.
Look through the code for lm (which calls ls.fit) or glm (which calls
glm.fit) and you will not see the R code equivalent of that formula.
You will see a QR decomposition of X or an expression like
chol(crossprod(X))
Note also that in the code shown below you evaluate residuals twice in
order to calculate the sum of squares of residuals. I don't know if
Y1 is a vector or a matrix. If it is a vector then the residual sum
of squares could be calculated as
sum((Y1-(t(XX1)%*%currbeta1)^2)
or, better,
sum((Y1 - crossprod(XX1, currbeta1))^2)
Even if Y1 is a matrix, an expression like
crossprod(Y1 - crossprod(XX1, currbeta1))
will get you the desired result more effectively.
So we can't tell you why you are getting the error message other than
"because the matrix you are trying to invert is computationally
singular?".
You will need to dig into the calculation to determine why this
particular calculation is failing and how to perform the calculation
more effectively. But that is what research is about, isn't it?
> Error in solve.default(diag(2) + ((1/currvar) * (XX1 %*% t(XX1)))) :
> system is computationally singular: reciprocal condition number = 0
> In addition: There were 50 or more warnings (use warnings() to see the first
> 50)
>
> The R code is part of a bigger program.
>
>
>
>
>
> ##sample from full conditional distribution of Si
> #Prob(Si = 1)
> for (j in 1:N)
> {
>
> numerat =
> currphi1*exp((-1/(2*currvar))*t(Y1-(t(XX1)%*%currbeta1))%*%(Y1-(t(XX1)%*%currbeta1)))
> denomin =
> currphi2*exp((-1/(2*currvar))*t(Y2-(t(XX2)%*%currbeta2))%*%(Y2-(t(XX2)%*%currbeta2)))
> sum=denomin + numerat
> ProbSi = numerat/sum
> arunofSi[j]=rbinom(1,1,ProbSi) #Generate 50 Bernoulli rvs and
> assign them to arunofSi array
>
> }
>
> N0=sum(arunofSi==0) #We check the number of zeros in the
> arrays
> N1= sum(arunofSi==1) #We check the number of ones in the array
> #The N0 and N1 values are the number of
> subjects in groups 0 and 1.
> #This is fed into the Dirichlet function
> below to create the currphi's
>
>
>
>
> Jim
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
