File: negexp.R

package info (click to toggle)
r-cran-mass 7.3-51.1-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster
  • size: 2,148 kB
  • sloc: ansic: 664; makefile: 2
file content (40 lines) | stat: -rw-r--r-- 1,415 bytes parent folder | download | duplicates (6)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
# file MASS/R/negexp.R
# copyright (C) 1994-9 W. N. Venables and B. D. Ripley
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 or 3 of the License
#  (at your option).
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  http://www.r-project.org/Licenses/
#

negexp.SSival <- function(mCall, data, LHS)
{
    x <- eval(mCall[["x"]], data)
    if(length(x) < 3L)
        stop("at least 3 distinct 'x' values are needed")
    y <- eval(LHS, data)
    mx <- mean(x)
    b <- as.vector(lsfit(cbind(x - mx,  - (x - mx)^2/2), y)$coef)
    rx <- range(x)
    xh <- mx + b[2L]/b[3L]
    if(prod(xh - rx) < 0)
        if(xh - rx[1L] > rx[2L] - xh)
            rx[2L] <- xh
        else rx[1L] <- xh
    x0 <- c(rx[1L], sum(rx)/2, rx[2L])
    dy <- diff(b[1L] + b[2L] * (x0 - mx) - (b[3L] * (x0 - mx)^2)/2)
    th <- (x0[2L] - x0[1L])/log(dy[1L]/dy[2L])
    b <- as.vector(lsfit(exp( - x/th), y)$coef)
    pars <- list(b[1L], b[2L], th)
    names(pars) <- mCall[c("b0", "b1", "th")]
    print(unlist(pars))
    pars
}