Feliz Navidad 2016 (en R)

Quisiera desearle a mis lectores una feliz Navidad 2016 ... al estilo R. Vacílense mi arbolito navideño 3d creado usando el paquete plot3D de R:

Árbol de navidad en perspectiva y animado

Aunque esto sería aún otra decoración navideña hecha en R, lo original de esta es que el árbol se muestra en perspectiva tridimensional. Hace un tiempo, yo mismo elaboré un código en R para visualizar un árbol navideño en dos dimensiones (ver esta entrada), el cual a su vez era basado en el código que apareció en el blog Wiekvoet. Sin embargo, esta vez el reto era crear un modelo tridimensional del árbol visualizado en perspectiva con una animación mostrando el árbol en rotación.

Un modelo matemático del árbol de Navidad

Mucho antes del nacimiento de Jesús, la mayoría de las culturas paganas al norte de Europa adoraban al sol y tenían rituales en diciembre para celebrar el renacimiento solar y atraer cosechas abundantes durante el año siguiente. Esta era la significancia simbólica el tronco de Yule y del árbol de Navidad en unas festividades de fuego y luz, abundantes en tonos rojos y guirnaldas verdes. En la Hispanoamérica, el motivo navideño principal según nuestras costumbres siempre ha sido más bien el Nacimiento del Niño Jesús, pero gracias a la transculturación y a la influencia de la cultura anglosajona, el árbol de navidad se ha convertido en otro elemento navideño importante. En todo caso, por mi parte siempre he disfrutado de las luces y las decoraciones cautivadoras del árbol navideño, que me traen recuerdos felices pasados junto a familiares y amigos en una época en la cual olvidamos viejos rencores. Dando rienda suelta al "geek" en mi - y porque evidentemente estoy "ladillao" - quise rendirle un tributo a este espíritu de la Navidad creando un modelo matemático del arbolito navideño.

Matemáticamente hablando, ¿qué es un pino navideño sino un fractal? Un fractal que consiste en un segmento de tronco con dos o tres ramas que se proyectan hacia afuera y una extensión del segmento de tronco que se proyecta en su misma dirección. Las ramas a su vez pueden ser consideradas como segmentos de troncos, cada una con más ramas extendiéndose hacia afuera, con extensiones de la rama proyectándose en la misma dirección que la rama madre y así sucesivamente. En consecuencia, la construcción del árbol navideño es un proceso recursivo. Siendo R un lenguaje de programación descendiente del LISP, la estructura de datos por excelencia para implementar este tipo de objetos recursivos es la lista, la cual desde luego es soportada en R base.

Hay que considerar varios aspectos del pino navideño para que nuestro modelo de árbol se parezca a tal y no a una palmera, un árbol de mango, o un arbolito bonzai. Por una parte, los pinos navideños son característicamente cónicos. Esto significa que el ángulo de separación de las ramas con respecto a los segmentos principales de tronco es inicialmente mayor, pero "se cierra" a medida que el pino se extiende hacia arriba. Para controlar el grado de "exterioridad" o de "estar hacia arriba" de una rama o segmento de tronco, se hace necesario incluir un parámetro de profundidad asociado a cada segmento del árbol. Inicialmente, cuando el árbol comienza con un sólo segmento pequeño y verde, la profundidad de este segmento es 1. A medida que nuevos segmentos se van agregando hacia arriba o el árbol se ramifica hacia afuera, las extensiones o ramas tienen una profundidad decreciente, hasta llegar al tope del pino o sus partes externas, cuando la profundidad se acerca a 0. Este parámetro de profundidad nos será útil para determinar en cuales partes del pino colocaremos los adornos y las luces.

