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