This shows you the differences between two versions of the page.
Both sides previous revision Previous revision | Last revision Both sides next revision | ||
the_gsru_pseudo-code [2012/05/14 11:16] andres |
the_gsru_pseudo-code [2012/05/14 11:19] andres |
||
---|---|---|---|
Line 1: | Line 1: | ||
This is from [[http://www.journalofdairyscience.org/article/S0022-0302(08)71471-1/abstract | Legarra & Misztal 2008]] | This is from [[http://www.journalofdairyscience.org/article/S0022-0302(08)71471-1/abstract | Legarra & Misztal 2008]] | ||
+ | * Fortran code: | ||
<code fortran> | <code fortran> | ||
Double precision:: xpx(neq),y(ndata),e(ndata),X(ndata,neq), & | Double precision:: xpx(neq),y(ndata),e(ndata),X(ndata,neq), & | ||
Line 25: | Line 26: | ||
enddo | enddo | ||
enddo | enddo | ||
+ | </code> | ||
+ | |||
+ | * R code: | ||
+ | <code rsplus> | ||
+ | #get diagonal of X’X | ||
+ | for (i in 1:neq) { | ||
+ | xpx[i]=crossprod(X[,i],X[,i]) | ||
+ | } | ||
+ | for (iter in 1:1000) { | ||
+ | #Gauss Seidel | ||
+ | for (i in 1:neq){ | ||
+ | lhs=xpx[i]+lambda | ||
+ | rhs=crossprod(X[,i],e) + xpx[i]*ahat[i] | ||
+ | val=rhs/lhs | ||
+ | e = e - X[,i]*(val - ahat[i]) | ||
+ | ahat[i]=val | ||
+ | } | ||
+ | } | ||
</code> | </code> |