Otro aspecto del modelo que es oportuno mencionar es que en cada iteración del algoritmo recursivo de crecimiento, el grosor de cada sección de tronco aumenta a medida que el algoritmo recorre todo el árbol. Finalmente, cuando ya no hay ramas o extensiones por recorrer porque se ha alcanzado un segmento terminal, se crean brotes nuevos para las ramas y la extensión con un grosor inicial de 1. En consecuencia, es necesario incluir este parámetro de grosor, el cual además será de utilidad para dibujar los segmentos del árbol como trozos lineales de determinado grosor. El grosor de un segmento también nos indica cuando ese segmento será dibujado de color marrón (representando el tronco y las ramas leñosas del pino) o verde oscuro (representando las hojas verdes). A partir de cierto nivel umbral del grosor de un segmento, cambia también su color (otro parámetro más). Por último, tenemos el punto inicial y terminal del segmento lineal de árbol, representados como n-plas en el espacio $\mathbb{R}^2$.

La estructura de datos para el árbolito de navidad sera una lista con los siguientes componentes o "slots":

El punto inicial del segmento de tronco/rama, dado como un vector cuyas coordenadas en $\mathbb{R}^3$ son $(x_0,y_0,z_0)$.
El punto terminal del segmento de tronco/rama, dado como un vector con coordenadas $(x_1,y_1,z_1)$. Ambos puntos determinan el vector dirección de la rama/tronco como $\vec{u}=(x_1-x_0,y_1-y_0,z_1-z_0)$. El vector dirección será utilizado para determinar la dirección de la extensión y las ramas hijas de cada segmento de árbol. También nos ayuda a determinar, en forma lineal parametrizada, los puntos donde el segmento se ramifica.
El parámetro de grosor, el cual también indica el grosor del segmento lineal cuando el árbol se dibuja utilizando la rutina segments3D del paquete plot3D.
El parámetro profundidad que indica cuan externo es el respectivo segmento del árbol, como ya se explicó.
Slots para las hasta tres ramas y una extensión (rama1, rama2, rama3 y extension), los cuales a su vez son listas recursivas como esta. Cuando se crea una rama o extensión como segmento adicional del árbol, estos slots se colocan como NULL para indicar que el segmento respectivo es un segmento terminal del árbol.

Crear la extensión de cada segmento no plantea mayores problemas. Uno simplemente toma como punto inicial de la extensión al punto final del segmento padre. El punto final de la extensión se determina sumándole a su punto inicial un múltiplo del vector dirección $\vec{u}=(x_1-x_0,y_1-y_0,z_1-z_0)$ por un escalar para hacer la extensión un poquito más corta que el tronco/rama original. Crear las ramas es un asunto más difícil: ¿cómo crearlas de modo que salgan del tronco principal en ángulos más o menos uniformes de modo que el árbol no tenga toda sus ramas hacia un lado y su forma sea aproximadamente cónica? ¿Cómo obtener los vectores directrices de esas ramas de modo que aparenten proyectarse hacia afuera desde el tronco? Es preciso echar mano de un poco de trigonometría y álgebra vectorial...

Si $\vec{u}$ es el vector directriz del segmento de tronco considerado, entonces lo que necesitamos es un vector perpendicular a $\vec{u}$ - llamemoslo $\vec{v}$ - de modo que podemos obtener el vector directriz de la rama - $\vec{b}$ - sumando $\vec{u}+\vec{v}$ y multiplicando $\vec{b}$ por un escalar para fijar su longitud a la magnitud deseada (ver Fig. 1). No obstante, considerando que tratamos con el espacio métrico $\mathbb{R}^3$, hay que tomar en cuenta que existen infinitos vectores perpendiculares al vector $\vec{u}$. Más aún, hay todo un plano bidimensional $\mathcal{V}$ compuesto de tales vectores que son todos perpendiculares a $\vec{u}$ (ver Fig. 2). Lo que necesitamos para encontrar $\vec{v}$ es una base ortonormal de dos vectores $\{\vec{v_1},\vec{v_2}\}$ que nos permitirá expresar $\vec{v}$ como combinación lineal de ellos.

