Hello! I know that probably my question is rather simple but I' m a very beginner R-user. I have to numerically integrate the product of two function A(x) and B(x). The integretion limits are [X*; +inf] Function A(x) is a pdf function while B(x)=e*x is a linear function whose value is equal to 0 when the x < X* Moreover I have to iterate this process for different value of X* and for different pdf of the same type. I know the comand INTEGRATE but I can' t make it work. Which is the best function to do this and how does it work? Thank you very much in advanced!! -- View this message in context: http://r.789695.n4.nabble.com/Numerical-integration-tp3634365p3634365.html Sent from the R help mailing list archive at Nabble.com.
Numerical integration
3 messages · nany23, Dennis Murphy
Hi:
You could write the function this way:
f <- function(x, xstar, k) dnorm(x) * k * x * (x >= xstar)
where the term in parentheses is a logical. For any x < xstar, f will
be zero by definition. Substitute your density in for dnorm().
To integrate over a grid of (xstar, k) values, you could try this:
# generate a grid of (xstar, k) pairs
pars <- expand.grid(xstar = seq(0, 2, by = 0.2), k = seq(0.5, 2, by = 0.5))
# We modify the function to do the integration and return its value.
# Note that since f = 0 for all x <= xstar, integrating from -Inf to Inf
# yields the same answer as integrating from xstar to Inf.
fun <- function(xstar, k) {
f <- function(x, xstar, k) dnorm(x) * k * x * (x >= xstar)
integrate(f, -Inf, Inf, xstar = xstar, k = k)$value
}
# Method 1: (outputs a data frame)
library(plyr)
out <- mdply(pars, fun)
# Method 2a: (outputs a matrix)
out <- cbind(pars, mapply(fun, xstar = pars[['xstar']], k = pars[['k']]))
# Method 2b: (outputs a data frame)
out <- data.frame(pars, value = mapply(fun, xstar = pars[['xstar']], k
= pars[['k']]))
HTH,
Dennis
On Wed, Jun 29, 2011 at 4:52 PM, nany23 <anna.botto at gmail.com> wrote:
Hello! I know that probably my question is rather simple ?but I' m a very beginner R-user. I have to numerically integrate ?the product of two function A(x) and B(x). The integretion limits are [X*; +inf] Function A(x) is a pdf function while B(x)=e*x is a linear function whose value is equal to 0 when the ?x < X* Moreover I have to iterate this process for different value of X* and for different pdf of the same type. I know the comand INTEGRATE but ?I can' t make it work. Which is the best function to do this and how does it work? Thank you very much in advanced!! -- View this message in context: http://r.789695.n4.nabble.com/Numerical-integration-tp3634365p3634365.html Sent from the R help mailing list archive at Nabble.com.
______________________________________________ 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 and provide commented, minimal, self-contained, reproducible code.
thank you very much for your suggestion!
I tried to do that with the psf I need to use: the 3 parameters Lognormal. I
did that with a single xstar and a single triplet of parameters to check it
works.[I put some numbers to make it woks , but actually they comes from
statistical analysis]
/# these are the 3 parameters
a<- 414.566
b<- 345.5445
g<- -0.9695679
xstar<- 1397.923
#I create a vector
pars <-expand.grid(xstar = xstar, a= a, b= b , g= g)
fun <- function(xstar, a,b,g,k) {
f <- function(x, xstar, a, b,g,k) f.lognorm(x) * k * x * (x >= xstar)
integrate(f, -Inf, Inf, xstar = xstar, a = a, b =b, g=g, k=k)$value
}
# Method 1: (outputs a data frame)
library(plyr)
out <- mdply(pars, fun)/
at this stage a warning message comes out: *Errore in k > -1e-07 : 'k' is
missing*
Any ideas of why this error come out ? I really don' know.....
Moreover , can I use the algorithm suggested iteratively by using a grid
of xstar values and different triplets of parameters ( instead of different
value of k) ?
Thank you so much for any help!!
--
View this message in context: http://r.789695.n4.nabble.com/Numerical-integration-tp3634365p3637092.html
Sent from the R help mailing list archive at Nabble.com.