[R] automated integration?

From: francogrex <francogrex_at_mail.com>
Date: Tue 23 Jan 2007 - 12:57:05 GMT

When I run the script below it works well it ouputs the result that I need, in this case were n=45 and e=5.02689 the result is: 6.669578



dens<-function(x){
n=45
e=5.02689
a1=0.0987 
b1=0.04261
a2=1.043
b2=1.222

p=0.121
if (n>200)  c((n=n/2),(e=e/2))
if (e>30)  c((n=n/2),(e=e/2))
lpp<-((lgamma(a1+n)))-(lgamma(a1)+lfactorial(n)+log((1+(e/b1))^a1)+log((1+(b1/e))^n))
lzz<-((lgamma(a2+n)))-(lgamma(a2)+lfactorial(n)+log((1+(e/b2))^a2)+log((1+(b2/e))^n))
pp<-exp(lpp)
zz<-exp(lzz)
qq= (p*pp)/((p*pp)+((1-p)*zz))
lgam<-(log((b1+e)^(a1+n)) + log(x^((a1+n)-1))) - (lgamma((a1+n)) + ((b1+e)*x))
lgom<-(log((b2+e)^(a2+n)) + log(x^((a2+n)-1))) - (lgamma((a2+n)) + ((b2+e)*x))
gam<-exp(lgam)
gom<-exp(lgom)
(qq*gam) + ((1-qq)*gom)
}

integ<-function(x){
n=45
e=5.02689
if (x>(n/e)) return (x=Inf)
if (x<0) return (x=0)
u<-integrate(dens,lower=0, upper=x)$value u-0.05
}

uni<-uniroot(integ,c(0,100), tol=1e-10)
uni$root

But what if n and e are vectors (like: n=c(31,22,47,38) and e=c(5.2,2.8,3.4,2.1)? can I have a string of results all at once instead of entering the values of n and e one by one it takes such a long time especially if the n and e contain a couple of 100 values!

Any suggestions are appreciated. Thanks

-- 
View this message in context: http://www.nabble.com/automated-integration--tf3063979.html#a8521241
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
R-help@stat.math.ethz.ch 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.
Received on Wed Jan 24 00:00:21 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 23 Jan 2007 - 13:30:28 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.