Relación geométrica entre el tronco de un árbol y una de sus ramas

Fig. 1 Relación geométrica entre una sección del tronco de un arbol ($\vec{u}$) y una de sus ramas ($\vec{b}$). $\theta$ es el ángulo con el cual $\vec{b}$ se proyecta hacia afuera desde el tronco principal.

Esto resulta ser más fácil de lo previsto. Repasemos: dos vectores se dicen ortogonales (perpendiculares) si su producto interno es nulo. Así que si tenemos un vector, digamos $(x,y,z)$, todo lo que tenemos que hacer es escoger cualquier coordenada no nula e intercambiarla con cualquiera de las restantes dos coordenadas, cambiando el signo de álguna de ellas y fijando la 3era coordenada que no tocamos como 0. Así, por ejemplo, $(-z,0,x)$, $(0,z,-y)$ y $(-y,x,0)$ serían todos vectores perpendiculares a $(x,y,z)$, siempre y cuando sean distintos del vector nulo. Resumiendo, a fin de encontrar una base ortonornmal para el plano perpendicular a $\vec{u}$, escogemos cualquier componente no-nula de este vector, escogemos alguna de las dos componentes restantes e intercambiamos lugares y cambiamos el signo de una de las componentes mientras igualamos la tercera componente a cero. Esto nos dará un vector perpendicular a $\vec{u}$. Para encontrar el segundo vector ortogonal, simplemente repetimos el proceso de arriba pero inercambiando con la primera componente no nula con la componente que habiamos hecho cero para el vector anterior. Finalmente, una vez que tenemos los dos vectores $\vec{v_1}$ and $\vec{v_2}$, debemos convertirlos a vectores unitarios dividiéndolos entre sus respectivas normas. El proceso se ilustra mejor en el código R al final de este post.

El vector v como combinación lineal de v1 y v2

Fig. 2 Encontrando al vector $\vec{v}$ como combianción lineal de la base ortonormal del plano $\mathcal{V}$ y luego utilizando este vector para encontrar al vector directríz $\vec{b}$ de la rama.

Con la base ortonormal del plano perpendicular al tronco del árbol ($\vec{u}$), sólo necesitamos dos componentes escalares $c_1$ y $c_2$ para representar el vector perpendicular $v$ que usaremos para construir la rama, como una combinación lineal de la base ortonormal de $\mathcal{V}$. ¿Cómo procederemos? Bueno, retrocedamos un poco y pensemos en lo que realmente necesitamos. Para las partes "más profundas" del árbol, necesitamos que cada segmento de tronco tenga tres ramificaciones y para las partes más superficiales, necesitamos sólo dos. Queremos que las ramas se distribuyan "uniformemente" alrededor del tronco: en el caso de tres ramas queremos un ángulo de $120^\circ$ entre las ramas (como el símbolo de estrella de 3 puntas de la Mercedes-Benz), en el caso de 2 ramas, un ángulo de separación entre ramos de $180^\circ\pm 20^\circ$ será adecuado.

Esto significa que queremos escoger los vectores $\vec{v}$ en el plano $\mathcal{V}$ de modo que su disposición sea como se muestra en la Fig. 3. Con esto, nos aseguraremos que nuestro pino navideño tenga una forma aproximadamente cónica.


Fig. 3a 3 vectores/ramas	Fig. 3b 2 vectores/ramas
Fig. 3 Disposición de dos o tres vectores en $\mathcal{V}$ para la generación de ramas.

