Universidade Federal do Paraná
Curso de Estatística
CE 083 -
Estatística Computacional I - 2014/2
Prof. Dr. Walmes Marques Zeviani
Tabela de conteúdo
##-----------------------------------------------------------------------------
## Função para calcular a hipotenusa de um triângulo retângulo.
hipo <- function(a, b){ ## Argumentos.
h <- sqrt(a^2+b^2) ## Processo/procedimento.
return(h) ## Retorno.
}
## Explorando o objeto.
class(hipo)
## [1] "function"
mode(hipo)
## [1] "function"
args(hipo)
## function (a, b)
## NULL
body(hipo)
## {
## h <- sqrt(a^2 + b^2)
## return(h)
## }
str(hipo)
## function (a, b)
## - attr(*, "srcref")=Class 'srcref' atomic [1:8] 4 9 7 1 9 1 4 7
## .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x30b0708>
## Usando a função.
hipo(3, 4)
## [1] 5
hipo(3, 4.6)
## [1] 5.492
hipo(3, 3)
## [1] 4.243
## Note que ela opera vetorialmente.
hipo(a=3, b=c(4, 4.6, 3))
## [1] 5.000 5.492 4.243
##-----------------------------------------------------------------------------
## Abrindo algumas funções.
body(plot) ## Indica que é uma função genérica.
## UseMethod("plot")
methods(generic.function="plot") ## As várias 'faces' da plot.
## [1] plot.acf* plot.data.frame* plot.decomposed.ts* plot.default
## [5] plot.dendrogram* plot.density* plot.ecdf plot.factor*
## [9] plot.formula* plot.function plot.hclust* plot.histogram*
## [13] plot.HoltWinters* plot.isoreg* plot.lm* plot.medpolish*
## [17] plot.mlm* plot.ppr* plot.prcomp* plot.princomp*
## [21] plot.profile.nls* plot.spec* plot.stepfun plot.stl*
## [25] plot.table* plot.ts plot.tskernel* plot.TukeyHSD*
##
## Non-visible functions are asterisked
## Abrindo a função.
plot.default
## function (x, y = NULL, type = "p", xlim = NULL, ylim = NULL,
## log = "", main = NULL, sub = NULL, xlab = NULL, ylab = NULL,
## ann = par("ann"), axes = TRUE, frame.plot = axes, panel.first = NULL,
## panel.last = NULL, asp = NA, ...)
## {
## localAxis <- function(..., col, bg, pch, cex, lty, lwd) Axis(...)
## localBox <- function(..., col, bg, pch, cex, lty, lwd) box(...)
## localWindow <- function(..., col, bg, pch, cex, lty, lwd) plot.window(...)
## localTitle <- function(..., col, bg, pch, cex, lty, lwd) title(...)
## xlabel <- if (!missing(x))
## deparse(substitute(x))
## ylabel <- if (!missing(y))
## deparse(substitute(y))
## xy <- xy.coords(x, y, xlabel, ylabel, log)
## xlab <- if (is.null(xlab))
## xy$xlab
## else xlab
## ylab <- if (is.null(ylab))
## xy$ylab
## else ylab
## xlim <- if (is.null(xlim))
## range(xy$x[is.finite(xy$x)])
## else xlim
## ylim <- if (is.null(ylim))
## range(xy$y[is.finite(xy$y)])
## else ylim
## dev.hold()
## on.exit(dev.flush())
## plot.new()
## localWindow(xlim, ylim, log, asp, ...)
## panel.first
## plot.xy(xy, type, ...)
## panel.last
## if (axes) {
## localAxis(if (is.null(y))
## xy$x
## else x, side = 1, ...)
## localAxis(if (is.null(y))
## x
## else y, side = 2, ...)
## }
## if (frame.plot)
## localBox(...)
## if (ann)
## localTitle(main = main, sub = sub, xlab = xlab, ylab = ylab,
## ...)
## invisible()
## }
## <bytecode: 0x24ea520>
## <environment: namespace:graphics>
plot.ecdf
## function (x, ..., ylab = "Fn(x)", verticals = FALSE, col.01line = "gray70",
## pch = 19)
## {
## plot.stepfun(x, ..., ylab = ylab, verticals = verticals,
## pch = pch)
## abline(h = c(0, 1), col = col.01line, lty = 2)
## }
## <bytecode: 0x2d1af38>
## <environment: namespace:stats>
## plot.ts
## As com asteriscos podem ser abertas assim.
getAnywhere(plot.histogram)
## A single object matching 'plot.histogram' was found
## It was found in the following places
## registered S3 method for plot from namespace graphics
## namespace:graphics
## with value
##
## function (x, freq = equidist, density = NULL, angle = 45, col = NULL,
## border = par("fg"), lty = NULL, main = paste("Histogram of",
## paste(x$xname, collapse = "\n")), sub = NULL, xlab = x$xname,
## ylab, xlim = range(x$breaks), ylim = NULL, axes = TRUE, labels = FALSE,
## add = FALSE, ann = TRUE, ...)
## {
## equidist <- if (is.logical(x$equidist))
## x$equidist
## else {
## h <- diff(x$breaks)
## diff(range(h)) < 1e-07 * mean(h)
## }
## if (freq && !equidist)
## warning("the AREAS in the plot are wrong -- rather use 'freq = FALSE'")
## y <- if (freq)
## x$counts
## else x$density
## nB <- length(x$breaks)
## if (is.null(y) || 0L == nB)
## stop("'x' is wrongly structured")
## dev.hold()
## on.exit(dev.flush())
## if (!add) {
## if (is.null(ylim))
## ylim <- range(y, 0)
## if (missing(ylab))
## ylab <- if (!freq)
## "Density"
## else "Frequency"
## plot.new()
## plot.window(xlim, ylim, "")
## if (ann)
## title(main = main, sub = sub, xlab = xlab, ylab = ylab,
## ...)
## if (axes) {
## axis(1, ...)
## axis(2, ...)
## }
## }
## rect(x$breaks[-nB], 0, x$breaks[-1L], y, col = col, border = border,
## angle = angle, density = density, lty = lty)
## if ((logl <- is.logical(labels) && labels) || is.character(labels))
## text(x$mids, y, labels = if (logl) {
## if (freq)
## x$counts
## else round(x$density, 3)
## }
## else labels, adj = c(0.5, -0.5))
## invisible()
## }
## <bytecode: 0x3439f58>
## <environment: namespace:graphics>
## getS3method(f="plot", class="histogram")
## Abrindo a função sequência.
methods(generic.function="seq")
## [1] seq.Date seq.default seq.POSIXt
getAnywhere(seq.default)
## A single object matching 'seq.default' was found
## It was found in the following places
## package:base
## registered S3 method for seq from namespace base
## namespace:base
## with value
##
## function (from = 1, to = 1, by = ((to - from)/(length.out - 1)),
## length.out = NULL, along.with = NULL, ...)
## {
## if ((One <- nargs() == 1L) && !missing(from)) {
## lf <- length(from)
## return(if (mode(from) == "numeric" && lf == 1L) {
## if (!is.finite(from)) stop("'from' cannot be NA, NaN or infinite")
## 1L:from
## } else if (lf) 1L:lf else integer())
## }
## if (!missing(along.with)) {
## length.out <- length(along.with)
## if (One)
## return(if (length.out) seq_len(length.out) else integer())
## }
## else if (!missing(length.out)) {
## len <- length(length.out)
## if (!len)
## stop("argument 'length.out' must be of length 1")
## if (len > 1L) {
## warning("first element used of 'length.out' argument")
## length.out <- length.out[1L]
## }
## length.out <- ceiling(length.out)
## }
## if (!missing(...))
## warning(sprintf(ngettext(length(list(...)), "extra argument %s will be disregarded",
## "extra arguments %s will be disregarded"), paste(sQuote(names(list(...))),
## collapse = ", ")), domain = NA)
## if (!missing(from) && length(from) != 1L)
## stop("'from' must be of length 1")
## if (!missing(to) && length(to) != 1L)
## stop("'to' must be of length 1")
## if (!missing(from) && !is.finite(from))
## stop("'from' cannot be NA, NaN or infinite")
## if (!missing(to) && !is.finite(to))
## stop("'to' cannot be NA, NaN or infinite")
## if (is.null(length.out))
## if (missing(by))
## from:to
## else {
## del <- to - from
## if (del == 0 && to == 0)
## return(to)
## n <- del/by
## if (!(length(n) && is.finite(n))) {
## if (length(by) && by == 0 && length(del) && del ==
## 0)
## return(from)
## stop("invalid (to - from)/by in seq(.)")
## }
## if (n < 0L)
## stop("wrong sign in 'by' argument")
## if (n > .Machine$integer.max)
## stop("'by' argument is much too small")
## dd <- abs(del)/max(abs(to), abs(from))
## if (dd < 100 * .Machine$double.eps)
## return(from)
## if (is.integer(del) && is.integer(by)) {
## n <- as.integer(n)
## from + (0L:n) * by
## }
## else {
## n <- as.integer(n + 1e-10)
## x <- from + (0L:n) * by
## if (by > 0)
## pmin(x, to)
## else pmax(x, to)
## }
## }
## else if (!is.finite(length.out) || length.out < 0L)
## stop("length must be non-negative number")
## else if (length.out == 0L)
## integer()
## else if (One)
## seq_len(length.out)
## else if (missing(by)) {
## if (missing(to))
## to <- from + length.out - 1L
## if (missing(from))
## from <- to - length.out + 1L
## if (length.out > 2L)
## if (from == to)
## rep.int(from, length.out)
## else as.vector(c(from, from + seq_len(length.out -
## 2L) * by, to))
## else as.vector(c(from, to))[seq_len(length.out)]
## }
## else if (missing(to))
## from + (0L:(length.out - 1L)) * by
## else if (missing(from))
## to - ((length.out - 1L):0L) * by
## else stop("too many arguments")
## }
## <bytecode: 0x375d7a0>
## <environment: namespace:base>
## Abrindo a matrix.
matrix
## function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)
## {
## if (is.object(data) || !is.atomic(data))
## data <- as.vector(data)
## .Internal(matrix(data, nrow, ncol, byrow, dimnames, missing(nrow),
## missing(ncol)))
## }
## <bytecode: 0x2485098>
## <environment: namespace:base>
runif ## Chama uma função escrita em C.
## function (n, min = 0, max = 1)
## .External(C_runif, n, min, max)
## <bytecode: 0x2e24420>
## <environment: namespace:stats>
##-----------------------------------------------------------------------------
## Incluindo o argumento ângulo.
hipo2 <- function(a, b, angle){
h2 <- a^2+b^2-2*cos(angle)
return(sqrt(h2))
}
hipo2(3, 4, angle=pi/2)
## [1] 5
hipo(3, 4)
## [1] 5
hipo2(3, 4, angle=pi/2+0.3)
## [1] 5.059
hipo2(3, 4, angle=pi/2-0.3)
## [1] 4.941
##-----------------------------------------------------------------------------
## Ângulo com valor default de 90 graus (pi/2 radianos).
hipo3 <- function(a, b, angle=pi/2){
h2 <- a^2+b^2-2*cos(angle)
return(sqrt(h2))
}
hipo3(3, 4)
## [1] 5
hipo3(3, 4, angle=pi/2+0.5)
## [1] 5.095
##-----------------------------------------------------------------------------
## Três implementações.
baskara1 <- function(a, b, c){
x1 <- (-b-sqrt(b^2-4*a*c))/(2*a)
x2 <- (-b+sqrt(b^2-4*a*c))/(2*a)
return(c(x1, x2))
}
baskara2 <- function(a, b, c){
s <- sqrt(b^2-4*a*c)
d <- 2*a
x1 <- (-b-s)/d
x2 <- (-b+s)/d
return(c(x1, x2))
}
baskara3 <- function(a, b, c){
x <- (-b+c(-1,1)*sqrt(b^2-4*a*c))/(2*a)
return(x)
}
a <- 2; b <- 1; c <- -3
curve(a*x^2+b*x+c, -3, 3)
abline(h=0, lty=2)
baskara1(a, b, c)
## [1] -1.5 1.0
baskara2(a, b, c)
## [1] -1.5 1.0
baskara3(a, b, c)
## [1] -1.5 1.0
##-----------------------------------------------------------------------------
## Tempo para 50 mil excussões de cada uma.
system.time(replicate(50000, baskara1(a, b, c)))
## user system elapsed
## 0.472 0.004 0.477
system.time(replicate(50000, baskara2(a, b, c)))
## user system elapsed
## 0.365 0.000 0.366
system.time(replicate(50000, baskara3(a, b, c)))
## user system elapsed
## 0.333 0.000 0.333
## As implementação diferentes implicam em custos diferentes.
##-----------------------------------------------------------------------------
## Parábola sem raízes.
a <- 2; b <- 1; c <- 3
curve(a*x^2+b*x+c, -5, 5)
abline(h=0, lty=2)
baskara(a, b, c)
## Error: could not find function "baskara"
##-----------------------------------------------------------------------------
## a==0, então não tem curvatura.
a <- 0; b <- 1; c <- 3
curve(a*x^2+b*x+c, -5, 5)
baskara(a, b, c)
## Error: could not find function "baskara"
## Como notificar da origem desses erros?
##-----------------------------------------------------------------------------
## Colocando mensagens de erro para testes feitos nos argumentos.
baskara4 <- function(a, b, c){
if (a==0){
stop("O argumento `a` tem que ser diferente de zero.")
}
delta <- b^2-4*a*c
if (delta<=0) {
stop("Função sem raíz pois `delta` é não positivo.")
}
x <- (-b+c(-1,1)*sqrt(delta))/(2*a)
return(x)
}
baskara4(2, 1, -3)
## [1] -1.5 1.0
baskara4(0, 1, -3)
## Error: O argumento `a` tem que ser diferente de zero.
baskara4(2, 1, 3)
## Error: Função sem raíz pois `delta` é não positivo.
##-----------------------------------------------------------------------------
##-----------------------------------------------------------------------------
## Função que calcula as raízes da equação e ainda os valores no ponto
## crítico (x, y). O ponto crítico sempre existe se a!=0. As raízes
## podem não existir.
baskara5 <- function(a, b, c){
if (a==0){
stop("O argumento `a` tem que ser diferente de zero.")
}
d <- 2*a
delta <- b^2-4*a*c
if (delta<=0){
warning("Função sem raíz pois `delta` é não positivo.")
r <- c(NA, NA)
} else {
r <- (-b+c(-1,1)*sqrt(delta))/d
}
x <- -b/d
y <- a*x^2+b*x+c
L <- list(x=x, y=y, r=r)
return(L)
}
a <- baskara5(2, 1, -3)
str(a)
## List of 3
## $ x: num -0.25
## $ y: num -3.12
## $ r: num [1:2] -1.5 1
curve(2*x^2+1*x-3, -2, 2)
abline(v=a$x, h=a$y, col=2, lty=2)
abline(v=a$r, h=0, col=3, lty=2)
##-----------------------------------------------------------------------------
## Argumentos como vetores, usa posição.
baskara6 <- function(abc){
if (length(abc)==3L){
x <- (-abc[2]+c(-1,1)*sqrt(abc[2]^2-4*abc[1]*abc[3]))/(2*abc[1])
return(x)
} else {
stop("O vetor `abc` deve ter comprimento 3.")
}
}
baskara6(c(2,1,-3))
## [1] -1.5 1.0
baskara6(c(2,1,-3,500))
## Error: O vetor `abc` deve ter comprimento 3.
baskara6(c(2,1))
## Error: O vetor `abc` deve ter comprimento 3.
##-----------------------------------------------------------------------------
## Argumentos como vetores nomeados.
baskara7 <- function(abc){
if (length(abc)==3L & all(names(abc)%in%c("a","b","c"))){
x <- (-abc["b"]+c(-1,1)*
sqrt(abc["b"]^2-4*abc["a"]*abc["c"]))/(2*abc["a"])
return(x)
} else {
stop("O vetor `abc` deve ter comprimento 3 e ser nomeado.")
}
}
baskara7(c(a=2, b=1, c=-3))
## [1] -1.5 1.0
baskara7(c(b=1, a=2, c=-3))
## [1] -1.5 1.0
baskara7(c(b=1, a=2, m=-3))
## Error: O vetor `abc` deve ter comprimento 3 e ser nomeado.
##-----------------------------------------------------------------------------
## Argumentos como vetores ou lista nomeados.
baskara8 <- function(abc){
if (is.vector(abc)) abc <- as.list(abc)
if (is.list(abc) & length(abc)==3L & all(names(abc)%in%c("a","b","c"))){
x <- with(abc,
(-b+c(-1,1)*sqrt(b^2-4*a*c))/(2*a))
return(x)
} else {
stop("O objeto `abc` deve ser vetor/lista nomeado de 3 elementos.")
}
}
baskara8(c(a=2, b=1, c=-3))
## [1] -1.5 1.0
baskara8(list(a=2, b=1, c=-3))
## [1] -1.5 1.0
baskara8(data.frame(a=2, b=1, c=-3))
## [1] -1.5 1.0
