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

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!

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