Nuestro problema es ahora generar 2 o 3 vectores en un plano con tales ángulos de separación. La buena noticia es que podemos definir los vectores sobre el plano cartesiano normal de dos dimensiones y, debido a que el álgebra lineal es tan cojonuda, podemos usar los componentes de esos vectores según la base canónica del plano $\{\hat{i},\hat{j}\}$ como nuestros escalares $c_1$ y $c_2$ usados en la definición de los vectores $\vec{v}$ como combinaciones lineales de $\vec{v_1}$ y $\vec{v_2}$. Esto es así porque, después de todo, $\{\vec{v_1},\vec{v_2}\}$ es una base ortonormal del plano $\mathcal{V}$. Así, por ejemplo, si el ángulo entre $(0,1)$ y $(-\tfrac{1}{2},\tfrac{\sqrt{3}}{2})$ es de $120^\circ$, entonces el ángulo entre $\vec{v_2}$ y $-\tfrac{1}{2}\vec{v_1} + \tfrac{\sqrt{3}}{2}\vec{v_2}$ también será de $120^\circ$. No sólo eso, sino que los dos vectores yacerán sobre el plano $\mathcal{V}$, que como es de recordar, es un subespacio de $\mathbb{R}$. De esta forma, podemos traducir los vectores bidimensionales en vectores ubicados en el plano $\mathcal{V}$ perpendicular al tronco $\vec{u}$.

Para escoger 2 or 3 vectors en el plano cartesiano bidimensional para luego convertirlos a vectores de tres coordenadas sobre el plano $\mathcal{V}$ debemos:

2 ramas: Consideramos $(1,0)$ como un primer vector y determinamos el segundo vector de norma unitaria rotando aquel en un ángulo aleatorio entre $160^\circ$ y $200^\circ$. Escogemos otro ángulo entre $0^\circ$ y $360^\circ$ para rotar el conjunto de los dos vecores.
3 ramas: Nuestros tres vectores son inicialmente $(1,0)$, $(\tfrac{\sqrt{3}}{2},-\tfrac{1}{2})$ y $(-\tfrac{\sqrt{3}}{2},-\tfrac{1}{2})$. Escogemos un ángulo entre $0^\circ$ y $360^\circ$ para rotar el conjunto de los tres vectores.

La razón por la cual subsecuentemente rotamos el conjunto de dos o tres vectores es para aleatorizar la orientación de las ramas de modo que las ramas se proyecten en varias direcciones diferentes. En ambos casos, para rotar un vector bidimensional en un ángulo de $\theta$ alrededor del orígen, utilizamos la siguiente transformación lineal:

\[ (x,y)\quad\mapsto (x\,cos(\theta)-y\,sen(\theta),x\,sen(\theta)+y\,cos(\theta)) \]

Hay algunos detalles pendientes por explicar antes de presentar el código en R. Como ya hemos mencionado, los troncos más superficiales del pino de navidad tienen dos ramificaciones mientreas que las partes más profundas, aquellas más cercanas al troncón principal desde donde comienza todo el pino navideo, tienen tres. Esto es controlado por el parámetro de profundidad. Cuando la $profundidad>0.8$, el segmento de tronco se ramifica en tres partes, en caso contrario, se ramifica en dos lugares. El parámetro de profundidad también controla donde a lo largo del segmento de tronco se escogen esos puntos de ramificación. Los segmentos de tronco a mayor profundidad tienen sus puntos de ramificación más hacia el punto terminal $(x_1,y_1,z_1)$ del segmento mientras que los segmentos de tronco más superficiales tienen los puntos de ramificación más cercanos al punto inicial $(x_0,y_0,z_0)$. Huelga decir que estos puntos de ramificación también se escogen aleatoriamente. Finalmete, la colocación de los adornos y las luces también es controlada por el parámetro de profundidad. Las partes más superficiales del árbol tienen más propensión a llevar luces y adornos, como ocurre con los pinos navideños de verdad. Nuevamente, los lugares donde el pino lleva adornos o luces se escogen de manera estocástica.

Script R para el pino navideño tridimensional

library(plot3D)

tronco <- list(
  x0=0, y0=0, z0=0,
  x1=0, y1=0, z1=6,
  extension=NULL,
  rama1=NULL, rama2=NULL,
  rama3=NULL, rama4=NULL,
  profundidad=1,
  grosor=1,
  color_segmento="brown"
)

