R/ainc_gamma.R

In fxjulio/slash: Slash reparametrized

ainc.gamma <- function( x, a){

erfc <- function(x) 2 * pnorm(x * sqrt(2), lower = FALSE)

mkgammak <- c(
1/12,
1/288,
-139/51840,
-571/2488320,
163879/209018880,
5246819/75246796800,
-534703531/902961561600,
-4483131259/86684309913600,
432261921612371/514904800886784000,
6232523202521089/86504006548979712000,
-25834629665134204969/13494625021640835072000,
-1579029138854919086429/9716130015581401251840000,
746590869962651602203151/116593560186976815022080000,
1511513601028097903631961/2798245444487443560529920000,
-8849272268392873147705987190261/299692087104605205332754432000000,
-142801712490607530608130701097701/57540880724084199423888850944000000,
2355444393109967510921431436000087153/13119320805091197468646658015232000000,
2346608607351903737647919577082115121863/155857531164483425927522297220956160000000,
-2603072187220373277150999431416562396331667/1870290373973801111130267566651473920000000,
-73239727426811935976967471475430268695630993/628417565655197173339769902394895237120000000,
34856851734234401648335623107688675640839679447003/2601648721812516297626647395914866281676800000000,
909773124599542506852275229422593983242880452145053/811714401205505084859513987525438279883161600000000,
-1527335577854677023023224272800947125313629267269390501/9740572814466061018314167850305259358597939200000000,
-183856455668177802003316143799518064719008299958634826921/14026424852831127866372401704439573476381032448000000000,
2583312098861137963745902036370496943872138148651712093816393/1178219687637814740775281743172924172016006725632000000000,
5180134290822682443757710427952467581918233549140896702364013/28277272503307553778606761836150180128384161415168000000000,
-527550309097873396592733540579928993424142983691519876840948418433873/14613128884259277641708402381685690086746366936130519040000000000,
-2114866241537081164613223324215572812504648703648482437460602956015127/701430186444445326802003314320913124163825612934264913920000000000,
180394412915538782140015777241228025103785450235726235175126981743099027459/260932029357333661570345232927379682188943128011546547978240000000000,
3226140192053936286912811949056082647586604417173687729452086326364208020303641/55891640688340870308367948893044727924871618020073270576939008000000000000
)


Pk <- list(c(1))
for( k in 2:length(mkgammak) ){
	temp <- c()
	temp[1] <- (2*k-1)*Pk[[k-1]][1] 
	if( k > 2 ){
		j <- 1:(2*k-4)
		temp[ j + 1 ] <- (2*k - 1 - j)*(Pk[[ k-1 ]][j+1] + Pk[[ k -1 ]][j] )
	}
	temp[2*k-2] <- 2*Pk[[ k-1 ]][2*k-3]
	temp[2*k-1] <- mkgammak[ k-1 ]

	Pk <- append(Pk, list(temp))
}

sapply( x, function(x){
	lambda <- x/a
	mu <- lambda - 1
	eta <- sign(mu)*sqrt(2*(mu-log(1+mu)))

	ck <- sapply( 1:length(mkgammak), function(k) {
	Qk <- (1+mu)*sum( Pk[[k]]*mu^(0:(2*k-2)) ) + mkgammak[k]*mu^(2*k)
	Ak <- 2^k*gamma(k+1/2)/gamma(1/2)
	(-1)^k*( Qk/mu^(2*k+1) - Ak/eta^(2*k+1) )  # ck
	})

	ck <- c( 1/mu - 1/eta  , ck )
	ck1 <- -eta*ck ## falta + dc[k-1]/deta
	
	k <- 0:(length(ck)-1)
	
	# dP/dx
	c(
	1/a/2/eta*(2-2/(1+mu))*( pi^(-1/2)*exp(-eta^2*a/2)*sqrt(a/2) - 
		sqrt(a/(2*pi))*( 1 -1/(mu+1)*sum( (-1)^k*c(1,mkgammak)*a^(-k) )*(mu+1)*eta/mu)*exp(-(1/2)*a*eta^2)  ),	
	1/2*erfc(-eta*sqrt(a/2)) - (2*pi*a)^(-1/2)*exp(-(1/2)*a*eta^2)*sum( ck*a^(-k)  )
	)
})

}