Rでべき分布(両対数)のCDFをきれいに書く。
a1<-0.5/runif(30000)
a2<-0.3/runif(30000)
a3<-0.1/runif(30000)
a4<-0.01/runif(30000)
plt1<-par()$plt
plt2<-plt1*c(1,1,1,1)
plt2[1]<-0.2
par(plt=plt2)plot(sort(a1),length(a1):1/length(a1),log="xy",type="l",col=1,tck=0.01,cex.axis=2,cex.lab=1.5,xaxt="n",yaxt="n",xlab="",ylab="",xlim=c(10^-3,10^3),lty=1,lwd=3)
points(sort(a2),length(a2):1/length(a2),type="l",col=1,lty=2,lwd=3)
points(sort(a3),length(a3):1/length(a3),type="l",col=1,lty=3,lwd=3)
points(sort(a4),length(a4):1/length(a4),type="l",col=1,lty=4,lwd=3)
j<-as.integer(exp(seq(log(length(a1)),log(1),length.out=20)))
points(sort(a1,decreasing=T)[j],(1:length(a1)/length(a1))[j],type="p",col=1,lty=4,lwd=3,pch=1)
j<-as.integer(exp(seq(log(length(a1)),log(1),length.out=20)))
points(sort(a2,decreasing=T)[j],(1:length(a2)/length(a2))[j],type="p",col=1,lty=2,lwd=3,pch=2)j<-as.integer(exp(seq(log(length(a1)),log(1),length.out=20)))
points(sort(a3,decreasing=T)[j],(1:length(a3)/length(a3))[j],type="p",col=1,lty=2,lwd=3,pch=3)
j<-as.integer(exp(seq(log(length(a1)),log(1),length.out=20)))points(sort(a4,decreasing=T)[j],(1:length(a4)/length(a4))[j],type="p",col=1,lty=2,lwd=3,pch=4)
curve(0.004/x,add=T,from=10^-0.3,to=10^1.3,lwd=1,lty=1)
axis(1,at=c(10^-3,10^(0),10^3),label=c(expression(10^-3),expression(10^0),expression(10^3)),cex.axis=2,tck=0.03)
axis(1,at=10^(min(log(axTicks(1),10)):max(log(axTicks(1),10))),label=F,tck=0.03)
axis(1,at=c(1:9%x%10^(min(log(axTicks(1),10)):max(log(axTicks(1),10)))),label=F,tck=0.01)axis(3,at=10^(min(log(axTicks(1),10)):max(log(axTicks(1),10))),label=F,tck=0.03)
axis(3,at=c(1:9%x%10^(min(log(axTicks(1),10)):max(log(axTicks(1),10)))),label=F,tck=0.01)
axis(2,at=c(10^0,10^(-2),10^(-4)),label=c(expression(10^0),expression(10^-2),expression(10^-4)),cex.axis=2,tck=0.03,las=2)
axis(2,at=c(10^(0:-5)),label=F,tck=0.03)
axis(2,at=c(1:9%x%10^(0:-5)),label=F,tck=0.01)
axis(4,at=c(10^(0:-5)),label=F,tck=0.03)
axis(4,at=c(1:9%x%10^(0:-5)),label=F,tck=0.01)
mtext(1,line=2.4,font=7,cex=1.5,text="x")
mtext(2,line=5.2,font=7,cex=1.5,text="CDF")
legend(x=10^-3,10^-2.9,lty=c(1,2,3,4,1),lwd=c(2,2,2,2,1),c("N=1","N=100","N=10000","N=900000",expression(""%prop%1/x)),pch=c(1,2,3,4,-1))
par(plt=plt1)
