Diyala

From InterSciWiki

back to Sample data for q-exponential

tip on Mac conversion from an excel column to a distribution for R

library(gsl) source("http://intersci.ss.uci.edu/wiki/Cosma/pareto.R") #use for plots source("http://intersci.ss.uci.edu/wiki/Cosma/disc.q.exp.estimation.R") #discrete estimation source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/tsal-v0.2.2.R") source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/disc.q.exp.estimation.R") source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/lnq-qexp-estimate.R") diyala=c(315000,250000,210000,210000,202500,135000,132600,122500,112500,108000,106000,100000,100000,90000,80000,72600,60000,52500,50000,50000,48000,45000,40000,40000,40000,40000,35000,30000,30000,29400,28900,25200,25000,24000,22500,21600,20000,19200,18000,18000,18000,16000,15000,14400,13600,13500,13000,12000,12000,10500,10200,10050,10000,9600,6600,6525,5000,3300,3220,2800,2800,2450,2250,1750,1700,1500,1500,1120,1050,950,900,700,600,560,420,400,400,360,350,320,300,300,270,270,240,180,180,135,120,100,90,60,40) #these data are square meters, not sizes or frequencies

###put your own data here and replicate: DISCRETE Disc.q.exp.estimation.R
rts=diyala
tsal.fit.disc(rts)$q
tsal.fit.disc(rts)$kappa

[edit] DISCRETE DISTRIBUTION MLE FIT (no xmin yet) Disc.q.exp.estimation.R

tsq=strtrim(as.character(tsal.fit.disc(rts)$q),5) #label for q
tsk=strtrim(as.character(tsal.fit.disc(rts)$kappa),5) #label for kappa
plot.survival.loglog(rts,ylab="Cumul. Probability, tsal.fit.disc",main="discrete MLE Survival function-Diyala data", sub=paste("zags=empiric,dash=Pareto,solid=Tsallis q=", tsq, "kappa=",tsk)) 
rts.tsal.fit.disc <- tsal.fit.disc (rts)#,xmin=1) 
curve(ptsal(x,rts.tsal.fit.disc$shape,rts.tsal.fit.disc$scale,lower.tail=FALSE),add=TRUE,col="blue") 
tsal.fit.disc(rts)$q     #
tsal.fit.disc(rts)$kappa #
#ADD PARETO TAIL FIT
rts.min <- 125000 #(eyeball estimate SET THIS CAREFULLY FOR POWER LAW)
rts.tailprob <- sum(rts>=rts.min)/length(rts) 
rts.pareto <- pareto.fit(rts,rts.min)
curve(rts.tailprob*ppareto(x,rts.min,rts.pareto$exponent,lower.tail=FALSE),from=rts.min,col="red",lty="dashed",add=TRUE)

[edit] CONTINUOUS DISTRIBUTION MLE FIT Tsal-v0.2.2.R

tsq=strtrim(as.character(tsal.fit(rts)$q),5)
tsk=strtrim(as.character(tsal.fit(rts)$kappa),5)
plot.survival.loglog(rts,ylab="Cumulative probability",main="continuous MLE Survival function-Diyala data", sub=paste("zags=empiric,dash=Pareto,solid=Tsallis q=", tsq, "kappa=",tsk)) 
rts.tsal <- tsal.fit(rts,xmin=1) 
curve(ptsal(x,rts.tsal$shape,rts.tsal$scale,lower.tail=FALSE),add=TRUE,col="blue")
tsal.fit(rts)$q     #
tsal.fit(rts)$kappa # 
#ADD PARETO TAIL FIT 
rts.min <- 125000 #(eyeball estimate SET THIS CAREFULLY FOR POWER LAW) 
rts.tailprob <- sum(rts>=rts.min)/length(rts) 
rts.pareto <- pareto.fit(rts,rts.min) 
curve(rts.tailprob*ppareto(x,rts.min,rts.pareto$exponent,lower.tail=FALSE),from=rts.min,col="red",lty="dashed",add=TRUE)