dibujar_arbol <- function(arbol) {
  #funcion para dibujar el arbol recursivamente
  if (is.null(arbol)) return()
  with(arbol,{
    segments3D(x0,y0,z0,x1,y1,z1,col=color_segmento,
      lwd=grosor,add=TRUE)
    dibujar_arbol(rama1)
    dibujar_arbol(rama2)
    dibujar_arbol(rama3)
    dibujar_arbol(extension)
  })
}

extender <- function(arbol) {
  #crear una extension del tronco principal
  #en la misma direccion que el vector u=p1-p0
  with(arbol,{
    factor_crecimiento <- runif(1,min=0.9,max=1)
    return(list(x0=x1,y0=y1,z0=z1,
         x1=x1+factor_crecimiento*(x1-x0),
         y1=y1+factor_crecimiento*(y1-y0),
         z1=z1+factor_crecimiento*(z1-z0),
         extension=NULL,
         rama1=NULL, rama2= NULL,
         rama3=NULL, rama4= NULL,
         profundidad=profundidad*0.9,
         grosor=1,
         color_segmento="darkgreen"))
  })
}

crear_rama <- function(arbol,c1,c2) {
  #esta funcion retorna ua rama para el arbol
  #c1 y c2 son los componentes para la base
  #vectorial {v1,v2}. v1 y v2 son dos vectores
  #ortonormales, ambos perpediculares al vector
  # direccion u=p1-p0 (p1 y p0 son los puntos
  #inicial y final para el tronco del arbol.
  #adonde es un escalar en (0,1), indicando
  #por donde, entre p0=(x0,y0,z0) y 
  #p1=(x1,y1,z1), crecera la rama.
  with(arbol,{
    #el vector u=(x1-x0,y1-y0,z1-z0) indica la 
    #direccion del tronco del arbol.
    #long_trozo es la longitud de ese troncon.
    u <- c(x1-x0,y1-y0,z1-z0)
    long_trozo <- sqrt(sum(u*u))
    #factor_crecimiento es el multiplo del
    #largo de la rama con respecto a long_tronco
    factor_crecimiento <- runif(1,min=0.7,max=0.8)
    donde <- runif(1,min=0.4,max=0.8)
    #crea dos vectores perpendiculares a u (v1 y v2)
    #de norma unitaria.
    i <- which(u!=0)[1]
    oi <- setdiff(1:3,i)
    v1 <- rep(0,3); v1[i] <- -u[oi[1]]; v1[oi[1]] <- u[i]
    v2 <- rep(0,3); v2[i] <- -u[oi[2]]; v2[oi[2]] <- u[i]
    v1 <- v1/sqrt(sum(v1*v1)); v2 <- v2/sqrt(sum(v2*v2))
    #El vector v es una combinacion lineal de v1 y v2.
    #Por lo tanto, yace sobre el plano perpendicular a u.
    v <- v1*c1+v2*c2
    #beta es el angulo de separacion entre la rama y el
    #tronco. El angulo es mayor para las partes mas
    #"profundas" del arbol.
    beta <- runif(1,min=0.7*profundidad,max=0.9*profundidad)    
    #calcula el escalar k que al ser multiplicado por v
    #en rama=u+kv, esa rama, formara un angulo beta
    #con el tronco.
    k <- tan(beta)*sqrt(sum(u*u))/sqrt(sum(v*v))
    #nuevo_u es el vector direccion de la rama,
    #ya  con la longitud definitiva de la rama.
    nuevo_u <- u+k*v
    nuevo_u <- nuevo_u / sqrt(sum(nuevo_u*nuevo_u))*
             stub_length*factor_crecimiento
    #el nuevo punto (x0,y0,z0) de la rama
    #estara entre p0 y p1 del tronco,
    #en el lugar indicado por "donde".
    nuevo_x0 <- x0+(x1-x0)*donde
    nuevo_y0 <- y0+(y1-y0)*donde
    nuevo_z0 <- z0+(z1-z0)*donde
    #El nuevo (x1,y1,z1) de la rama es su punto
    #terminal.
    nuevo_x1 <- nuevo_x0+nuevo_u[1]
    nuevo_y1 <- nuevo_y0+nuevo_u[2]
    nuevo_z1 <- nuevo_z0+nuevo_u[3]
    list(
      x0=nuevo_x0,y0=nuevo_y0,z0=nuevo_z0,
      x1=nuevo_x1,y1=nuevo_y1,z1=nuevo_z1,
      extension=NULL,
      rama1=NULL, rama2= NULL,
      rama3=NULL, rama4= NULL,
      profundidad=profundidad*0.8,
      grosor=1,
      color_segmento="darkgreen")    
  })
}

