2 Ajustament de poblacions

2.1 Llegir fitxers de dades

Llegir fitxers de dades en python és molt senzill. De fet, inclou un paquet natiu que s’anomena csv i que permet llegir fitxers en format csv (comma separated values) directament (o també ho podríem fer en un altre format, per descomptat). Donat que usarem les funcionalitats del scipy ens interessarà que ens carregui les dades com una array:

from numpy import *
from scipy.io import write_array
from scipy.io import read_array
data = zeros((3,3)) # creem una matriu qualsevol.
#Escrivim les dades:
write_array("fitxer.txt", data)
#Les llegim:
data = read_array(file("fitxer.txt"))

Observem que data ens inclou les dades que volíem. Anem a fer-ho de manera més sofisticada. Anem a triar dues poblacions d’una mateixa província espanyola i a veure la seva evolució. Per exemple, per veure dos comportaments diferents triarem, de la província de Madrid les localitats de “Pozuelo de Alarcón” i de “Olmedo de la Sierra”. Ens descarreguem les dades històriques de població que podem trobar a

http://www.ine.es/jaxiBD/tabla.do?per=12&type=db&divi=DPOH&idtab=30http://www.ine.es

Any	Olmeda de las Fuentes	Pozuelo de Alarcón
1900	492	1.873
1910	496	2.436
1920	455	2.677
1930	372	4.064
1940	355	2.517
1950	420	5.105
1960	332	9.412
1970	181	16.784
1981	152	31.228
1991	120	48.297

Copiem i enganxem les dades (o bé ens les descarreguem) en un fitxer anomenat “dades.csv” (separació amb espais o el que triem). L’obrim amb un editor de text i mirem que no hi quedin símbols estranys. Comentem la primera fila (que conté els noms de les columnes amb un símbol #). Anem a llegir-ho i a representar-ho

from numpy import *
from scipy.io import write_array
from scipy.io import read_array
data = read_array(file("dades.csv"), columns=(0,1,2))
# Indiquem quin es el separador i les línies que llegim

Observem que ara data inclou les sèries de poblacions. Anem-ho a representar:

import pylab as p
anys= data[:,0] # indiquem aixi la primera columna. : indica tots els valors.
olmeda= data[:,1]
pozuelo= data[:,2]
p.plot(anys,pozuelo,'ro')
p.show()
p.plot(anys,olmeda,'bo')
p.show()

Per veure-ho millor anem a fer una gràfica més sofisticada amb la comanda figure. Com que hi haurà més instruccions anem a editar-ho tot en un fitxer de text. Fem:

! gedit poblacio.py &

o el nostre editor preferit. En el fitxer insertem el següent:

import pylab as p
from numpy import *
from scipy.io import write_array
from scipy.io import read_array
data = read_array(file("dades.csv"), columns=(0,1,2))
anys= data[:,0] # indiquem aixi la primera columna. : indica tots els valors.
olmeda= data[:,1]
pozuelo= data[:,2]
# Iniciem una nova figura
p.figure(1)
p.subplot(211)
p.title("Pozuelo de Alarcon")
p.plot(anys,pozuelo,'ro')

# Aquest espai es important
p.subplot(212)
p.title("Olmeda de la Sierra")
p.plot(anys,olmeda,'bo')
p.show()

Un cop l’hem guardat al mateix directori on estem treballant fem:

execfile("poblacio.py")

2.2 Ajustament exponencial d’algunes gràfiques

En aquesta secció veurem com ajustar una d’aquestes poblacions, per exemple la de Pozuelo de Alarcón, a una exponencial i a veure’n l’error. Per això usarem la instrucció polyfit que és semblant a la similar del Matlab.

logpoblacio=log(pozuelo) #calculem el logaritme.
(ar,br)=polyfit(anys,logpoblacio,1) # ajustem per polinomi de grau 1
# aixi pozuelo= exp(ar*poblacio+br)
xr=polyval([ar,br],anys) # Avaluem el polinomi en la regressio
pozueloajustat=exp(xr) # Aquest es el nostre valor ajustat.

print('Ajustament logistic usant el polyfit')
print('regressio: a=%.2f b=%.2f' % (ar,br))

#Anem a dibuixar-ho
p.title('Ajustament de la poblacio de pozuelo')
p.plot(anys,pozuelo,'r.--')
p.plot(anys,pozueloajustat,'g.-')
p.legend(['Poblacio real', 'ajustament'])
# afegim una llegenda.
p.show()

Observem com el creixement exponencial és molt més acusat després de 1940 ja que entre 1930 i 1940 decreix (sobren les explicacions). Anem a mirar d’ajustar-ho per a aquests anys. Per això només cal que triem convenientment les dades

import pylab as p
from numpy import *
from scipy.io import write_array
from scipy.io import read_array
data = read_array(file("dades.csv"),columns=(0,1,2))
anys= data[4:,0] # indiquem aixi la primera columna. 4: indica \ge 4
pozuelo= data[4:,2]
logpoblacio=log(pozuelo) #calculem el logaritme.
(ar,br)=polyfit(anys,logpoblacio,1) # ajustem per polinomi de grau 1
# aixi pozuelo= exp(ar*poblacio+br)
xr=polyval([ar,br],anys) # Avaluem el polinomi en la regressio
pozueloajustat=exp(xr) # Aquest es el nostre valor ajustat.

print('Ajustament logistic usant el polyfit')
print('regressio: a=%.2f b=%.2f' % (ar,br))

#Anem a dibuixar-ho
p.title('Ajustament de la poblacio de pozuelo')
p.plot(anys,pozuelo,'r.--')
p.plot(anys,pozueloajustat,'g.-')
p.legend(['Poblacio real', 'ajustament'])
# afegim una llegenda.
p.show()

Per calcular l’error podem fer-ho un xic millor usant el paquet stats del scipy:

from scipy import stats
(a_s,b_s,r,pvalor,estderr)=stats.linregress(anys,logpoblacio)
# a_s i b_s son els nous parametres de regressio
# r valor r (quadrat del coeficient de correlacio)
# pvalor (doncs aixo)
# estderr error estandard
#
#
print('Regressio lineal usant scipy.stats')
print('Parametres regressio: a=%.2f b=%.2f, std error= %.3f' % (a_s,b_s,estderr))

Comproveu que l’error és petit i l’ajustament molt bo (mireu el p-valor). També ho podríeu fer amb rpy.

2.3 Exercicis

Trieu una població en què pogueu usar les dades de l’INE i feu el mateix ajustament. Quina és la població que s’aproxima més a una exponecial que creix i a una que decreix?
(*) El creixement d’articles a wikipedia sembla ser també exponencial. Consulteu per exemple la pàgina:

http://en.wikipedia.org/wiki/Wikipedia:Size_of_WikipediaWikipedia

Mireu de reproduir els càlculs de la figura. A l’hora de dibuixar us poden anar bé les gràfiques d’exemple de:

http://matplotlib.sourceforge.net/screenshots.htmlMatplotlib Screeenshots