[Rd] How does .Fortran "dqrls" work?

From: yangleicq <yanglei_cq_at_126.com>
Date: Wed, 25 Apr 2012 23:00:21 -0700 (PDT)


ink1">Hi, all.
I want to write some functions like glm() so i studied it.
In glm.fit(), it calls a fortran subroutine named "dqrfit" to compute least ink2">squares solutions
 to the system

              x * b = y

To learn how "dqrfit" works, I just follow how glm() calls "dqrfit" by my own example, my codes are given below:

> qr <-
> matrix(c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14,2.3,1.7,1.3,1.7,1.7,1.6,1,1.7,1.7,1.7),ncol=2)
> qr

      [,1] [,2]

 [1,] 4.17  2.3
 [2,] 5.58  1.7
 [3,] 5.18  1.3
 [4,] 6.11  1.7
 [5,] 4.50  1.7
 [6,] 4.61  1.6
 [7,] 5.17  1.0
 [8,] 4.53  1.7
 [9,] 5.33  1.7

[10,] 5.14 1.7
> n=10
> p=2
> y <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
> ny=1L
> tol=1e-07
> coefficients=double(p)
> residuals=double(n)
> effects=double(n)
> rank=integer(1L)
> pivot=1:n
> qraux=double(n)
> work=double(2*n)
>
>
>

> fittt<-.Fortran("dqrls", qr =qr, n = n,
+                 p = p, y = y, ny = ny, tol = tol,
coefficients=coefficients,
+                 residuals = residuals, effects = effects, 
+                 rank = rank, pivot = pivot, qraux = qraux, 
+                 work = work, PACKAGE = "base")

>
> fittt$coefficients

[1] 0 0

but when i use lm() which also calls "dqrls" internally to fit this model, it gives reasonable result.

> lm(y~qr)

Call:
lm(formula = y ~ qr)

Coefficients:

(Intercept)          qr1          qr2  
    11.1766      -0.8833      -1.2628  


when I change the coefficients to be c(1,1), the output from "dqrls", fittt$coefficients also equals to c(1,1). That means the .Fortran("dqrls", qr=qr,n=n,p=p,...) did nothing to the coefficients! I don't know why, is there anything I did wrong or missed? How can I get the result from "dqrls" as what lm() or glm() gets from "dqrls"?

Thanks in advance. Best Regards.

--
View this message in context: http://r.789695.n4.nabble.com/How-does-Fortran-dqrls-work-tp4588973p4588973.html
Sent from the R devel mailing list archive at Nabble.com.

______________________________________________
R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Thu 26 Apr 2012 - 11:19:54 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

All messages

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 27 Apr 2012 - 13:00:51 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive