# R/gammaIntegral.R In rGammaGamma: Gamma convolutions for methylation array background correction.

#### Documented in gammaIntegral

```gammaIntegral <- function(total, params, offset=50, minx=1) {

# require('gsl') # for hyperg_1F1 in the conditional expectation of the signal

## this bit is the most obvious "farm me out to C++" piece of all...
if(length(total) > 1) return(sapply(total, gammaIntegral, params=params))

g = params
a = params
d = params
b = params
bg.mean = d * b
bg.sd = sqrt( d * b * b )
## cat('total =',total,'... bg.mean =',bg.mean,'... bg.sd =',bg.sd,"\n")
if(total > ( bg.mean + ( 3 * bg.sd ) )) {
return(pmax(total - bg.mean, minx))
} else {
## FIXME: need to write this as a function object for C++ to integrate it
res = try(
integrate(
function(x) {
(exp(x*((1/b)-(1/a)))*(total**(1-g-d))*((total-x)**(d-1))*(x**(g-1)))*
(1/(beta(g,d)*hyperg_1F1(g, g+d, total*((1/b)-(1/a)), strict=F)))*x
},
0, total
)\$value # else will return a list with value, abs.error, subdivisions, ...
) # i.e., integrate(PrSignalGivenTotal, /* from */ 0, /* to */ total);
if(class(res) == 'try-error') {
return(pmax(total-bg.mean, minx)+offset)
} else {
return(pmax(res, minx)+offset)
}
}

}
```

## Try the rGammaGamma package in your browser

Any scripts or data that you put into this service are public.

rGammaGamma documentation built on May 2, 2019, 6:38 a.m.