Jednoduché počty s R

Nejen programováním živ je člověk, dalo by se parafrázovat oblíbené přísloví. Někdy je skutečně nutné dělat něco jiného. Pokud by se to jiné odehrávalo na poli matematiky, a to zejména lineární algebry nebo statistiky, a chtěli jste si ušetřit manuální práci bez toho, že by jste si přidělali práci s instalováním Excelu nebo MathLabu, a nezatěžovali svůj účet nebo svědomí, podívejte se po prostředí R (šířeno pod GNU licencí). V podstatě se jedná jen o interpret jednoduchého jazyka, který se stejně učit nemusíte, a o sadu knihoven funkcí většinou statisticky orientovaných, které také používat nemusíte (Nejedná se jen o nějakou sadu knihoven. Je velmi pravděpodobné, že při troše vůle a hledání v nich najdete řešení jakékoliv statistické úlohy, která vás napadne.). Kromě toho R ještě disponuje celkem slušnou schopností grafického výstupu (X11, WIN nebo LaTeX, postscript).

Prostředí, tj. interpret, si můžete nainstalovat z binárek nebo přeložit. K tomu budete potřebovat kromě překladače C, např. gcc, ještě překladač Fortranu. Mimochodem bylo to poprvé, co jsem viděl Fortran v akci. Existují verze jak pro Linux a ostatní Unixy, tak pro WIN32 a Mac. A pokud si rozumíte s Emacsem, tak si nezapomeňte nainstalovat podporu SCC

Je nediskutovatelnou pravdou, že pravý programátor něco takového jako R nepotřebuje. Jelikož je R psané převážně v C, tak vše co dokáže v Rku, v C dokáže také. Na druhou stranu, pokud nemá perfektní třídy pro matice, nebo briskně neovládá STL, tak napsat kód pro násobení matic, nemluvě o něčem složitějším, nějakou tu půlhodinku zabere. O to větší důvod podívat se po něčem jiném. Nezávislí experti doporučují R. Zkuste si:

> x<-c(1,2,3,4)
> y<-c(4,5,6,7)
> x%*%y
     [,1]
[1,]   60

V příkladu se počítá skalární součin dvou čtyř-prvkových vektorů. Hodnoty vektorů jsou uložené v proměnných x a y. Pro skalární součet se v R používá symbol %*%, obvyklá * je rezervována pro složkový součin. Jak jednoduché. Symbol ← se používá ve smyslu přiřazovacího příkazu. Funkce matrix použijeme jako konstruktor (funkce vytvářející určitý objekt) matice. Jejím prvním parametrem je vektor všech prvků v matici (nemusí se jednat pouze o čísla). V matici tyto prvky se poté uspořádají po sloupcích. Další argumenty mají význam počet řádků a počet sloupců. Vyzkoušejte si:

> sin(1+diag(3))
          [,1]      [,2]      [,3]
[1,] 0.9092974 0.8414710 0.8414710
[2,] 0.8414710 0.9092974 0.8414710
[3,] 0.8414710 0.8414710 0.9092974

Zkuste se zamyslet, jak dlouho byste psali program pro výpočet sinu každého prvku jednotkové matice o třech řádcích a třech sloupcích zvětšeného o hodnotu jedna. V jakémkoliv jiném jazyce by byl tento výraz nesmyslný (R není určeno pro kritické aplikace a tak není třeba bazírovat na shodě typů). Funkce diag(n) vrací jednotkovou matici o n řádcích a n sloupců, a této matici přičítám konstantu. Parametrem funkce sin je pouze hodnota, nikoliv matice atd. R si „domyslí“ jak jsme to mysleli, zde vytvoří matici n,n prvků, kdy každý prvek matice obsahuje hodnotu jedna, s kterou již může provést složkový součet s jednotkovou maticí. A jelikož funkce sin neakceptuje matici, tak provede operaci sin s každým prvkem matice. Zatímco zmíněný výraz nemá smysl, následující

sum((x-mean(x))^2)/(length(x) - 1)

má. Vyjadřuje výpočet rozptylu na souboru hodnot, uložených ve vektoru x. Hezké na tom je to, že se jedná, v podstatě, o pouhý přepis vzorečku. Jen pro úplnost dodám, že v R naleznete funkce pro výpočet inverzních matic, determinantů a pro řešení lineárních rovnic. Zatím jsem tady R prezentoval jen jako chytřejší kalkulačku, R je však také plnohodnotným programovacím jazykem. Důkazem nechť je následující program řešící úlohu lineárního programování (jde o nalezení maxima funkce funkce na množině S dané omezujícími podmínkami x>=0; Ax<=b):

## Řešení úlohy LP simplexovou metodou, tj. úlohy najít
## maximum funkce c.x, přičemž x >= 0 a A.x <= b. Předpokládáme,
## že úloha LP je nedegenerovaná. Postup řešení: Vybrané statě
## z operačního výzkumu, RNDr. Ing. Jaroslav Kvaňa, 1987

simplex <- function (c, A, x, b, debug.print=F)
  {
    m <- length (b)
    n <- ncol(A)

    if (m != nrow (A))
      stop ("Pocet radku matice A je jiny nez delka vektoru b")

    if (length (c) != n)
      stop ("Pocet sloupcu matice A je jiny nez delka vektoru c")

    A  <- matrix (c (A, diag (m)), m, m+n)
    x  <- c (rep (0,n), b)
    c  <- c (c, rep (0, m))
    cn <- rep (0, m)
    ai <- (n+1):(m+n)
    d  <- NA
    iterace <- 0

    repeat {
      iterace <- iterace + 1
      for (j in 1:(n+m))
        d[j] <- cn %*% A[,j] - c[j]

      M <- min (d)
      if (M >= 0) {
        x <<- x[1:n]
        return (as.vector(cn%*%b))
      }

      for (j in 1:(m+n))
        if (d[j] == M)
          break

      bj <- b / A [,j]
      amin <- NA

      for (ip in 1:m)
        if (!((bj[ip] == Inf)|(bj[ip] <= 0)))
          if (is.na(amin)) {
            amin <- bj[ip]
            i <- ip
          }
          else
            if (amin > bj[ip]) {
              i <- ip
              amin <- bj[ip]
            }

      if (is.na(amin))
          stop ("Nebyl nalezen ridici prvek")

      for (k in 1:m)
        if (k != i) {
          koef  <- - A[k,j] / A[i,j]
          A[k,] <- A[i,]*koef + A[k,]
          b[k]  <- b[i]*koef + b[k]
        } else {
          b[k]  <- b[k]  / A[k,j]
          A[k,] <- A[k,] / A[k,j]
        }

      ai[i] <- j
      cn[i] <- c[j]
      x[ai] <- b

      if (debug.print) {
        Ta <- matrix(NA,m+1,m+n+3)
        Ta[1:m,3:(m+n+3)]<- matrix(c(A,b),m,m+n+1)
        Ta[1:m,1] <- cn
        Ta[1:m,2] <- ai
        Ta[m+1,3:(m+n+3)] <- c(d,cn%*%b)

        cat("
krok",iterace,",x =(",x,"), [",i,j,"]
")
        print(Ta)
      }
    }
  }

## Vyzkoušejte
## A <- matrix (c(5,6,-1,2,12,1),3,2)
## b <- c(20,72,3)
## c <- c(16,6)
##
## pak L <- simplex (c,A,x,b)
## L ==> 64 a x ==> (4,0)

Jak vidno, R může být velice užitečným pomocníkem tam, kde je třeba snadno a rychle realizovat složitější výpočty a matematické operace, stejně jako je třeba Perl zažitým pomocníkem pro zpracování textových reportů.