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<lig...@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@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.
______________________________________________
R-help@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
and provide commented, minimal, self-contained, reproducible code.
