Hlavní navigace

Jednoduché počty s R

15. 11. 1999
Doba čtení: 5 minut

Sdílet

Linux nacházel téměř odjakživa uplatnění na akademické a vědecké půdě. Přesto by se možná mohlo zdát, že aplikací z kategorie vědeckých není zase až tak mnoho. Ale kdo hledá, najde. Třeba právě 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í

UX DAy - tip 2

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

Byl pro vás článek přínosný?

Autor článku

Pavel Stěhule je odborníkem na relační databázový systém PostgreSQL, pracuje jako školitel a konzultant.