べき分布の指数の推定関数
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)