# W&J Example page 67. por Eduardo S. Cypriano

library(zipfR)

lum<-seq(0.0,5,0.0001)
a <- -0.5


#Luminosity function
lf<-function(x,alpha) return (x^alpha * exp(-x) / gamma(alpha+1))
clf<-function(x,alpha)  return ( Igamma(alpha+1,x)/ gamma(alpha+1)) 



dev.set(2)

plot(lum,lf(lum,a), type="l",
                    ylim=c(0,1),
                    xlab="Luminosidade (L*)",
                    ylab="Densidade de Probabilidade")

lines(lum,clf(lum,a))

# n'th Order Statistics

g_i<-function(y,alpha,i,N) return (factorial(N) / (factorial(i-1)*factorial(N-i)) * 
                             clf(y,alpha)^(i-1) * (1-clf(y,alpha))^(N-i) * lf(y,alpha)   
                            )


if (length(dev.list())<2) {
    X11()
 }

dev.set(3)
plot(lum,lf(lum,a), type="l",
                    ylim=c(0,1),
                    xlab="Luminosidade (L*)",
                    ylab="Densidade de Probabilidade")

lines(lum,g_i(lum,a,10,10),lty=("dashed"),col="red")
lines(lum,g_i(lum,a,100,100),lty=("longdash"),col="blue")


#lines(lum,g_i(lum,-0.5,9,10),lty=("dotted"),col="green")