Hoy cambiamos tópico, veremos como se puede hacer mapas en R. Queremos tratar de presentar como una mapa el contenido del página de wikipedia https://en.wikipedia.org/wiki/List_of_countries_by_real_population_density_based_on_food_growing_capacity un listado de países por su densidad de población real, es decir, no en su relación al area (que esta contando vinedos y olivares iguales como desiertos y montañas), pero en su relación a la capacidad de producir alimentos. Saltamos sobre las problemas en como medir eso, y concentramos en su presentación.
Primero, como importar los datos de wikipedia en una hoja de cálculo, que podemos despues leer en R:
Abrir una cuenta en google docs, si no lo tienes, y abrir alla una nueva hoja de cálculo. En la celda A1 tienes que escribir:
=ImportHtml("https://en.wikipedia.org/wiki/List_of_countries_by_real_population_density_based_on_food_growing_capacity","table",5)
Los tres argumentos a la función ImportHtml es primero, el url que queremos importar, segundo, que queremos importar de alla una table, y tercero, que la tabla que nos interesa es el quinto en de esa página. Tenemos que experimentar un poco para encontrar que cinco es el número correcto.
Eso esta hecho, y la table que yo he importado así esta publicado ahora con la url:
https://docs.google.com/spreadsheet/pub?key=0AmA1eTG-BTE9dGI1N0IyVTVOZEk2Ri10N0xkbzFRS3c&single=true&gid=0&output=csv
(para obtener este enlace, hay que especificar que el formato publicado debe ser como csv, el default es como página web!)
que podemos usar ahora para importarlo en R, usando read.csv (o talves read.csv2).
En R:
> require(RCurl)
Loading required package: RCurl
Loading required package: bitops
Note: no visible binding for global variable '.Data'
Note: no visible binding for global variable '.Data'
> myURL <- getURL("https://docs.google.com/spreadsheet/pub?key=0AmA1eTG-BTE9dGI1N0IyVTVOZEk2Ri10N0xkbzFRS3c&single=true&gid=0&output=csv")
> DR <- read.csv(textConnection(myURL))
> dim(DR)
[1] 234 8
> summary(DR)
Rank Country Population..July.2005.est..
- : 1 Afghanistan : 1 10006835: 1
1 : 1 Albania : 1 10079380: 1
10 : 1 Algeria : 1 10241138: 1
100 : 1 * American Samoa (US)*: 1 103176 : 1
101 : 1 Andorra : 1 10364388: 1
102 : 1 Angola : 1 1041806 : 1
(Other):228 (Other) :228 (Other) :228
Land.Area..kmÂ.. X..of.arable.land..2005.est.. Arable.Land..kmÂ..
102 : 2 0% : 15 0 : 15
0.44 : 1 20% : 4 20 : 8
10 : 1 1.64% : 3 10 : 6
10000 : 1 5.71% : 3 40 : 5
1001 : 1 10% : 2 30 : 3
100250 : 1 16.67% : 2 53 : 3
(Other):227 (Other):205 (Other):194
Population.Density..pop.per.kmÂ..
3 : 9
14 : 5
12 : 3
141 : 3
15 : 3
17 : 3
(Other):208
Real.Population.Density..pop.per.kmÂ..of.arable.land.
- : 15
381 : 3
1434 : 2
219 : 2
225 : 2
248 : 2
(Other):208
> colnames(DR)
[1] "Rank"
[2] "Country"
[3] "Population..July.2005.est.."
[4] "Land.Area..kmÂ.."
[5] "X..of.arable.land..2005.est.."
[6] "Arable.Land..kmÂ.."
[7] "Population.Density..pop.per.kmÂ.."
[8] "Real.Population.Density..pop.per.kmÂ..of.arable.land."
Bién, antes del análisis talves debemos camiar algunos nombres!
(Encontramos un problema en trabajar de esta manera: Editar los datos en la hoja de google aparwece como dificil, porque todas las celdas estan como formulas, no como valores, y no encontré una manera de cambiar esto. Por esa razon, más abajo reimportamos los datos, primero exportando la hoja de google como una hoja excel, y editar eso en la maquina local).
Ahora, para dibujar una mapa mundi:
> require(maps)
Loading required package: maps
> require(mapdata)
Loading required package: mapdata
map("worldHires", xlim=c(-170, 180), ylim=c(-70, 70), col="grey90", fill=TRUE)
Ahora rellemos los datos, despues de bajar la hjoja electrónca a un fichero local, y editar un poco:
DR <- read.csv("wikip_listado_paises_dens_pob_real.csv",header=TRUE, colClasses=c("numeric", "character", rep("numeric", 6)))
> colnames(DR) <- c("Rank", "Country", "Pop", "Area","X.Arable" , "Arable", "Pop.dens", "Real.Pop.dens")
Para hacer el mapa mundi temático usamos el paquete rworldmap:
> require(rworldmap)
Loading required package: rworldmap
### Welcome to rworldmap ###
For a short introduction type : vignette('rworldmap')
> spDR <- joinCountryData2Map(DR, joinCode="NAME", nameJoinColumn="Country")
191 codes from your data successfully matched countries in the map
42 codes from your data failed to match with a country code in the map
53 codes from the map weren't represented in your data
mapCountryData(spDR, nameColumnToPlot="Real.Pop.dens", mapTitle="Real Population Density", catMethod="quantiles", colourPalette="heat")
Y aqui se debe ver, que paises/regiones tiene realmente baja densidad de población; Norteamerica, argentina, Australia, Rusia, alto: grandes partes de asia. Grandes partes de Africa es intermedio: con excepción de pocos países, sobrepoblación no parece ser un problema de Africa.
Estadística, programacíon y cosas afines...
Enlaces, ideas y ejemplos ...
sábado, 5 de octubre de 2013
jueves, 26 de septiembre de 2013
Continuamos explorando las posibilidades de hacer plots en julia. Hoy el paquete Gaston, que implementa un interfaz a gnuplot, que tiene que ser instalado en el sistema separadamente.
julia> Pkg.add("Gaston")
julia> using Gaston
julia> gaston_demo()
WARNING: contains(collection, item) is deprecated, use in(item, collection) instead
in contains at reduce.jl:238
Se abre rápidamente muchos ventanas de gnuplot con plots. Mostramos unos ejemplos, tomados del manual de Gaston.
julia> t = 0:0.0001:.15
0.0:9.999999999999999e-5:0.15
julia> length(t)
1501
julia> carrier = cos(2pi*t*200) ;
julia> modulator = 0.7+0.5*cos(2pi*t*15) ;
julia> am = carrier.*modulator ;
julia> plot(t,am,"color","blue","legend","AM DSB-SC","linewidth",1.5,
t,modulator,"color","black","legend","Envelope",
t,-modulator,"color","black","title","AM DSB-SC example",
"xlabel","Time (s)","ylabel","Amplitude",
"box","horizontal top left")
0
julia> set_filename("ejemplo_plot_2.png")
"ejemplo_plot_2.png"
julia> printfigure("png")
1
Y ahora el plot esta escrito a fichero, y aparece como abajo:
julia> Pkg.add("Gaston")
julia> using Gaston
julia> gaston_demo()
WARNING: contains(collection, item) is deprecated, use in(item, collection) instead
in contains at reduce.jl:238
Se abre rápidamente muchos ventanas de gnuplot con plots. Mostramos unos ejemplos, tomados del manual de Gaston.
julia> t = 0:0.0001:.15
0.0:9.999999999999999e-5:0.15
julia> length(t)
1501
julia> carrier = cos(2pi*t*200) ;
julia> am = carrier.*modulator ;
julia> plot(t,am,"color","blue","legend","AM DSB-SC","linewidth",1.5,
t,modulator,"color","black","legend","Envelope",
t,-modulator,"color","black","title","AM DSB-SC example",
"xlabel","Time (s)","ylabel","Amplitude",
"box","horizontal top left")
0
julia> set_filename("ejemplo_plot_2.png")
"ejemplo_plot_2.png"
julia> printfigure("png")
1
Y ahora el plot esta escrito a fichero, y aparece como abajo:
Continuamos con otro ejemplo tomado del manual, un histograma de datos simulados de une densidad jicuadrada con dos grados de libertad (el man ual lo llama "densidad Rayleigh"):
julia> y = sqrt( randn(10000).^2 + randn(10000).^2 ) ;
julia> histogram(y,"bins",25,"norm",1,"color","blue","title","Rayleigh density (25 bins)")
1
julia> set_filename("ejemplo_plot_3.png")
"ejemplo_plot_3.png"
julia> printfigure("png")
1
Hoy tratarémos de hacer plots en julia. Existen varios paquetes para hacer plots, pero esta muy temprano en su desarrollo por el momento. Primero veremos a Winston, que implemanta plots en el estili de octave (o Matlab...).
julia> using Winston
Warning: using Winston.plot in module Main conflicts with an existing identifier.
julia> x = [-pi:pi/60:pi];
julia> length(x)
121
julia> y=sin(x);
Terminar la line de comando con ; inhibe que se imprime todo el resultado (como en octave). Pero en julia no es necesario usar ; en las definiciones de funciones, solo en uso interactivo.
julia> p = FramedPlot()
FramedPlot(...)
julia> add(p, Curve(x,y))
julia> setattr(p, "xlabel", "x")
"x"
julia> setattr(p, "ylabel", "\\sin(x)")
"\\sin(x)"
julia> setattr(p, "title", "curva de sinus")
"curva de sinus"
julia> file(p,"plot_ejemplo_1.png")
Por el momento parece que no hay manera de construir el plot interactivamente, es necesrio escribirlo a un fichero como hemos hecho arriba.
Para uso estadístico, sin embargo, el paquete Gadfly parece más interesante, como implementa el "Grammar of Graphics". Pero por el momento solo me da mensajes de error ...
julia> using Winston
Warning: using Winston.plot in module Main conflicts with an existing identifier.
julia> x = [-pi:pi/60:pi];
julia> length(x)
121
julia> y=sin(x);
Terminar la line de comando con ; inhibe que se imprime todo el resultado (como en octave). Pero en julia no es necesario usar ; en las definiciones de funciones, solo en uso interactivo.
julia> p = FramedPlot()
FramedPlot(...)
julia> add(p, Curve(x,y))
julia> setattr(p, "xlabel", "x")
"x"
julia> setattr(p, "ylabel", "\\sin(x)")
"\\sin(x)"
julia> setattr(p, "title", "curva de sinus")
"curva de sinus"
julia> file(p,"plot_ejemplo_1.png")
Para uso estadístico, sin embargo, el paquete Gadfly parece más interesante, como implementa el "Grammar of Graphics". Pero por el momento solo me da mensajes de error ...
martes, 24 de septiembre de 2013
julia tiene un sistema de paquetes de extensión, y para uso para análisis de datos tenemos el
paquete DataFrames, que se base en data.frames como originalmente implementado en R/S. Para usarlo tenemos que escribir:
julia> using DataFrames
paquete DataFrames, que se base en data.frames como originalmente implementado en R/S. Para usarlo tenemos que escribir:
julia> using DataFrames
y ahora tenemos acceso a las definiciones del paquete. Primero, implementa un valor especial NA "no accesible" para representar valores de variables desconocidos:
julia> NA
NA
julia> NA+1
NA
julia> sin(NA)
NA
Hacemos un vector de datos: (el tipo DataArray)
julia> dv =DataArray([1,4,2,9,4,5,2])
7-element DataArray{Int64,1}:
1
4
2
9
4
5
2
julia> length(dv)
7
julia> size(dv)
(7,)
Podemos asignar NA a uno de los elementos:
julia> dv[2]=NA
NA
julia> dv'
1x7 DataArray{Int64,2}:
1 NA 2 9 4 5 2
Note el uso de ' para transposición! (Y aqui, para ahorrar espacio).
julia> mean(dv)
NA
¿Podemos calcular el promedio de los elementos que no es NA? Usamos la función
removeNA:
julia> removeNA(dv)
6-element Array{Int64,1}:
1
2
9
4
5
2
julia> mean(removeNA(dv))
3.8333333333333335
Ahora construimos un DataFrame, una cierta matriz de filas y coumnas, donde cada columna es un
DataArray:
julia> df=DataFrame(A=1:5, B=["M","M","H","H","M"])
5x2 DataFrame:
A B
[1,] 1 "M"
[2,] 2 "M"
[3,] 3 "H"
[4,] 4 "H"
[5,] 5 "M"
julia> describe(df)
A
Min 1.0
1st Qu. 2.0
Median 3.0
Mean 3.0
3rd Qu. 4.0
Max 5.0
NAs 0
NA% 0.0%
B
Length 5
Type ASCIIString
NAs 0
NA% 0.0%
Unique 2
julia> size(df)
(5,2)
julia> length(df)
2
Veremos unos conjuntos de datos clasicos contenido en R:
julia> using RDatasets
julia> anscombe = data("datasets","anscombe")
11x9 DataFrame:
x1 x2 x3 x4 y1 y2 y3 y4
[1,] 1 10 10 10 8 8.04 9.14 7.46 6.58
[2,] 2 8 8 8 8 6.95 8.14 6.77 5.76
[3,] 3 13 13 13 8 7.58 8.74 12.74 7.71
[4,] 4 9 9 9 8 8.81 8.77 7.11 8.84
[5,] 5 11 11 11 8 8.33 9.26 7.81 8.47
[6,] 6 14 14 14 8 9.96 8.1 8.84 7.04
[7,] 7 6 6 6 8 7.24 6.13 6.08 5.25
[8,] 8 4 4 4 19 4.26 3.1 5.39 12.5
[9,] 9 12 12 12 8 10.84 9.13 8.15 5.56
[10,] 10 7 7 7 8 4.82 7.26 6.42 7.91
[11,] 11 5 5 5 8 5.68 4.74 5.73 6.89
julia> describe(anscombe)
Min 1.0
1st Qu. 3.5
Median 6.0
Mean 6.0
3rd Qu. 8.5
Max 11.0
NAs 0
NA% 0.0%
x1
Min 4.0
1st Qu. 6.5
Median 9.0
Mean 9.0
3rd Qu. 11.5
Max 14.0
NAs 0
NA% 0.0%
x2
Min 4.0
1st Qu. 6.5
Median 9.0
Mean 9.0
3rd Qu. 11.5
Max 14.0
NAs 0
NA% 0.0%
x3
Min 4.0
1st Qu. 6.5
Median 9.0
Mean 9.0
3rd Qu. 11.5
Max 14.0
NAs 0
NA% 0.0%
x4
Min 8.0
1st Qu. 8.0
Median 8.0
Mean 9.0
3rd Qu. 8.0
Max 19.0
NAs 0
NA% 0.0%
y1
Min 4.26
1st Qu. 6.3149999999999995
Median 7.58
Mean 7.500909090909093
3rd Qu. 8.57
Max 10.84
NAs 0
NA% 0.0%
y2
Min 3.1
1st Qu. 6.695
Median 8.14
Mean 7.500909090909091
3rd Qu. 8.95
Max 9.26
NAs 0
NA% 0.0%
y3
Min 5.39
1st Qu. 6.25
Median 7.11
Mean 7.500000000000001
3rd Qu. 7.98
Max 12.74
NAs 0
NA% 0.0%
y4
Min 5.25
1st Qu. 6.17
Median 7.04
Mean 7.50090909090909
3rd Qu. 8.190000000000001
Max 12.5
NAs 0
NA% 0.0%
julia> size(anscombe)
(11,9)
julia> length(anscombe)
9
julia> typeof(anscombe)
DataFrame (constructor with 22 methods)
julia> colnames(anscombe)
9-element Array{Union(ASCIIString,UTF8String),1}:
""
"x1"
"x2"
"x3"
"x4"
"y1"
"y2"
"y3"
"y4"
Eso es "anscombes quartet", un conjunto de data clasico, construido para mostrar que pares de variables
(x,y) puede tener exactamente los mismos valores de muchas estadisticos descriptivos, como correlación, coeficientes de regresión, etc, pero, sim embargo, tener plots muy diferentes!
Podemos investigar eso con julia?
julia> using GLM
julia> using Distributions
julia> colnames(anscombe)
9-element Array{Union(ASCIIString,UTF8String),1}:
""
"x1"
"x2"
"x3"
"x4"
"y1"
"y2"
"y3"
"y4"
julia> a1 = lm(:(y1 ~ x1), anscombe)
Formula: y1 ~ x1
Coefficients:
2x4 DataFrame:
Estimate Std.Error t value Pr(>|t|)
[1,] 3.00009 1.12475 2.66735 0.0257341
[2,] 0.500091 0.117906 4.24146 0.00216963
julia> a2 = lm(:(y2 ~ x2), anscombe)
Formula: y2 ~ x2
Coefficients:
2x4 DataFrame:
Estimate Std.Error t value Pr(>|t|)
[1,] 3.00091 1.1253 2.66676 0.0257589
[2,] 0.5 0.117964 4.23859 0.00217882
Aqui, la "Formula" en :(y2 ~ x2) se interpreta como en R.
En otro blog veremos como hacer plots!
viernes, 20 de septiembre de 2013
Continuamos el estudio de julia! julia tiene un sistema de paquetes, como R. Tratarémos de instalar el paquete para modelos lineales generalizados, , GLM:
julia> Pkg.add("GLM")
INFO: Initializing package repository /home/kjetil/.julia.
INFO: Cloning METADATA from git://github.com/JuliaLang/METADATA.jl
INFO: Cloning cache of Blocks from git://github.com/tanmaykm/Blocks.jl.git
INFO: Cloning cache of GLM from git://github.com/JuliaStats/GLM.jl.git
INFO: Cloning cache of DataFrames from git://github.com/JuliaStats/DataFrames.jl.git
INFO: Cloning cache of Distributions from git://github.com/JuliaStats/Distributions.jl.git
INFO: Cloning cache of GZip from git://github.com/kmsquire/GZip.jl.git
INFO: Cloning cache of NumericExtensions from https://github.com/lindahua/NumericExtensions.jl.git
WARNING: gnome-keyring:: couldn't connect to: /run/user/kjetil/keyring-q2W4Tr/pkcs11: No such file or directory
INFO: Cloning cache of Stats from git://github.com/JuliaStats/Stats.jl.git
INFO: Installing Blocks v0.0.0
INFO: Installing GLM v0.2.0
INFO: Installing DataFrames v0.3.11
INFO: Installing Distributions v0.2.7
INFO: Installing GZip v0.2.4
INFO: Installing NumericExtensions v0.2.13
INFO: Installing Stats v0.2.7
INFO: REQUIRE updated.
julia> Pkg.add("GLM")
INFO: Initializing package repository /home/kjetil/.julia.
INFO: Cloning METADATA from git://github.com/JuliaLang/METADATA.jl
INFO: Cloning cache of Blocks from git://github.com/tanmaykm/Blocks.jl.git
INFO: Cloning cache of GLM from git://github.com/JuliaStats/GLM.jl.git
INFO: Cloning cache of DataFrames from git://github.com/JuliaStats/DataFrames.jl.git
INFO: Cloning cache of Distributions from git://github.com/JuliaStats/Distributions.jl.git
INFO: Cloning cache of GZip from git://github.com/kmsquire/GZip.jl.git
INFO: Cloning cache of NumericExtensions from https://github.com/lindahua/NumericExtensions.jl.git
WARNING: gnome-keyring:: couldn't connect to: /run/user/kjetil/keyring-q2W4Tr/pkcs11: No such file or directory
INFO: Cloning cache of Stats from git://github.com/JuliaStats/Stats.jl.git
INFO: Installing Blocks v0.0.0
INFO: Installing GLM v0.2.0
INFO: Installing DataFrames v0.3.11
INFO: Installing Distributions v0.2.7
INFO: Installing GZip v0.2.4
INFO: Installing NumericExtensions v0.2.13
INFO: Installing Stats v0.2.7
INFO: REQUIRE updated.
Peo ... hasta hoy no existe mucho documentación sobre como actualmente _usar_ loe métodos de GLM, entonces regresarémos a este tópico otro día!
Ha aparecido un nuevo lenguage de programación dinámica para cosas numericos y estadísticos. julia, mire http://julialang.org/
Parece muy interesante, con sintaxis parecido a matlab/octave pero velocidad de c! Comenzamos a explorarlo aqui:
kjetil@kjetil-HP-Pavilion-dv4-Notebook-PC:~$ julia
_
_ _ _(_)_ | A fresh approach to technical computing
(_) | (_) (_) | Documentation: http://docs.julialang.org
_ _ _| |_ __ _ | Type "help()" to list help topics
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 0.2.0
_/ |\__'_|_|_|\__'_| |
|__/ | x86_64-linux-gnu
julia> 2+2
4
julia> x=2
2
julia> y=4x^3-2x+1
29
Parece muy interesante, con sintaxis parecido a matlab/octave pero velocidad de c! Comenzamos a explorarlo aqui:
kjetil@kjetil-HP-Pavilion-dv4-Notebook-PC:~$ julia
_
_ _ _(_)_ | A fresh approach to technical computing
(_) | (_) (_) | Documentation: http://docs.julialang.org
_ _ _| |_ __ _ | Type "help()" to list help topics
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 0.2.0
_/ |\__'_|_|_|\__'_| |
|__/ | x86_64-linux-gnu
julia> 2+2
4
julia> x=2
2
julia> y=4x^3-2x+1
29
Mire el sintaxis interesante, multiplicación con un número literal sin *.
Soporte muy directo de matrices, con notación natural:
julia> B = [1 0 0 ; 1 2 3 ; 0 3 5]
3x3 Array{Int64,2}:
1 0 0
1 2 3
0 3 5
julia> rank(B)
3
julia> B * B^-1
3x3 Array{Float64,2}:
1.0 0.0 0.0
0.0 1.0 8.88178e-16
3.55271e-15 -3.55271e-15 1.0
Pero se ve que la implementación hasta hoy falta en robustez:
julia> A = [1 2 3; 2 3 4;3 4 5]
3x3 Array{Int64,2}:
1 2 3
2 3 4
3 4 5
julia> size(A)
(3,3)
julia> A^-1
3x3 Array{Float64,2}:
-4.5036e15 9.0072e15 -4.5036e15
9.0072e15 -1.80144e16 9.0072e15
-4.5036e15 9.0072e15 -4.5036e15
julia> A^-1 * A
3x3 Array{Float64,2}:
0.0 0.0 0.0
-4.0 0.0 8.0
2.0 0.0 0.0
julia> rank(A)
2
julia devuelve un inverso para A aun que no existe! Esperamos que se corrige esto .
julia> whos()
A 3x3 Array{Int64,2}
B 3x3 Array{Int64,2}
Base Module
Core Module
Main Module
x Int64
y Int64
Observamos que whos() da la información que whos
da en matlab/octave. Pero, es una función, no un comando.
julia> tri = [1 2 3; 0 4 5;0 0 6]
3x3 Array{Int64,2}:
1 2 3
0 4 5
0 0 6
julia> tri'
3x3 Array{Int64,2}:
1 0 0
2 4 0
3 5 6
Notación natural para la transpuesta, como en octave.
Finalizamos hoy con unos ejemplos de estadística:
julia> x = [ 1., 3.1, 2.2, 5, 4, 2, 6, 8, 7, 9, 7, 8, 6.9, 5]
14-element Array{Float64,1}:
1.0
3.1
2.2
5.0
4.0
2.0
6.0
8.0
7.0
9.0
7.0
8.0
6.9
5.0
julia> mean(x)
5.3
julia> std(x)
2.52373349806012
julia> median(x)
5.5
julia> quantile(x, [ .22, .44, .89])
3-element Array{Float64,1}:
2.974
5.0
8.0
Suscribirse a:
Entradas (Atom)