[R] Is it possible to recursively update a function?
Uwe Ligges
ligges at statistik.tu-dortmund.de
Sat Mar 6 18:45:00 CET 2010
It is not obvious to me if you can calculate f1, ..., f4, ...
automatically or have to set them manually. Looks like you need to do it
manually. In that case I'd suggest to write something along the lines:
dose <- c(-1.47, -1.1, -0.69, -0.42, 0.0, 0.42)
f1 <- function(x) exp(-x)
f2 <- function(x) f1(x) * (1 - 0.2^x)
f3 <- function(x) f2(x) * (1 - 0.3^x)
f4 <- function(x) f3(x) * 0.3^x
calcid <- function(f, dose){
denom <- integrate(f, 0, 100)$value
alpha <- integrate(function(x) x * f(x), 0, 100)$value / denom
which.min(abs(-0.5 * log(0.2^(-1/alpha) - 1) - dose))
}
calcid(f2, dose)
calcid(f3, dose)
calcid(f4, dose)
If the hierarchy becomes deeper (i.e. 10 functions or so, it might make
sense to reduce the number of hierarchical calls for efficiency, but
that is not necessary for the current setup.
Uwe Ligges
On 05.03.2010 17:34, Ming Zhong wrote:
> I was trying to replicate one CRM simulation. The following code works but
> seems redundant so I want to create a loop.
>
> ####O'Quigley CRM example 1#####
> f1<-function(x) exp(-x)
> dose<-c(-1.47,-1.1,-.69,-.42,0.0,.42)
> p<-c(0.05,0.1,0.2,0.3,0.5,0.7)
> y<-c(0,0,1,0,1,1,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,1,0,1,1)
>
> f2<-function(x) f1(x)*(1-0.2^x)
> denom2<-integrate(f2,0,100)$value
> alpha2<-integrate(function(x) x*f2(x),0,100)$value/denom2
> id<-which.min(abs(-.5*log(0.2^(-1/alpha2)-1)-dose)) ### id is used to
> locate the next prob.
>
> f3<-function(x) f2(x)*(1-.3^x)
> denom3<-integrate(f3,0,100)$value
> alpha3<-integrate(function(x) x*f3(x),0,100)$value/denom3
> id<-which.min(abs(-.5*log(0.2^(-1/alpha3)-1)-dose))
>
> f4<-function(x) f3(x)*.3^x
> denom4<-integrate(f4,0,100)$value
> alpha4<-integrate(function(x) x*f4(x),0,100)$value/denom4
> id<-which.min(abs(-.5*log(0.2^(-1/alpha4)-1)-dose))
>
> 2010/3/5 Uwe Ligges<ligges at statistik.tu-dortmund.de>
>
>>
>>
>> On 05.03.2010 01:40, Carl Witthoft wrote:
>>
>>> My foolish move for this week: I'm going to go way out on a limb and
>>> guess what the OP wanted was something like this.
>>>
>>> i=1, foo = x*exp(-x)
>>>
>>> i=2, foo= x^2*exp(-x)
>>> i=3, foo = x^3*exp(-x)
>>> .
>>> .
>>> .
>>>
>>>
>>> In which case he really should create a vector bar<-rep(na,5) ,
>>> and then inside the loop,
>>>
>>> bar[i]<-x^i*foo(x)
>>>
>>
>> Since in this case foo(x) is independent of i, you are wasting resources.
>> Moreover you could calculate it for a whole matrix at once. Say you want to
>> calculate this for i=1, ..., n with n=5 for some (here pseudo random x),
>> then you could do it simpler after defining some data as in:
>>
>> set.seed(123)
>> x<- rnorm(10)
>> n<- 5
>>
>>
>> using the single and probably most efficient line:
>>
>> outer(x, 1:n, "^") * exp(-x)
>>
>> or if x is a length 1 vector then even simpler:
>>
>> set.seed(123)
>> x<- rnorm(1)
>> n<- 5
>>
>> x^(1:5) * exp(-x)
>>
>> But we still do not know if this is really the question ...
>>
>> Uwe Ligges
>>
>>
>>
>>
>>> Carl
>>>
>>>
>>>
>>> quoted material:
>>> Date: Thu, 04 Mar 2010 11:37:23 -0800 (PST)
>>>
>>>
>>> I need to update posterior dist function upon the coming results and
>>> find the posterior mean each time.
>>>
>>>
>>> On Mar 4, 1:31 pm, jim holtman<jholt..._at_gmail.com> wrote:
>>> > What exactly are you trying to do? 'foo' calls 'foo' calls 'foo' ....
>>> > How did you expect it to stop the recursive calls?
>>> >
>>> >
>>> >
>>> >
>>> >
>>> > On Thu, Mar 4, 2010 at 2:08 PM, Seeker<zhongm..._at_gmail.com> wrote:
>>> > > Here is the test code.
>>> >
>>> > > foo<-function(x) exp(-x)
>>> > > for (i in 1:5)
>>> > > {
>>> > > foo<-function(x) foo(x)*x
>>> > > foo(2)
>>> > > }
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html<http://www.r-project.org/posting-guide.html>
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>
>
More information about the R-help
mailing list