
R : Copyright 2004, The R Foundation for Statistical Computing
Version 2.0.1  (2004-11-15), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> dat <- read.table("/tmp/tlm.dat", header = TRUE)
> dat
             r N.fix
1  -0.52407085    76
2   0.10536052    45
3  -0.17435339    50
4   0.19415601    42
5   0.48701498    51
6  -0.50681760    83
7  -0.17435339    50
8   0.55278982    42
9   0.15219182    73
10  0.49899117    85
11  0.10821358   140
12 -0.83034830   156
13 -0.30748470    68
14 -0.22314355    50
15  0.04879016    40
16 -0.04879016    42
17  0.75377180    40
18 -0.12516314    85
19 -0.36624439    75
> plot(r ~ N.fix, dat)
> tlm <- deriv( ~ 1-((N/k)^theta), c("k", "theta"), function(N, k, theta){})
> tlm
function (N, k, theta) 
{
    .expr1 <- N/k
    .expr2 <- .expr1^theta
    .value <- 1 - .expr2
    .grad <- array(0, c(length(.value), 2), list(NULL, c("k", 
        "theta")))
    .grad[, "k"] <- .expr1^(theta - 1) * (theta * (N/k^2))
    .grad[, "theta"] <- -(.expr2 * log(.expr1))
    attr(.value, "gradient") <- .grad
    .value
}
> fm1 <- nls(r ~ tlm(N.fix, k, theta), dat, start=c(k = 65, theta = 1), alg = "plinear", trace = TRUE)
2.347334 : 65.0000000  1.0000000  0.3836488 
2.233973 : 56.633165 -0.431863 -1.282247 
2.187513 : 54.8608041 -1.3747423 -0.4908351 
2.167521 : 52.8247351 -2.1635728 -0.3592099 
2.155828 : 51.4475043 -2.8951236 -0.3003772 
2.149315 : 50.4174840 -3.5536098 -0.2682023 
2.146427 : 49.7548459 -4.0636580 -0.2500007 
2.145514 : 49.3861622 -4.3828441 -0.2405601 
2.145301 : 49.2043941 -4.5476404 -0.2362051 
2.145260 : 49.1221764 -4.6225967 -0.2343549 
2.145252 : 49.0870152 -4.6544862 -0.2335967 
2.145251 : 49.0724421 -4.6676430 -0.2332897 
> newdat <- list(N.fix = seq(40, 160))
> lines(newdat$N.fix, predict(fm1, newdat))
> q('no')

Process R finished at Wed Feb 23 11:20:23 2005
