Sample data for q-exponential
From InterSciWiki
Back to Using the new discrete estimator and producing sampling distribution plots --
Uruk longitudinal scaling
library(gsl)
rts=c(98,87,65,40,33,20,18,15,14,13,12,11,9,8)
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")
source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/yule.R")
source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/lnorm.R")
source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/disclnorm.R")
source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/exp.R")
source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/discexp.R")
source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/powerexp.R")
source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/discpowerexp.R")
source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/powerexp-exponential-integral.R")
source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/weibull.R")
source("http://intersci.ss.uci.edu/wiki/pub/Shalizi/discweib.R")
tsal.fit.disc(rts)$q
tsal.fit.disc(rts)$kappa
###put your own data here and replicate
rts=LateUrukSiteSizeTimes10
LateUrukSiteSizeTimes10=c(2500,500,255,250,240,140,135,115,110,100,100,100,90,79,78,78,77,67,66,65,60,60,
60,60,58,58,55,55,53,52,50,50,49,48,48,48,47,44,42,42,40,40,40,40,40,38,36,35,
34,30,30,30,30,30,29,29,28,26,26,26,26,24,20,20,20,20,20,20,19,19,18,18,17,17,
17,17,17,16,16,15,15,14,12,11,11,10,10,10,10,10,10,10,10,10,9,9,9,9,9,9,8,8,8,
8,7,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,4,4,4,3,3,3,2,2,2,2,2,2,1,1,1,1)
rts=LateUrukSiteSizeTimes10
references at Using the new discrete estimator and producing sampling distribution plots
Contents |
[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-Warka Late Uruk data", sub=paste("zags=empirical,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 #[1] 1.463075
tsal.fit.disc(rts)$kappa #[1] 24.61816
#ADD PARETO TAIL FIT
rts.min <- 40 #(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-Warka Late Uruk data", sub=paste("zags=empirical,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 #[1] 1.485
tsal.fit(rts)$kappa #[1] 23.37
#ADD PARETO TAIL FIT
rts.min <- 40 #(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)
#ADD LOGNORMAL fit & OTHERS: lnorm.R, exp.R, discpowerexp.R, yule.R, Weibull.R -- TO SEE OTHER heavy tail distributions
#Thanks to Bela Nagy for solution to comparing to Log-normal in three steps.
#Step 1. Fix a bug in the function lnorm.fit.moments in the file
#http://intersci.ss.uci.edu/wiki/pub/Shalizi/lnorm.R
lnorm.fit.moments <- function(x, threshold = 0) {
x <- x[x>=threshold]
x <- x-threshold
LogData <- log(x)
M = mean(LogData)
#SD = sd(LogData) # replace this line with
SD = stats::sd(LogData)
Lambda = lnorm.loglike.shift(x,M,SD,threshold)
fit <- list(type="lnorm", meanlog=M, sdlog=SD, loglike=Lambda)
return(fit)
}
# Step 2. Use the fixed function by calling lnorm.fit
rts.lnorm <- lnorm.fit(rts)
# Step 3. Add the new curve to the existing plot
curve(plnorm(x,meanlog=rts.lnorm$meanlog,sdlog=rts.lnorm$sdlog,lower.tail=FALSE),from=rts.min,col="green",lty="dashed",add=TRUE)
[edit] Your Data
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)
###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 ###please leave your data and results on this page
###put your own data here and replicate: CONTINUOUS Tsal.R tsal.fit(rts)$q tsal.fit(rts)$kappa ###please leave your data and results on this page
###put your own data here and replicate: Lng(x) by x lnq-qexp-estimate.R - might be a better estimate for very small samples rts=c( ) tsal.qlog.estimate(rts)$q #####This is not the right command tsal.qlog.estimate(rts)$kappa #####This is not the right command ###please leave your data and results on this page
[edit] SKIP ALL BELOW
# Plot the density of the discrete MLE's values of q. # of the variation comes from sampling. plot(density(many.ml.estimates.q.small)) # To overlay the other densities, we need to store them as separate objects qlog.temp <- density(many.qlog.estimates.q.small) cmle.temp <- density(many.ml.estimates.q.small.continuous) # Now do the nice plotting - start wth the discrete MLE's sampling distribution plot(density(many.ml.estimates.q.small),xlim=c(0,10),main="Sampling distribution of q estimates: discrete data",xlab="black discrete MLE, green continuous MLE, red q-logarithm",sub="n=50")#change|to 10 from 10 # add lines showing the other two - note using $x to extract the component named "x", # and $y to extract the component named "y". Use help(lines) to see more about ways # of adding lines to plots. lines(qlog.temp$x,qlog.temp$y,col="red") lines(cmle.temp$x,cmle.temp$y,col="green")
# add a vertical line at the true q value. Use help(abline) to see more about ways of adding # simple lines to plots. abline(v=1.7,col="blue")
# Now re-do this with just the discrete MLE and the q-logarithm estimator. plot(density(many.ml.estimates.q.small),xlim=c(0.5,5),main="Sampling distribution of q estimates: discrete data",xlab="black discrete MLE, red q-logarithm",sub="n=50") lines(qlog.temp$x,qlog.temp$y,col="red") abline(v=1.7,col="blue")
##you can skip this part, its for examining the variation in the estimates, not all the elements are here (can be done in the previous page) #Since the true q is unchanging, all INSERT YOUR Q=1.46,KAPPA=24.62 ##num=length(rts) ##many.ml.estimates.q.small <- replicate(1e2,tsal.fit.disc(rtsal.disc(num,q=1.46,kappa=24.62,precompute=500))$q) ##many.qlog.estimates.q.small <- replicate(1e2,qlog.estimate(rtsal.disc(num,q=1.46,kappa=24.62,precompute=500))$q) ##many.ml.estimates.q.small.continuous <- replicate(1e2,tsal.fit(rtsal.disc(num,q=1.46,kappa=24.62,precompute=500))$q)
