sábado, 5 de octubre de 2013

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.

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:


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 ...



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

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.

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
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