CODEWITHSUNDEEP
rintegralnumerical-integrationcalculus
DPdl
DPdl

Reputation: 755

Integrating in R

I am trying to compute L(psi) as given here in R. I have the following values. 

nb <- 100
tb <- 25
ns <- 90
ts <- 15
A0 <- 1
S_norm <- 0.4

R <-tb/ts
y_meas <- (ns-nb/R)/A0
sigma_meas = sqrt(ns+(nb+1)/R^2)/A0

I am very confused on how I can integrate L(psi), say from -10 to 10. Because I am integrating with respect to log A.

Upvotes: 1

Views: 628

Answers (2)

Severin Pappadeux
Severin Pappadeux

Reputation: 20080

On top of an excellent answer by @SandipanDey, if you could extend limits to -Infinity...+Infinity, there is a better way to integrate functions with e-x2 kernel: Gauss-Hermite quadrature, and there is an R package for that.

Simple example:

library(gaussquad)

n.quad <- 128 # integration order

# get the particular (weights,abscissas) as data frame
# with 2 observables and n.quad observations
rule <- ghermite.h.quadrature.rules(n.quad, mu = 0.0)[[n.quad]]

# test function - integrate 1 over exp(-x^2) from -Inf to Inf
# should get sqrt(pi) as an answer
f <- function(x) {
    1.0
}

q <- ghermite.h.quadrature(f, rule)
print(q - sqrt(pi))

Upvotes: 1

Sandipan Dey
Sandipan Dey

Reputation: 23101

You can substitute for logA and for a fixed value of psi you can integrate as follows:

enter image description here

psi <- 5
integrate(function(x) exp(-0.5*(((x/A0)/S_norm)^2 + ((psi-y_meas*A0/exp(x))/sigma_meas)^2)),
                                                                                   -10, 10)
# 0.1775989 with absolute error < 6.6e-05

Upvotes: 2

Related Questions