Hoy vamos a ver como podemos juntar lo bueno de R, algunas de sus librerías, con Python usando rpy2.
Pero, lo primero de todo, ¿qué es rpy2? rpy2 es una interfaz que permite que podamos comunicar información entre R y Python y que podamos acceder a funcionalidad de R desde Python. Por tanto, podemos estar usando Python para todo nuestro análisis y en el caso de que necesitemos alguna librería estadística especializada de R podremos acceder a la misma usando rpy2.
Para poder usar rpy2 necesitarás tener instalado tanto Python (CPython versión >= 2.7.x) como R (versión >=3), además de las librerías R a las que quieras acceder. Conda permite realizar todo el proceso de instalación de los intérpretes de Python y R, además de librerías, pero no he trabajado con Conda y R por lo que no puedo aportar mucho más en este aspecto. Supongo que será parecido a lo que hacemos con Conda y Python.
Para este microtutorial voy a hacer uso de la librería extRemes de R que permite hacer análisis de valores extremos usando varias de las metodologías más comúnmente aceptadas.
Como siempre, primero de todo, importaremos la funcionalidad que necesitamos para la ocasión.
In [1]:
# Importamos pandas y numpy para manejar los datos que pasaremos a R
import pandas as pd
import numpy as np
# Usamos rpy2 para interactuar con R
import rpy2.robjects as ro
# Activamos la conversión automática de tipos de rpy2
import rpy2.robjects.numpy2ri
rpy2.robjects.numpy2ri.activate()
import matplotlib.pyplot as plt
%matplotlib inline
En el anterior código podemos ver una serie de cosas nuevas que voy a explicar brevemente:
import rpy2.robjects as ro, esto lo explicaremos un poquito más abajo.import rpy2.robjects.numpy2ri, importamos el módulo numpy2ri. Este módulo permite que hagamos conversión automática de objetos numpy a objetos rpy2.rpy2.robjects.numpy2ri.activate(), hacemos uso de la función activate que activa la conversión automática de objetos que hemos comentado en la línea anterior.Para evaluar directamente código R podemos hacerlo usando rpy2.robjects.r con el código R expresado como una cadena (rpy2.robjects lo he importado como ro en este caso, como podéis ver más arriba):
In [2]:
codigo_r = """
saluda <- function(cadena) {
return(paste("Hola, ", cadena))
}
"""
ro.r(codigo_r)
Out[2]:
En la anterior celda hemos creado una función R llamada saluda y que ahora está disponible en el espacio de nombres global de R. Podemos acceder a la misma desde Python de la siguiente forma:
In [3]:
saluda_py = ro.globalenv['saluda']
Y podemos usarla de la siguiente forma:
In [4]:
res = saluda_py('pepe')
print(res[0])
En la anterior celda véis que para acceder al resultado he tenido que usar res[0]. En realidad, lo que nos devuleve rpy2 es:
In [5]:
print(type(res))
print(res.shape)
En este caso un numpy array con diversa información del objeto rpy2. Como el objeto solo devuelve un string pues el numpy array solo tiene un elemento.
Podemos acceder al código R de la función de la siguiente forma:
In [6]:
print(saluda_py.r_repr())
Hemos visto como acceder desde Python a nombres disponibles en el entorno global de R. ¿Cómo podemos hacer para que algo que creemos en Python este accesible en R?
In [7]:
variable_r_creada_desde_python = ro.FloatVector(np.arange(1,5,0.1))
Veamos como es esta variable_r_creada_desde_python dentro de Python
In [8]:
variable_r_creada_desde_python
Out[8]:
¿Y lo que se tendría que ver en R?
In [9]:
print(variable_r_creada_desde_python.r_repr())
Pero ahora mismo esa variable no está disponible desde R y no la podríamos usar dentro de código R que permanece en el espacio R (vaya lío, ¿no?)
In [10]:
ro.r('variable_r_creada_desde_python')
Vale, tendremos que hacer que sea accesible desde R de la siguiente forma:
In [11]:
ro.globalenv["variable_ahora_en_r"] = variable_r_creada_desde_python
print(ro.r("variable_ahora_en_r"))
Ahora que ya la tenemos accesible la podemos usar desde R. Por ejemplo, vamos a usar la función sum en R que suma los elementos pero directamente desde R:
In [12]:
print(ro.r('sum(variable_ahora_en_r)'))
print(np.sum(variable_r_creada_desde_python))
Perfecto, ya sabemos, de forma muy sencilla y básica, como podemos usar R desde Python, como podemos pasar información desde R hacia Python y desde Python hacia R. ¡¡¡Esto es muy poderoso!!!, estamos juntando lo mejor de dos mundos, la solidez de las herramientas científicas de Python con la funcionalidad especializada que nos pueden aportar algunas librerías de R no disponibles en otros ámbitos.
Vamos a empezar importando la librería extRemes de R:
In [13]:
# Importamos la librería extRemes de R
from rpy2.robjects.packages import importr
extremes = importr('extRemes')
En la anterior celda hemos hecho lo siguiente:
from rpy2.robjects.packages import importr, La función importr nos servirá para importar las librerías Rextremes = importr('extRemes'), de esta forma importamos la librería extRemes de R, sería equivalente a hacer en R library(extRemes).Leemos datos con pandas. En el mismo repo donde está este notebook está también un fichero de texto con datos que creé a priori. Supuestamente son datos horarios de velocidad del viento por lo que vamos a hacer análisis de valores extremos de velocidad del viento horaria.
In [15]:
data = pd.read_csv('datasets/Synthetic_data.txt',
sep = '\s*', skiprows = 1, parse_dates = [[0, 1]],
names = ['date','time','wspd'], index_col = 0)
In [16]:
data.head(3)
Out[16]:
Extraemos los máximos anuales los cuales usaremos posteriormente dentro de R para hacer cálculo de valores extremos usando la distribución generalizada de valores extremos (GEV):
In [17]:
max_y = data.wspd.groupby(pd.TimeGrouper(freq = 'A')).max()
Dibujamos los valores máximos anuales usando Pandas:
In [18]:
max_y.plot(kind = 'bar', figsize = (12, 4))
Out[18]:
Referenciamos la funcionalidad fevd (fit extreme value distribution) dentro del paquete extremes de R para poder usarla directamente con los valores máximos que hemos obtenido usando Pandas y desde Python.
In [19]:
fevd = extremes.fevd
Como hemos comentado anteriormente, vamos a calcular los parámetros de la GEV usando el método de ajuste GMLE (Generalised Maximum Lihelihood Estimation) y los vamos a guardar directamente en una variable Python.
Veamos la ayuda antes:
In [20]:
print(fevd.__doc__)
Y ahora vamos a hacer un cálculo sin meternos mucho en todas las opciones posibles.
In [21]:
res = fevd(max_y.values, type = "GEV", method = "GMLE")
¿Qué estructura tiene la variable res que acabamos de crear y que tiene los resultados del ajuste?
In [22]:
print(type(res))
In [23]:
print(res.r_repr)
Según nos indica lo anterior, ahora res es un vector que está compuesto de diferentes elementos. Los vectores pueden tener un nombre para todos o algunos de los elementos. Para acceder a estor nombres podemos hacer:
In [24]:
res.names
Out[24]:
Según el output anterior, parece que hay un nombre results, ahí es donde se guardan los valores del ajuste, los estimadores. Para acceder al mismo podemos hacerlo de diferentes formas. Con Python tendriamos que saber el índice y acceder de forma normal (__getitem__()). Existe una forma alternativa usando el método rx que nos permite acceder directamente con el nombre:
In [25]:
results = res.rx('results')
In [26]:
print(results.r_repr)
Parece que tenemos un único elemento:
In [27]:
results = results[0]
results.r_repr
Out[27]:
Vemos ahora que results tiene un elemento con nombre par donde se guardan los valores de los estimadores del ajuste a la GEV que hemos obtenido usando GMLE. Vamos a obtener finalmente los valores de los estimadores:
In [28]:
location, scale, shape = results.rx('par')[0][:]
print(location, scale, shape)
Usamos la antigua función mágica rmagic que ahora se activará en el notebook de la siguiente forma:
In [29]:
%load_ext rpy2.ipython
Veamos como funciona la functión mágica de R:
In [30]:
help(rpy2.ipython.rmagic.RMagics.R)
A veces, será más simple usar la función mágica para interactuar con R. Veamos un ejemplo donde le pasamos a R el valor obtenido de la función fevd del paquete extRemes de R que he usado anteriormente y corremos cierto código directamente desde R sin tener que usar ro.r.
In [31]:
%R -i res plot.fevd(res)
En la anterior celda de código le he pasado como parámetro de entrada (- i res) la variable res que había obtenido anteriormente para que esté disponible desde R. y he ejecutado código R puro (plot.fevd(res)).
Si lo anterior lo quiero hacer con rpy2 puedo hacer lo siquiente:
CUIDADO, la siguiente celda de código puede provocar que se reinicialice el notebook y se rompa la sesión. Si has hecho cambios en el notebook guárdalos antes de ejecutar la celda, por lo que pueda pasar...
In [35]:
ro.globalenv['res'] = res
ro.r("plot.fevd(res)")
Out[35]:
Lo anterior me bloquea el notebook y me 'rompe' la sesión (en windows, al menos) ya que la ventana de gráficos se abre de forma externa... Por tanto, una buena opción para trabajar de forma interactiva con Python y R de forma conjunta y que no se 'rompa' nada es usar tanto rpy2 como su extensión para el notebook de Jupyter (dejaremos de llamarlo IPython poco a poco).
Vamos a combinar las dos formas de trabajar con rpy2 en el siguiente ejemplo:
In [32]:
metodos = ["MLE", "GMLE"]
tipos = ["GEV", "Gumbel"]
Lo que vamos a hacer es calcular los parámetros del ajuste usando la distribución GEV y Gumbel, que es un caso especial de la GEV. El ajuste lo calculamos usando tanto MLE como GMLE. Además de mostrar los valores resultantes del ajuste para los estimadores vamos a mostrar el dibujo de cada uno de los ajustes y algunos test de bondad. Usamos Python para toda la maquinaria de los bucles, usamos rpy2 para obtener los estimadores y usamos la función mágica de rpy2 para mostrar los gráficos del resultado.
In [33]:
for t in tipos:
for m in metodos:
print('tipo de ajuste: ', t)
print('método de ajuste: ', m)
res = fevd(max_y.values, method = m, type = t)
if m == "Bayesian":
print(res.rx('results')[0][-1][0:-2])
elif m == "Lmoments":
print(res.rx('results')[0])
else:
print(res.rx('results')[0].rx('par')[0][:])
%R -i res plot.fevd(res)
Espero que este microtutorial os valga, al menos, para conocer rpy2 y la potencia que os puede llegar a aportar a vuestros análisis 'pythónicos'. Como resumen: