[R] Integrating a function

[David Winsemius]

Does this do any better, or perhaps offer other ideas? I could not 
understand what the "Hct" and step()'s were doing from the description 
of your desired result, so they were omitted from geq(). I also could 
not understand why you were rounding off the integrals or why there wee 
two integral terms. I think the problem may be (in part) that 
integrate() returns a list unless you specify the value column. When I 
tested your function or my earliest efforts at fixing it I continued to 
get a single value despite passing it a sequence or a matrix.

surviv<-function(x){k=.001;alpha=0.3;tau=70; 
exp(-k*x)/(1+exp(alpha*(x-tau)))}
geq<-function(n){ integrate(surviv, 0,n)$value }

Nmax=250
gxt<-matrix(nrow=Nmax)
gxt<-as.matrix(1:250)
# there may be better methods o making this matrix

geq.results<-geq.results<-apply(gxt,1,geq)
geq.results
plot(gxt, geq.results)

I suspect that I have not understood what you wanted from geq(.) since 
the plot looks rather "boring". I am guessing that Hct0 is a hematocrit 
value. And geq basically returns a value very close to its input up 
until tau and then stays at approximately tau for the rest of time.

-- 
David Winsemius, MD

Robert Kalicki wrote:
> Hello everybody!
>
> I have problems with integrating my function.
>
> My primary function is a “survival function” of the following type:
>
> surviv <- exp(-k*x)/(1+exp(alpha*(x-tau)))
>
> I would like to integrate this function over a defined range and obtain a
> vector with all the values from integrate(surviv, 0, x) for x <- 1:N. 
>
> I unfortunately obtain just a scalar value corresponding to the last
> integrated x.
>
> How do I have to proceed if I want to obtain a vector with all the results?
>
> I have tried to solve this problem by adding a “while loop”. It works but
> the expression becomes to complex if I want to perform further optimization.
>
> Have any one a solution for this problem?
>
> Thanks a lot in advance.
>
>  
>
> Example
>
> # General equation
>
> beta1 <- 0.001
> beta2 <- 0
> Hct0 <- 0.27
> tau <- 70
> k <- 0.001
> alpha <- 0.3
> T2 <-40
>
> step <- function(x){
>    if(x>=0) step <- 1
>    else step <- 0
>
> }
>    Nmax <- 250
>    n <- 1
>       while(n<Nmax){
>          surviv <- function(n){
>             surviv <- exp(-k*n)/(1+exp(alpha*(n-tau)))
>          }
>          geq <- function(n){
>             geq <-
> Hct0+beta1*as.integer(integrate(surviv,0,n)[1])+step(n-T2)*(beta2-beta1)*as.
> integer(integrate(surviv,0,n-T2)[1])
>          }
>             plot(n, geq(n), type="p", pch=19, xlim=c(0,250),
> ylim=c(0.25,0.35), xlab="time in days", ylab="Hct")
>
>             par(new=T)
>             print(c(n, surviv(n), geq(n),
> as.integer(integrate(surviv,0,n)[1]))                
>        n <- n + 1
>        }
>   
> Robert M. Kalicki, MD
>
> Department of Nephrology and Hypertension