Hlavní navigace

Jednoduché počty s R

Pavel Stěhule

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í

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

Našli jste v článku chybu?

19. 1. 2009 11:44

Rko je bezvadný, sám ho často používám, nicméně, několik připomínek:
Rko je freewarová verze S+ (jazyk je S+).
Při exportu grafů někdy docela skřípu zubama. Nabízí to pestrou paletu grDevices, nicméně exportovat graf na A4 v 1000dpi (požaduje elzevír) to je vopravdu na mašli, export do eps mi zmršil legendy grafů. Export do PDF je moc pěknej.
Co dělá Rko zajímavym je především miliarda balíčků (od gis, přes statistiky, účetnictví, připojeni na MySQL...). To ušetří dost práce.
Taky se mi ho…


120na80.cz: 5 nejčastějších mýtů o kondomech

5 nejčastějších mýtů o kondomech

Podnikatel.cz: Změny v cestovních náhradách 2017

Změny v cestovních náhradách 2017

Podnikatel.cz: K EET. Štamgast už peníze na stole nenechá

K EET. Štamgast už peníze na stole nenechá

Vitalia.cz: Taky věříte na pravidlo 5 sekund?

Taky věříte na pravidlo 5 sekund?

DigiZone.cz: TV Philips a Android verze 6.0

TV Philips a Android verze 6.0

Root.cz: Certifikáty zadarmo jsou horší než za peníze?

Certifikáty zadarmo jsou horší než za peníze?

Vitalia.cz: Láska na vozíku: Přitažliví jsme pro tzv. pečovatelky

Láska na vozíku: Přitažliví jsme pro tzv. pečovatelky

Podnikatel.cz: Udávání kvůli EET začalo

Udávání kvůli EET začalo

Vitalia.cz: Proč vás každý zubař posílá na dentální hygienu

Proč vás každý zubař posílá na dentální hygienu

Podnikatel.cz: Chaos u EET pokračuje. Jsou tu další návrhy

Chaos u EET pokračuje. Jsou tu další návrhy

Vitalia.cz: Spor o mortadelu: podle Lidlu falšovaná nebyla

Spor o mortadelu: podle Lidlu falšovaná nebyla

Vitalia.cz: Paštiky plné masa ho zatím neuživí

Paštiky plné masa ho zatím neuživí

Vitalia.cz: Jmenuje se Janina a žije bez cukru

Jmenuje se Janina a žije bez cukru

Lupa.cz: Proč firmy málo chrání data? Chovají se logicky

Proč firmy málo chrání data? Chovají se logicky

120na80.cz: Pánové, pečujte o svoje přirození a prostatu

Pánové, pečujte o svoje přirození a prostatu

Lupa.cz: Co se dá měřit přes Internet věcí

Co se dá měřit přes Internet věcí

Lupa.cz: Teletext je „internetem hipsterů“

Teletext je „internetem hipsterů“

DigiZone.cz: Rádio Šlágr má licenci pro digi vysílání

Rádio Šlágr má licenci pro digi vysílání

Vitalia.cz: Jsou čajové sáčky toxické?

Jsou čajové sáčky toxické?

Měšec.cz: mBank cenzuruje, zrušila mFórum

mBank cenzuruje, zrušila mFórum