不完全ベータ関数(1-alpha,0,x)の漸近展開の実装の例 後半は部分分数分解やテーラー展開で積分
bx_real<-function(x,alpha){
z<-x
if(length(x[abs(x)>1])>=1){
n<--log(0.00001,base=min(x[abs(x)>1]))
str1<-0
for(i in 0:n){
str1<-str1+x[abs(x)>1]^-i/(i+alpha)
}
stra<-alpha*gamma(2-alpha)*gamma(-alpha)/(gamma(1-alpha)^2)*str1
strb<-cos(pi*alpha)*gamma(2-alpha)*gamma(alpha)*x[abs(x)>1]^alpha
z[abs(x)>1]<-1/(alpha-1)*x[abs(x)>1]^-alpha*(stra+strb)
}
if(length(x[abs(x)<1])>=1){
n<-log(0.00001,base=max(x[abs(x)<1]))
str1<-0
for(i in 1:n){
str1<-str1+(x[abs(x)<1]^i)/(i-alpha)
}
z[abs(x)<1]<-x[abs(x)<1]^(-alpha)*str1
}
if(length(x[abs(x)==1])>=1){
z[abs(x)==1]<-Inf
}
-z
}
bx_int<-function(x,alpha0,N){
alpha<-abs(alpha0)
if(alpha0<0){
str1<-log(x)
for(k in 1:alpha){
str1<-str1+ choose(alpha,k)*(N^(-k))/k*x^k
}
}else{
if(alpha==0){
str1<-log(x)
}else{
str1<-log(x)-log(abs(x+N))
if(alpha>=2){
for(k in 1:(alpha-1)){
str1<-str1-N^k/(-k)*1/(x+N)^k
}
}
}
}
str1
}
i_bx0<-function(z,alpha){
ans1<-NA
f12<-function(x){bx(x,alpha)-z}
try(x1<-uniroot(f12,c(0,1)))
try(ans1<-x1$root)
ans2<-NA
try(x2<-uniroot(f12,c(1.001,10^6)))
try(ans2<-x2$root)
c(ans1,ans2)
}
th_y<-function(tt,r,alpha,x0,Y){
t0<-bx(1+x0/Y,alpha)
Y*(-1+i_bx0(t0-r*(tt),alpha))
}
th_t_real<-function(xx,r,alpha,x0,Y){
t0<-bx_real(1+x0/Y,alpha)
td<-bx_real(1+xx/Y,alpha)
tt<-(td-t0)/r
tt
}
th_t_int<-function(xx,r,alpha,x0,Y){
t0<-bx_int(x0,alpha,Y)
td<-bx_int(xx,alpha,Y)
tt<-(td-t0)/r
tt
}
th_t<-fuction(xx,r,alpha,x0,Y){
if(!is.integer(alpha)){
th_t_real(xx,r,alpha,x0,Y)
}else{
th_t_int(xx,r,alpha,x0,Y)
}
}