crecer_arbol <- function(arbol) {
  if (is.null(arbol) ) return(NULL)
  arbol$grosor <- arbol$grosor*1.4
  if (arbol$grosor>2.0) arbol$color_segmento <- 'brown4'
  if (is.null(arbol$extension)) {
    arbol$extension <- extender(arbol)
    if (arbol$profundidad<0.5) {
      #hacer dos ramas con un agulo de mas de 160 grados
      #entre ellas y rota todo el conjunto
      rot1 <- runif(1,min=0,max=2*pi)
      rot2 <- (runif(1,min=8*pi/9,max=10*pi/9)+rot1)%%(2*pi)
      arbol$rama1 <- crear_rama(arbol,cos(rot1),sin(rot1))
      arbol$rama2 <- crear_rama(arbol,cos(rot2),sin(rot2))
    } else { 
      #hacer 3 ramaes con angulos de 120 grados etre ellas
      #y rota todo el conjunto
      rot <- runif(1,min=0,max=2*pi)
      arbol$rama1 <- crear_rama(arbol,-sin(rot),cos(rot))
      arbol$rama2 <- crear_rama(arbol,
                        sqrt(3)/2*cos(rot)+sin(rot)/2,
                        sqrt(3)/2*sin(rot)-cos(rot)/2)
      arbol$rama3 <- crear_rama(arbol,
                        -sqrt(3)/2*cos(rot)+sin(rot)/2,
                        -sqrt(3)/2*sin(rot)-cos(rot)/2)
    }
  } else {
    arbol$extension <- crecer_arbol(arbol$extension)
    arbol$rama1 <- crecer_arbol(arbol$rama1)
    arbol$rama2 <- crecer_arbol(arbol$rama2)
    if (arbol$profundidad>=0.5) 
      arbol$rama3 <- crecer_arbol(arbol$rama3)
  }
  return(arbol)
}

crear_adornos <- function(arbol) {
  if (is.null(arbol)) return()
  po <- (1-arbol$profundidad)^6
  adorno <- sample(c(T,F),size=1,prob=c(po,1-po))
  co <- sample(c("red","darkgoldenrod4"),size=1,
               prob=c(0.6,0.4))
  if (adorno)
    adornos <<- rbind(adornos,
      data.frame(x=arbol$x1,y=arbol$y1,z=arbol$z1,color=co))
  crear_adornos(arbol$rama1)
  crear_adornos(arbol$rama2)
  crear_adornos(arbol$rama3)
  crear_adornos(arbol$extension)
}

crear_luces1 <- function(arbol) {
  if (is.null(arbol)) return()
  po <- (1-arbol$profundidad)^4
  light <- sample(c(T,F),size=1,prob=c(po,1-po))
  if (light)
    luces1 <<- rbind(luces1,
      data.frame(x=arbol$x1,y=arbol$y1,z=arbol$z1))
  crear_luces1(arbol$rama1)
  crear_luces1(arbol$rama2)
  crear_luces1(arbol$rama3)
  crear_luces1(arbol$extension)
}

