Blackouts.q

From InterSciWiki

Jump to: navigation, search

blackouts <- c(7500000,2100000,2085000,1660000,1600000,1500000,1500000,1500000,1300000,1140000,899000,877000,875000,870000,725000,660000,650000,600000,600000,598000,570000,557354,500000,500000,500000,490000,464000,460000,404000,400000,400000,385000,375000,363476,360000,350000,350000,320831,315000,312000,300000,300000,300000,300000,290000,284000,272000,272000,258000,257718,250000,246000,246000,242910,240000,238000,235000,234000,230000,219000,210882,207200,206000,203000,200000,191000,190000,173000,173000,166000,164500,163000,160000,160000,160000,158000,148000,147000,146000,145000,145000,142000,133000,130000,130000,128000,126000,124000,122000,120000,120000,115000,114500,114000,113200,112000,106850,100000,100000,100000,100000,95630,95000,94285,92000,92000,91000,90000,88000,82500,81000,80000,80000,75000,75000,74000,71000,71000,71000,70000,70000,70000,70000,66005,65000,65000,63500,62000,60000,60000,60000,60000,60000,59000,58000,56000,56000,55000,55000,53000,51000,51000,50462,50000,50000,50000,

   + 50000,50000,50000,50000,50000,50000,48000,46000,45000,43696,43000,40911,40000,40000,40000,39500,38500,37000,36073,35000,33000,32000,32000,30500,30001,30000,29900,29000,29000,26334,25000,25000,25000,25000,25000,25000,24506,24000,20000,20000,19000,18819,18351,18000,18000,17000,15000,15000,14273,12000,11529,11000,10300,10000,10000,9000,8000,7500,5300,4150,2900,2000,1800,1646,1000) 
   blackouts.min <- 230000 
   blackouts.tailprob <- sum(blackouts>=blackouts.min)/length(blackouts) 
   blackouts.pareto <- pareto.fit(blackouts,blackouts.min) 
   blackouts.random <- blackouts.min*((1-ppareto(blackouts,blackouts.min,blackouts.pareto$exponent,lower.tail=FALSE))^(1/(1-blackouts.pareto$exponent)))
   blackouts.tsal <- tsal.fit(blackouts,xmin=1000)
   #blackouts.tsal <- tsal.fit(ifelse(blackouts<=blackouts.min,blackouts,blackouts.random),xmin=1000)
   plot.survival.loglog(blackouts,ylab="Cumulative probability",main="Survival function of blackouts data", sub="zagged=empirical, solid=Tsallis, dashed=Pareto, ") 
   curve(ptsal(x,blackouts.tsal$shape,blackouts.tsal$scale,lower.tail=FALSE),add=TRUE,col="blue") 
   curve(blackouts.tailprob*ppareto(x,blackouts.min,blackouts.pareto$exponent,lower.tail=FALSE),from=blackouts.min,col="red",lty="dashed",add=TRUE) 
   blackouts.tsal
   
   $type
   [1] "tsal"
   
   $q
   [1] 1.565541
   
   $kappa
   [1] 113940.9
   
   $shape
   [1] 1.768219
   
   $scale
   [1] 201472.4
   
   $loglike
   [1] -2788.139
   
   $n
   [1] 211
   
   $xmin
   [1] 1000
   
   $method
   [1] "mle.equation"
Personal tools