べき分布の指数の推定関数
powerFit<-function(x){
D<-0
P<-0
alpha123<-0
N123<-100
xminmin<-0
xminb<-0
kazu<-0;
Dp<-function(x){
i<-1:100
2*sum(((-1)^(i-1))*exp(-2*i*i*x*x))
}
h1b<-x
ALPHA<-0
XMIN<-0
PVALE<-0
i<-1
xminxmin<-exp(seq(min(log(h1b)),max(log(h1b)),length.out=N123) )
h1c<-x
while(length(h1c)>=2){
xmin<-xminxmin[i]
h1c<-h1b[h1b>=xmin]
alpha<-1+length(h1c)*1.0/sum(log(h1c/xmin))-1
xminb[i]<-xmin
D[i]=max(abs(xmin^alpha/(sort(h1c))^(alpha)-length(h1c):1/length(h1c)))
alpha123[i]<-alpha
cat(D[i],alpha,length(h1c),"\n")
P[i]<-Dp(sqrt(length(h1c))*D[i]);
cat(D[i],alpha123[i],length(h1c),length(alpha123),length(D),length(P),"\n")
i<-i+1;
}
Db<-D[!is.na(D) & D != Inf & !is.na(alpha123) & alpha123 != Inf & !is.na(P) & P != Inf & !is.na(xminb) & xminb != Inf]
alpha123b<-alpha123[!is.na(D) & D != Inf & !is.na(alpha123) & alpha123 != Inf & !is.na(P) & P != Inf & !is.na(xminb) & xminb != Inf]
Pb<-P[!is.na(D) & D != Inf & !is.na(alpha123) & alpha123 != Inf & !is.na(P) & P != Inf & !is.na(xminb) & xminb != Inf]
xminbb<-xminb[!is.na(D) & D != Inf & !is.na(alpha123) & alpha123 != Inf & !is.na(P) & P != Inf & !is.na(xminb) & xminb != Inf]
xmin<-xminbb[which(Db==min(Db))]
alpha<-alpha123b[which(Db==min(Db))]
pb<-Pb[which(Db==min(Db))]
c(alpha,xmin,pb)
#data.frame(Db,alpha123b,xminbb)
}
95%信頼区間
powerFit2<-function(x){
D<-0
P<-0
alpha123<-0
N123<-100
xminmin<-0
xminb<-0
kazu<-0;
Dp<-function(x){
i<-1:100
2*sum(((-1)^(i-1))*exp(-2*i*i*x*x))
}
h1b<-x
ALPHA<-0
XMIN<-0
PVALE<-0
i<-1
xminxmin<-exp(seq(min(log(h1b)),max(log(h1b)),length.out=N123) )
h1c<-x
while(length(h1c)>=2){
xmin<-xminxmin[i]
h1c<-h1b[h1b>=xmin]
alpha<-1+length(h1c)*1.0/sum(log(h1c/xmin))-1
xminb[i]<-xmin
D[i]=max(abs(xmin^alpha/(sort(h1c))^(alpha)-length(h1c):1/length(h1c)))
alpha123[i]<-alpha
cat(D[i],alpha,length(h1c),"\n")
P[i]<-Dp(sqrt(length(h1c))*D[i]);
cat(D[i],alpha123[i],length(h1c),length(alpha123),length(D),length(P),"\n")
i<-i+1;
}
Db<-D[!is.na(D) & D != Inf & !is.na(alpha123) & alpha123 != Inf & !is.na(P) & P != Inf & !is.na(xminb) & xminb != Inf]
alpha123b<-alpha123[!is.na(D) & D != Inf & !is.na(alpha123) & alpha123 != Inf & !is.na(P) & P != Inf & !is.na(xminb) & xminb != Inf]
Pb<-P[!is.na(D) & D != Inf & !is.na(alpha123) & alpha123 != Inf & !is.na(P) & P != Inf & !is.na(xminb) & xminb != Inf]
xminbb<-xminb[!is.na(D) & D != Inf & !is.na(alpha123) & alpha123 != Inf & !is.na(P) & P != Inf & !is.na(xminb) & xminb != Inf]
xmin<-xminbb[which(Db==min(Db))]
alpha<-alpha123b[which(Db==min(Db))]
pb<-Pb[which(Db==min(Db))]
#x<-runif(10000)^(-1/1.5)
#q<-powerFit(x)
x1<-x[x>=xmin]
#alpha_real=1+length(x1)*(sum(log(x1/min(x1))))^-1-1
xu<-qnorm(0.025,mean=alpha+1,sd=(alpha+1-1)/sqrt(length(x1)))-1
xb<-qnorm(1-0.025,mean=alpha+1,sd=(alpha+1-1)/sqrt(length(x1)))-1
c(alpha,xmin,pb,xu,xb)
#data.frame(Db,alpha123b,xminbb,xu,xb)
}
x<-runif(10000)^(-1/1.3)
q<-powerFit2(x)