crear_luces2 <- function(arbol) {
  if (is.null(arbol)) return()
  po <- (1-arbol$profundidad)^4
  light <- sample(c(T,F),size=1,prob=c(po,1-po))
  if (light)
    luces2 <<- rbind(luces2,
      data.frame(x=(arbol$x1+arbol$x0)/2,
                 y=(arbol$y1+arbol$y0)/2,
                 z=(arbol$z1+arbol$z0)/2))
  crear_luces2(arbol$rama1)
  crear_luces2(arbol$rama2)
  crear_luces2(arbol$rama3)
  crear_luces2(arbol$extension)
}

dibujar_adornos <- function() {
  with(adornos,
    {points3D(x=x,y=y,z=z,
              pch=19,cex=1.2,col=as.character(color),
              colkey=FALSE,add=TRUE)})
}

dibujar_luces1 <- function() {
  with(luces1,
    {points3D(x=x,y=y,z=z,pch="+",cex=0.8,col="white",
              colkey=FALSE,add=TRUE)})
}

dibujar_luces2 <- function() {
  with(luces2,
    {points3D(x=x,y=y,z=z,pch="+",cex=0.8,col="yellow",
              colkey=FALSE,add=TRUE)})
}

set.seed(20161224)
pino_navidad <- tronco
for (i in 1:5) pino_navidad <- crecer_arbol(pino_navidad)
adornos <- NULL;
luces1 <- NULL;
luces2 <- NULL;
crear_adornos(pino_navidad)
crear_luces1(pino_navidad)
crear_luces2(pino_navidad)

png("arbol%02d.png")
for (i in 0:35) {
    perspbox(x=c(-25,25),y=c(-25,25),z=c(0,40),bty="n",
             phi=8,theta=i*10,col="white",alpha=0)
    dibujar_arbol(pino_navidad)
    dibujar_adornos()
    switch((i%%4)+1,{},
     {dibujar_luces1()},
     {dibujar_luces2()},
     {dibujar_luces1(); dibujar_luces2()})
}
graphics.off()

Bibliogrfía

Apostol, T. (1985). Calculus, Vol. II (2nda ed). (F. Veléz trad.). Caracas: Editorial Reverté.
carlop. (April 26, 2012). Respuesta a "How to find vector perpendicular to another vector?" en math.stackexchange. [Recuperado 20/12/2016 desde http://math.stackexchange.com/questions/137362/how-to-find-perpendicular-vector-to-another-vector].
Cooney, B. (1967). Christmas. New York: Thomas Y. Crowell Company.
Duineveld, K. (2013, December). "Merry Christmas". Entrada de blog. [Recuperado 20/12/2016 desde http://wiekvoet.blogspot.com/2013/12/merry-christmas.html].
R DEVELOPMENT CORE TEAM (2009). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0. http://www.R-project.org.
Soetaert, K. (2016). plot3D: Plotting multi-dimensional data. Paquete R. Versión 1.1. https://cran.r-project.org/web/packages/plot3D/index.html.

Estimado lector:

Si te gustó o te pareció útil este contenido, compártelo en las redes sociales y dale tu voto positivo en el botón "me gusta" de G+, para que otros puedan encontrar el contenido también.

Abajo encontrarás un enlace para descargar una versión pdf imprimible de este artículo más el código fuente en R pero esencialmente es la misma información que estoy compartiendo contigo arriba. El enlace de llevará a un gateway de pago en Bitcoin que es muy seguro y anónimo (no recogeré tu información confidencial sensible). El monto es modesto - 0,0005 BTC o aproximadamente $0,45 al cambio de hoy. Considéralo como una donación y no una transacción comercial, que me ayudará a seguir produciendo contenido y a sobrevivir en este hueco infernal en que se ha convertido Venezuela. !Muchas gracias de antemano!