# 3. přednáška # rm(list=ls()) # podejme data data(EU2010) head(EU2010) # které země patří do V4? (V4 = EU2010[,"země"] %in% c("CZ","SK","PL","HU")) # omezme se na země V4 (V4_2010 = EU2010[V4,]) # zpřístupněme informace o zemích V4 attach(V4_2010) # # POZOR, lesy JSOU V TISÍCÍCH ha!!! # přepočítejme tedy proměnnou lesy na km^2 (lesy = 10*lesy) # names(lesy) = země lesy # # průměr NA JEDNU ZEMI mean(lesy) # Giniho koeficient # rozdíly plochy lesů pro všechy dvojice outer(lesy,lesy,"-") abs(outer(lesy,lesy,"-")) # střední diference (Delta = mean(abs(outer(lesy,lesy,"-")))) # Giniho koeficeint (G = Delta/(2*mean(lesy))) # měří nerovnoměrnost plochy lesů mezi čtyřmi zeměmi library(reldist) ## ## pokud nebude balíček reldist, pro případ stejných ## vah si pomůžeme pomocí ## gini = function(x) mean(abs(outer(x,x,"-")))/(2*mean(x)) ## gini(lesy) # jaká část lesů je v jednotlivých zemích: (podílLesy = lesy/sum(lesy)) sort(podílLesy) ## kumulativní relativní plocha lesů: cumsum(sort(podílLesy)) # plot((0:4)/4,c(0,cumsum(sort(podílLesy))),pch=3,asp=1, xlab="státy",ylab="lesy") lines(c(0,1),c(0,1)) lines(c(0,1,1,0,0),c(0,0,1,1,0),lty=3) lines((0:4)/4,c(0,cumsum(sort(podílLesy)))) # # Co se tu dá vyčíst? # # Jakou část lesů má na svém území "chudší" polovina zemí, # tj. Slovensko a Maďarsko? cumsum(sort(podílLesy)) # SK HU CZ PL # 0.1214819 0.2492948 0.4158465 1.0000000 # odpověď: čvrtinu lines(c(0.5,0.5,0),c(0,cumsum(sort(podílLesy))[2],cumsum(sort(podílLesy))[2]), lty=2,col="red") # polygon(c(x=(0:4)/4,0),y=c(0,cumsum(sort(podílLesy)),0),col="yellow") # Velikost žluté plochy vztažená k velikosti "dolného trojúhelníku" # je rovna Giniho koeficientu koncentrace gini(lesy) # až dosud porovnáváme čtyři státy, jejich lesní bohazství # nově se budeme zajímat o zalesnění celého území V4, # nakolik nerovnoměrně lesy pokrývají území V4 # # jaká část rozlohy jednotlivých států je zalesněna (zalesnění = lesy/rozloha) # zajímá nás variační řada sort(zalesnění) # v jakém pořadí vzít zalesnění, # abych dostal variační řadu? (ord = order(zalesnění)) # kontrola zalesnění[ord] # doplňme název také k rozloha names(rozloha) = names(lesy) # jak přibývají (RELATIVNĚ) kilometry čtvereční rozlohy? (rozloha[ord]/sum(rozloha)) # kumulativní podíly jsou tedy cumsum((rozloha[ord]/sum(rozloha))) # jak přibývají [RELATIVNĚ] kilometry čtvereční lesů? (lesy[ord]/sum(lesy)) # kumulativní podíly jsou tedy cumsum((lesy[ord]/sum(lesy))) # plot(c(0,cumsum(rozloha[ord]))/sum(rozloha), c(0,cumsum(lesy[ord]))/sum(lesy),pch=3,asp=1, xlab="území",ylab="lesy") lines(c(0,cumsum(rozloha[ord]))/sum(rozloha), c(0,cumsum(lesy[ord]))/sum(lesy)) lines(c(0,1),c(0,1)) lines(c(0,1,1,0,0),c(0,0,1,1,0),lty=3) # jak rychle stoupají jednotlivé úsečky? # osa x (rozloha/sum(rozloha))[ord] # osa y (lesy/sum(lesy))[ord] # polygon(x=c(0,cumsum(rozloha[ord]),0)/sum(rozloha), y=c(0,cumsum(lesy[ord]),0)/sum(lesy),col="yellow") # obdoba výpočtu na slidu 66 gini(zalesnění,rozloha) ## ## pokud není balíček reldist, pak gini <- function(yPrum,x=rep(1,length(yPrum))){ permutace <- order(yPrum); k <- length(yPrum) yPrum <- yPrum[permutace]; x <- x[permutace] Xrel <- cumsum(x)/sum(x); Yrel <- cumsum(x*yPrum)/sum(x*yPrum) gini <- sum(Xrel[-k]*Yrel[-1]-Xrel[-1]*Yrel[-k]) return(gini=gini) } gini(zalesnění,rozloha) ## # průměrné zalesnění velého území V4 (yMean = sum(lesy)/sum(rozloha)) names(yMean) = "V4" yMean weighted.mean(zalesnění,w=rozloha) # (X = cumsum(rozloha[ord])) (Y = cumsum(lesy[ord])) (k = length(X)) yMean = Y[k]/X[k] # výpočet průměru jinak (Delta = 2*sum(X[-k]*Y[-1]-X[-1]*Y[-k])/X[k]^2) (G = Delta/(2*yMean)) # k rovnoměrnosti přispívá skutečnost, # že uvnitř každé země uvažujeme rovnoměrné # zalesnění # prijem=c(200,150,80,70,60,60,20,20,18,18,15,15,10,10) # gini(prijem) (Skup = factor(c(rep("A",2),rep("B",4),rep("C",8)))) data.frame(Skup,prijem) (ni = table(Skup)) (prumery = tapply(prijem,Skup,mean)) (GB = gini(prumery,ni)) # nerovnoměrnost průměrů (Gi = tapply(prijem,Skup,gini)) # # jiná míra nerovnoměrnosti - Theilův index # library(ineq) Theil(prijem) mean((prijem/mean(prijem))*log(prijem/mean(prijem))) weighted.mean(log(prijem/mean(prijem)),w=prijem/mean(prijem)) ## ## nebude-li knihovna, pak ## Theil = function(x) mean((x/mean(x))*log(x/mean(x))) # každý násobek vah je stejně dobrý: weighted.mean(log(prijem/mean(prijem)),w=prijem) # # (prumer = mean(prijem)) (n = sum(ni)) (index=Theil(prijem)) (indexy=tapply(prijem,Skup,Theil)) (TB = sum(ni*prumery/prumer*log(prumery/prumer))/n) (TW = sum(ni*prumery/(n*prumer)*indexy)) # kontrola c(index,TB+TW) c(TB,TW)/index*100 #