Objetivo:

“The most damaging phrase in any language is: ‘IT’s ALWAYS BEEN DONE THAT WAY’” (Grace Hopper)

Relação de Recorrência:

Em Matemática, uma relação de recorrência é uma equação em que o n-ésimo termo de uma sequência de números é igual a algum tipo de combinação dos termos anteriores. Frequentemente, apenas k termos anteriores da sequência compõe a equação. Nas recorrências lineares, n-ésimo termo é função linear dos k termos anteriores, onde k é denominada a ordem da relação. Um exemplo é a recorrência para os números de Fibonacci:

\[ F_n = F_{n-1} + F_{n-2} . \]

Esse é um exemplo de uma recorrência linear com coeficientes constantes, os quais não dependem de n. Para essas recorrências, pode-se expressar o termo geral como uma expressão fechada de n. Resolver uma relação de recorrência significa obter uma solução fechada: um função não recursiva de n.

As relações de recorrência são compostas pelas condições iniciais, que devem ser conhecidas, e pela equação de recorrência.

Função Recursiva:

Uma função é dita recursiva quando ela chama a si mesma uma ou mais vezes. A ideia é simples e busca-se dividir um problema original um subproblemas menores de mesma natureza (divisão) e depois combinar as soluções obtidas para gerar a solução do problema original de tamanho maior (conquista). Os subproblemas são resolvidos recursivamente do mesmo modo em função de instâncias menores, até se tornarem problemas triviais que são resolvidos de forma direta, interrompendo a recursão. Em resumo, para resolver um problema X usando uma função f() procede-se da seguinte maneira:

  1. Divida o problema do tipo X original em um ou mais problemas do tipo X.

  2. Com f(), chama-se f() em cada um dos subproblemas.

  3. Com f(), juntam-se os resultados obtidos em (2) para solucionar o problema original.

Seja por exemplo a soma dos \(n\) primeiros números naturais (\(S_n\)). A equação de recorrência é \(S_n = n + S_{n-1}\) e a condição de parada é \(S_1 = 1\). Assim, o problema é dividido em subproblemas:

  1. \(S_5 = 5 + S_4\)

  2. \(S_4 = 4 + S_3\)

  3. \(S_3 = 3 + S_2\)

  4. \(S_2 = 2 + S_1\)

  5. \(S_1 = 1\) (condição de parada)

Essa técnica usa o conceito de recursão para executar tarefas iterativas (passos a. a b.) que chama a si mesma recursivamente, dividindo o problema em componentes menores e chamando a função em cada um desses componentes. Os resultados de cada componente são usados para resolver o problema original. Esse tipo de função exige uma condição de parada para evitar um looping contínuo (passo final e.). O código abaixo implementa a ideia detalhada anteriormente:

  soma <- function(n){
        if (n == 1){
            return(1)
        }else{
          return(n + soma(n - 1))
        }
}

Para calcular a soma dos cinco primeiros números naturais (\(S_5\)) usamos o comando soma(5) que nos fornece o valor de 15. O exemplo acima foi utilizado para exemplificar o processo de recursão, pois para calcular \(S_5\) deveríamos usar o comando sum(1:5), aproveitando-se das propriedades de vetorização do R.

Seguem dois outros exemplos que, por sua simplicidade, ajudam a entender o uso de recursão na construção de funções, embora essa não seja a abordagem mais adequada para obtenção da solução:

Para calcular a potência \(B^n\), com \(n \in \mathbb{N}\), a equação de recursão é \(f(n) = B f(n-1)\) (\(B^n = B \times B^{n-1}\))e a condição de parada é \(f(0) = 1\) (\(B^0 = 1\)). Verifique o código para sua implementação e detalhe passo a passo o comportamento da recursão para calcular, por exemplo \(4^4\):

potencia <- function(base, n) {
    if (n == 0) return(1)
    return(base * potencia(base, n-1))
}

O comando potencia(4, 4) fornece o valor 256, que pode ser obtido com mais eficiência pelo comando do R Base 4**4 (está implementado em C++). Verifique em detalhes também o código abaixo, que utiliza recursão para calcular a soma \(\sum_{i = 1}^n x_i^2\) do vetor vet = c(x1, x2, ..., xn). Verifique o código para sua implementação e detalhe-o passo a passo para calcular, por exemplo, a soma do quadrado dos dez primeiros números naturais.

soma.serie <- function(vet){
    if(length(vet)<=1) return(vet^2)
     return(vet[1]^2 + soma.serie(vet[-1]))
}

Assim, executamos o comando soma.serie(series), em que series <- c(1:10) e seu resultado é 385. Claro que para calcular essa soma, é mais adequado usarmos o comando sum(series^2). PAREI AQUI!!!!!

A recursão é uma maneira elegante para resolver muitos problemas, sendo uma abordagem utilizada quando o loop aumenta a necessidade de alocação de memória.

Exemplos

Os exemplos abaixo ilustram o que significa uma função chamar a si mesma.

Fatorial:

O fatorial é definido pela reção de recorrência

\[ n! = n (n-1)! \text{, para } \forall n>0, \, n \in \mathbb{N} \text{, com a condição inicial de } 0! = 1 . \]

Essa é uma recorrência linear com coeficientes polinomiais de ordem 1 do tipo \(f(n) = n f(n-1)\). O código baixo implementa função recursiva para o cálculo do fatorial do número natural \(n\):

            # ------- Fatorial -------


fatorial <- function(n){
    if(n == 0) {
        1
    } else {
        n * fatorial(n-1)
    }
}

Assim, calula-se facilmente o fatorial de 10, por exemplo:

fatorial(5)
## [1] 120

Essa função é construída como exemplo, pois, em R, podemos calcular esse fatorial de maneira vetorizada, usando, por exemplo, um produtório:

 prod(1:5)
## [1] 120

Executamos o comando fatorial(5) para encontrar o valor de \(5!\), a função retorna um valor como 5*fatorial(5-1), ou seja, dentro da função fatorial chama-se novamente a mesma função com o inteiro anterior, \(4\). Esse procedimento continuará até a função fatorial(0) retornar o valor de (condição de parada). Finalmente, o produto de cada valor levará à resposta 120.

Fonte: Blog Learn R programming.
Fonte: Blog Learn R programming.

Números de Fibonacci:

Os números de Fibonacci formam uma recorrência de ordem 2. Eles são um exemplo canônico de uma recorrência linear homogênea, com coeficientes constantes. A sequência de Fibonacci é definida usando a equação de recorrência \(F_n = F_{n-1} + F_{n-2}\), com condições iniciais \(F_0 = 0\) e \(F_1 = 1\) (ou \(F_1 = 1\) e \(F_2 = 1\)). Os termos iniciais da sequência de Fibonacci são: \(0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, ...\). É frequente a ocorrência dessa interessante sucessão de Fibonacci na natureza. Esse clipe é especial e muito interessante.

Com o código da função fibo() abaixo, podemos determinar o (n+1)-ésimo termo inicial da sequência de Fibonacci (considerando que o primeiro termo é zero):

fibo <- function(n) {
         # Calcula f(n)
        if(n == 0 || n == 1) {
            return(n)
        } else {
            return(fibo(n - 1) + fibo(n - 2))
        }
}

Note Que poderíamos implementar o código de maneira ligeiramente diferente. Quando a condição inicial pro verdadeira (if(n == 0 || n == 1)), então o comando return(n) retornará o valor de n e encerrará a execução da função. Assim, na realidade, não necessitamos da condição else, pois as linhas finais sempre serão executadas enquanto a recorrência não chegar em 1 e 0 (último passo da recorrência). Dessa maneira, podemos utilizar o código abaixo, que está mais “limpo”:

fibo <- function(n) {
      # Calcula f(n)
     if(n == 0 || n == 1) return(n) # condição de parada
       # equação de recorrência
     return(fibo(n - 1) + fibo(n - 2))
}

O comando fibo(10) fornece o valor do 11º elemento (considerando ser zero o 1º elemento) da sequência de Fibonacci que é 55. A estrutura da função é bem simples. A função tem uma única condicional (if(n == 0 || n == 1)) que define a condição de parada do código. Enquanto o resultado essa condição for FALSE a função retorna a soma fibo(n - 1) + fibo(n - 2). Perceba que essa recursão continua até chegar nos dois primeiros termos, os quais acionam a condição de parada. Abaixo um diagrama para visualizar as chamadas recursivas da função, ao executar o comando fibo(5).

Árvore de Recorrência (Fonte: DAVIES, 2016.)
Árvore de Recorrência (Fonte: DAVIES, 2016.)

É importante notar que a função fibo() irá funcionar desde que o usuário forneça um número inteiro não negativo para o argumento n. Se n for negativo, a condição de parada nunca será atendida e a recursão será repetida indefinidadmente (loop infinito). Salienta-se que em casos mais simples, como nesse exemplo de números de Fibonacci, a abordagem por recursão frequentemente exige um custo computacional maior do que uma abordagem por looping iterativo ou usando dos recursos de vetorização oferecidos pelo R.

Conversão Decimal em Binário:

Para registrar uma quantidade inteira de objetos, seres ou coisas, somos ensinados a recorrer ao uso dos símbolos 0, 1, 2, 3, 4, 5, 6, 7, 8 e 9, chamados algarismos, que agrupados segundo poucas regras, formam um numeral do Sistema de Numeração Decimal. Representações de quantidades por meio de registros escritos formados por algarismos são chamadas de numerais. Um número (natural) é a expressão abstrata de uma quantidade. O termo algarismo refere-se a cada um dos símbolos que são combinados para representar números em um dado sistema de numeração, e cada uma dessas combinações de algarismos é chamada um numeral. Um mesmo número admite diferentes representações por meio de diferentes numerais. Dessa maneira conceitua-se Sistema de Numeração como todo conjunto de regras para produção sistemática de numerais.

Os sistemas babilônio, maia e indo-arábico de bases 60, 20 e 10, respectivamente, mantêm em comum a característica de que o valor dos seus símbolos (ou algarismos) depende da posição no agrupamento. A natureza posicional, bem como a certeza sobre como se organiza cada numeral de um sistema de numeração posicional, permitem a busca pela representação de um número em diferentes bases. Vejamos por exemplo o número expresso nas bases 2, 5, 7 e 9:

  1. Sistema decimal: \(2023 = 2\times 10^3 + 0 \times 10^2 + 2 \times 10^1 + 3 \times 10^0\)

  2. Sistema binário: \(2023 = 1 \times 2^{10} + 1 \times 2^{9} + 1 \times 2^{8} + 1 \times 2^{7} + 1 \times 2^{6} + 1 \times 2^{5} + 0 \times 2^{4} + 0 \times 2^{3} + 1 \times 2^{2} + 1 \times 2^{1} + 1 \times 2^{0}\) = \([11111100111]_2\).

  3. Sistema na base 5: \(2023 = 3 \times 5^{4} + 1 \times 5^{3} + 0 \times 5^{2} + 4 \times 5^{1} + 3 \times 5^{0} = [31043]_5\)

  4. Sistema na base 7: \(2023 = 5 \times 7^{3} + 6 \times 7^{2} + 2 \times 7^{1} + 0 \times 7^{0} = [5620]_7\)

  5. Sistema na base 9: \(2023 = 2 \times 9^{3} + 6 \times 9^{2} + 8 \times 9^{1} + 7 \times 9^{0} = [2687]_9\)

O sistema binário é a base da tecnologia digital, qualquer dispositivo que tenha circuitos integrados (chips) só é possível graças a esse sistema numérico. Por conta do sistema binário ser o sistema numérico mais simples, porque só usa dois dígitos, é possível armazenar e manusear os números em forma física, por exemplo, como a eletricidade que viaja através de um fio ou cabo, também como a luz que viaja através de uma fibra ótica, como campos magnéticos em um disco rígido ou ainda como sinais eletromagnéticos que viajam através do ar. Voce pode conhecer um pouco mais sobre os sistema binário e seus usos nos posts Como funciona o sistema binário?;

Para converter um número decimal no seu equivalente binário, devemos fazer sucessivas divisões sobre a base, que é 2, até que não seja mais possível dividir. Em seguida, construímos o número binário com o último resultado (valor menor que a base) seguido do resto das divisões anteriores. Por exemplo, o decimal \(25\) é \([11001]_2\), pois:

  • $25/2 = 12 , com resto \(1\).

  • \(12/2 = 6\), com resto \(0\).

  • \(6/2 = 3\), com resto \(0\)

  • \(3/2 = 1\), com resto \(1\)

  • 0 binário será composto pelo resultado final (\(1\)) seguido dos restos anteriores em ordem inversa (do último resto obtido para o primeiro), ou seja \(11001\). Esse procedimento de transformação esta detalhado no vídeo a seguir (bem curto) e no post Os números binários

detalhadamente aqui neste vídeo.

Assim, um procedimento recursivo para converter um número no sistema decimal em binário está codificado abaixo:

dec2bin <- function(n){
 bin <- ''
  if(n > 1) bin <- dec2bin(as.integer(n/2)) # a condição de parada é n<=1
  bin <- paste0(bin, n %% 2)
  return(as.numeric(bin))
}

O binário de \(25\) é então \(11001\). A condição de parada (if(n > 1) bin <- dec2bin(as.integer(n/2))) retorna o último resultado. Enquanto o algoritmo trabalha recursivamente, os restos estão sendo concatenando na string bin que é construída ao longo de todo o procedimento, iniciando com o último resultado concatenando pos restos anteriores em ordem inversa. O algoritmo é bem simples, mas é necessário que verifique cada um de seus passos, para \(n = 8\), por exemplo. é executada internamente a função dec2bin(8/2), que chama por sua vez a função dec2bin(4/2), que chama a função dec2bin(2/2), que chama a função dec2bin(1/2). Para \(1\), a recursão para e será retornada a string "1", para \(2\) a função retornará a string "10", para \(4\), a string "100" e, para \(8\), a string 1000, correspondendo portanto ao número \(1000\) (return(as.numeric(bin))), ou seja, \(1 \times 2^3 + 0 \times 2^2 + 0 \times 2 ^1 + 0 \times 2^0\).

Por outro lado, para converter um número binário em decimal pegamos cada um dos seus dígitos e o multiplicamos pela base, que é 2, elevado à potência correspondente de acordo com a sua posição e, em seguida somamos os resultados. Observe esse procedimento no vídeo abaixo e, com mais detalhes, no post Sistema de numeração binária.

Discutimos em sala procedimento para a conversão do numeral binário para decimal. Cada dígito (\(0\) ou \(1\)) do número binário deve ser multiplicado pela potência de \(2\), com expoente correspondente a sua posição, como no exemplo mostrado no diagrama abaixo:

Conversão para decimal (Fonte: AGARWAL, 2023.)
Conversão para decimal (Fonte: AGARWAL, 2023.)

O código da função bin2dec() aceita como argumento de entrada váriavel string ou numérica, já que é usada a função strsplit() usada para separar cada carcetere de strings.

bin2dec <- function(binario){
          # transforma digitos do binario em vetor
        X.dig <- as.numeric(strsplit(as.character(binario),"")[[1]])
          # vetor com os expoentes das potências de 2
        potencias <- sort(seq_along(X.dig), decreasing = T)-1
          # conversão binario para decimal
        decimal <- sum(X.dig*2**potencias)
        return(decimal)
}

Dessa maneira se entramos com o número (ou string) binário X <- 11111100111 como argumento da função bin2dec(binario = X), teremos como saída o decimal \(2023\). Essa função não é vetorizada. Uma maneira de trabalhar com ela de maneira vetorizada seria usando a função sapply (ou lapply), como por exemplo para transformar o vetor de binários (11111100111, 110100, 1010), teríamos o seguinte resultado:

    sapply(c(11111100111, 110100, 1010), bin2dec)
## [1] 2023   52   10

Há pacotes no R com funções que realizam, de maneira mais elaborada, as conversões que discutimos aqui.

Algoritmo Quicksort:

Um exemplo clássico de uso de funções recursivas é no algoritmo Quick Sort. Ele é um algoritmo eficiente para ordenar vetores numéricos, do menor para o maior valor. Para fins de compreensão desse procedimento de ordenção utilizaremos o conceito de pivô que é o elemento de um array que, de acordo com algum critério, é selecionado primeiro pelo algoritmo. O procedimento Quicksort baseia-se em particionamentos sucessivos, de acordo com uma rotina que escolhe um pivô e o posiciona no array de uma maneira que os elementos menores ao pivô, formarão a partição esquerda (estão à esquerda do pivô) e os elementos maiores que o pivô, que formarão a partição direita (estão à direita do pivô). Seja por exemplo o vetor (9, -3, 5, 2, 6, 8, -6, 1, 3). Os passos de um procedimento de ordenação por Quick Sorting estão detalhados abaixo, podendo ser acompanhados pelo diagrama “Árvore de recursão do Quicksort”:

  1. Escolher como pivô um dos elementos do vetor. O recomendado é que seja escolhido o elemento mais central como pivô. Algumas implementações selecionam literalmente o elemento mais central, enquanto outros selecionam o primeiro, o último elemento ou mesmo uma escolha aleatória. No exemplo, utilizaremos o último elemento como pivô (3).

  2. Comparar todos os valores com o pivô (3), formando duas partições do vetor: a dos elementos menores que o pivô ((-3, 2, -6, 1)) e a dos elementos maiores que o pivô ((8, 5, 9, 6)).

  3. Aplica-se repetidamente a função nos dois subvetores das partições à esquerda e à direita, particionando cada um deles em cada recursão. As duas primeiras partições ficarão então ordenadas dos menores valores para os maiores, ou sejam, (-6, -3, 1, 2) e (5, 6, 8, 9).

  4. Juntam-se as partições ordenadas com o pivô, ficando (-6, -3, 1, 2), 3 e (5, 6, 8, 9).

Árvore de recursão do Quicksort (Fonte: BIANCHINI et al., 2016.)
Árvore de recursão do Quicksort (Fonte: BIANCHINI et al., 2016.)

A abordagem de função recursiva é utilizada no código da função qsort(), em que, por facilidade de implementação utilizamos como pivô o primeiro elemento de cada partição. Verifique esse código atentamente e acompanhe, no diagrama acima, os resultados de cada partição.

qsort <- function(x) {
        # condição de parada
        if (length(x) <= 1) return(x)
  
        # determinação do pivô
        pivo <- x[1]
        
          # vetor sem o pivô
        resto <- x[-1]
        
          # partição à esquerda
        sv1 <- resto[resto < pivo]
    
          # partição à direita
        sv2 <- resto[resto >= pivo]
        
          # ordenação primeira partição à esquerda
        sv1 <- qsort(sv1)
    
          # ordenação primeira partição à direita 
        sv2 <- qsort(sv2)
        
          # composição das partições e pivô
        return(c(sv1,pivo,sv2))
}

Note que, sem a condição de parada (if (length(x) <= 1) return(x)), a função qsort() chamaria a si mesmo repetidamente, em um loop infinito. Tendo como entrada o vetor xis <- c(9, -3, 5, 2, 6, 8, -6, 1, 3), o comando qsort(xis) apresenta o vetor ordenado desejado, ou seja, (-6, -3, 1, 2, 3, 5, 6, 8, 9). Essa função foi apresentada para exemplicar o uso de recursão em funções em R. A função sort() do R é muito mais rápida e eficiente, pois está escrita em C++.

Cuidados

A técnica de recursão não é usada comumente em análise estatística, mas é importante estar arento a ela. Ela é uma abordagem poderosa, especialmente você não sabe quantas vezes uma função necessita ser chamada para concluir uma tarega. Além disso, a recursão oferece soluções mais rápidas e eficientes para muitos algoritmos de ordenação e busca.

Entretanto, ela apresenra dois incovenientes em potencial:

  • É uma técnica bastante abstrata. Pode-se perceber que ela é apenas o inversa da prova por indução matemática, mas muitos desenvolvedores acham ela difícil.

  • É uma técnica tem um custo computacional elevado, o que pode ser um problema no R quando aplicada a grandes problemas.

Exercícios

Abaixo seguem considerações teóricas e computacionais sobre os conteúdos dos exercícios propostos nesse estudo dirigido.

Suavização Exponencial:

É difícl converter um loop num funcional quando a relação entre os elementos não é independente ou é definida recursivamente. Por exemplo, a suavização exponencial ou (média móvel exponencialmente ponderada), definida formalmente a seguir, funciona tomando uma média ponderada entre o valor atual e os valores anteriores.

Seja \(x_1, x_2, \dots, x_m\) os valores de interesse, observados a intervalos regulares de tempo (contínuo ou discreto). A estatística amostral da média móvel ponderada (EWMA - Exponentially Weighted Moving Average) é definida como:

\[ y_i = \lambda x_i + (1 - \lambda) Y_{i - 1} \text{, com } 0 < \lambda \leq 1 \text{ e } Y_0 = \mu_0 \]

A média móvel exponencialmente ponderada (EWMA) é amplamente usada em modelagem de séries temporais e em previsões e em monitoramento de média de característica de qualidade de processo industrial (\(x_i\) pode ser uma medida individual da característica de qualidade ou uma média de tamanho \(n\), obtida a intervalos regulares).

A expressão é recorrente e pode ser desdobrada, ou seja: \[\begin{align} Y_1 & = \lambda \bar{X}_1 + (1 - \lambda) \mu_0 \\ Y_2 & = \lambda \bar{X}_2 + (1 - \lambda) Y_1 = \lambda \bar{X}_2 + \lambda (1- \lambda)\bar{X}_1 + (1 - \lambda)^2 \mu_0 \\ Y_3 & = \lambda \bar{X}_3 + (1-\lambda)Y_2 = \lambda \bar{X}_3 + \lambda (1-\lambda)\bar{X}_2 + \lambda (1-\lambda)^2 \bar{X}_1 + (1-\lambda)^3 \mu_0 \\ \vdots & = \quad \quad \quad \quad \vdots \\ Y_i & = \lambda \bar{X}_i + (1- \lambda) Y_{i-1} = \lambda \bar{X}_i + \lambda (1- \lambda) \bar{X}_{i-1} + \lambda (1- \lambda)^2 \bar{X}_{i-2}\dots + \lambda(1-\lambda)^{i-1}\bar{X}_1 + (1 - \lambda)^i \mu_0 \end{align}\]

Continuando a substituir recursivamente \(Y_{i-j} \text{, } j = 2, 3, \dots, i\), obtemos:

\[\begin{align} Y_i & = \lambda \sum_{j=0}^{i-1}(1-\lambda)^j \bar{X}_{i-j}+(1-\lambda)^i Y_0 \end{align}\]

Os pesos \(\lambda(1-\lambda)^j \text{, } j = 0, 1, \dots, i-1\) decrescem geometricamente à taxa \((1-\lambda)\), com a idade da média amostral. Se \(\lambda = 0,2\), então o peso associado à média amostral \(i\) é 0.2 \(0,2\) e os pesos dados às médias precedentes são 0.16, 0.128, 0.1024, e assim por diante.

Além disso, os pesos até a i-ésima amostra tem soma:

\[ \sum_{j=0}^{i-1} \lambda (1-\lambda)^j= \lambda\left[\frac{1 - (1-\lambda)^j}{1 - (1-\lambda)} \right]=1-(1-\lambda)^i \]

ou seja, para \(i \to \infty\) os pesos tem soma um. Como esses pesos decrescem geometricamente, \(Y_i\) é também denominada média móvel geométrica.

A função ewma() implementa a média móvel exponencialmente ponderada - EWMA (ou suavização exponencial, no contexto de séries temporais), comn um loop em for. Nessa implementação, adotou-se \(Y_0 = \mu_0\), em que \(\mu_0\) é o valor-alvo da média da característica de qualidade de processo industrial em monitoramento. No contexto de séries temporais, é usual adotar-se \(Y_0 = X_1\), ou seja, pode-se revisar essa linha no código.

ewma <- function(x, lambda, mi0){
        Y <- numeric(length(x) + 1)
        Y[1] <- mi0
        for(i in seq_along(Y)[-1]){
            Y[i] <- lambda*x[i - 1] + (1 - lambda)*Y[i - 1]
        }
        return(Y)
}

Não podemos eliminar o loop em for, porque a saída na posição i depende da da entrada (x[i]), quanto da saída na posição i-1 (Y[i-1]). Pode-se eliminar o loop em for nesse caso, resolvendo, por meio de função recursiva, a relação de recorrência \(y_i = \lambda x_i + (1 - \lambda) Y_{i - 1}\), com critério de parada \(Y_0 = \mu_0\). Isso pode produzir uma função mais simples

Nosso objetivo é aplicar essa função nas observações individuais de característica de qualidade de um processo industrial e que estão disponíveis no conjunto de dados Tab7.1-Costa.txt. São 31 observações, armazenadas no objeto dados e cujos valores são os seguintes:

dados
##  [1] 100.23 100.19 102.02  99.59  99.81  99.86  99.60 100.35  99.38 100.83
## [11]  99.73  98.00  99.72 101.34  98.77  99.94  99.47  99.95 100.33  99.57
## [21] 100.82 100.23 100.68 101.64 100.86  99.28 101.41 100.21 101.85 101.43
## [31] 103.00

Sejam as observações individuais do objeto dados, de uma característica de qualidade cuja média tem valor-alvo 100 e que deseja-se podenderar com fator \(\lambda = 0.2\). Aplicando-se esses argumentos na função ewma() obtemos a média exponencialmente ponderada de cada uma dessas observações, cujo gráfico é dado a seguir:

   # Dados para construção do gráfico de EWMA
alvo <- 100
peso <- 0.2

  # Cálculo da média exponencialmente ponderada
medias.ew <- ewma(x = dados, lambda = 0.2, mi0 = alvo)
medias.ew
##  [1] 100.00000 100.04600 100.07480 100.46384 100.28907 100.19326 100.12661
##  [8] 100.02128 100.08703  99.94562 100.12250 100.04400  99.63520  99.65216
## [15]  99.98973  99.74578  99.78463  99.72170  99.76736  99.87989  99.81791
## [22] 100.01833 100.06066 100.18453 100.47562 100.55250 100.29800 100.52040
## [29] 100.45832 100.73666 100.87532 101.30026
  # Gráfico de EWMA - esboço
plot(medias.ew, type = "o", pch = 19, xlab = "Amostra")
Gráfico de EWMA de processo.

Gráfico de EWMA de processo.

Salienta-se que na expressão de \(Y_i\), percebe-se que, no somatório, os expoentes de \((1-\lambda)\) é um vetor \((0, 1, \dots, i-1)\) e que os índices do termo de \(X_j\) formam o vetor \((i, i-1, \dots, 1)\), ambos os vetores de tamanho \(i\). A média móvel exponencialmente ponderada (EWMA) é usada em modelagem de séries temporais e em previsões. Dessa maneira a EWMA pode ser considerada uma média ponderada de todas as observações passadas e a atual e, por esse motivo, o gráifco de EWMA é insensível a desvios da hipótese de normalidade.

Enunciado:

Baseando-se na equação de recorrência advinda das expressão de \(Y_i\), elabore uma função em R para cálculo da média nóvel exponencialmente ponderada das observações individuais de característica de qualidade, armazenadas no conjunto de dados Tab7.1-Costa.txt. No código deve-se utilizar obrigatoriamente a abordagem recursiva. Compare (e comente) o desempenho de sua função com aquela baseada na abordagem iterativa (loop em for) que foi apresentada nesse texto (função ewma()).

Coeficientes Binomiais:

Os coeficientes binomiais são números inteiros positivos que ocorrem como componentes no teorema binomial, um teorema importante com aplicações em vários algoritmos de aprendizagem de máquinas e em Estatística. Um coeficiente binomial é indexado por um par de inteiros \(n\) e \(k\) e é denotado por \({\displaystyle{\tbinom{n}{k}}}\). Ele é o coeficiente do termo \(x^k\) na expansão polinomial de \((1 + x)^n\). Este coeficiente pode ser pode ser calculado pela fórmula multiplicativa:

\[ {\displaystyle {\binom {n}{k}}={\frac {n\times (n-1)\times \cdots \times (n-k+1)}{k\times (k-1)\times \cdots \times 1}} \text{,}} \]

podendo ser expresso de forma mais compacta com a notação fatorial:

\[ {\displaystyle {\binom {n}{k}}={\frac {n!}{k!(n-k)!}} \text{.}} \] Os coeficientes binomiais podem também ser calculados de uma maneira puramente aditiva, conforme a relação de recorrência indicada abaixo, sendo um exemplo de recursão multinomial:

\[ \binom{n}{k} = \binom{n - 1}{k - 1} + \binom{n - 1}{k} \text{, com } \binom{n}{0} = \binom{n}{n} = 1. \tag{*} \]

O objetivo primeiro é desenvolver uma função recursiva em R, com dois argumentos (C(n, k)) que tenha como saída o coeficiente binomial $ $. O valor de C(n, k) será calculado recursivamente usando a fórmula padrão para os coeficientes binomiais, ou seja, C(n, k) = C(n - 1, k - 1) + C(n - 1, k), com C(n, 0) = C(n, n) = 1. Following is a simple recursive implementation that simply follows the recursive structure mentioned above. Abaixo encontra-se o código da função coefBin(), cuja implementação recursiva é bastante simples, seguindo àquela mecionada estrutura recursiva:

coefBin <- function(n, k){
              # condições de parada
            if(k > n) return(0)
            if(k == 0 || k == n) return(1)

              # recursão
            valor <- coefBin(n - 1, k - 1) + coefBin(n - 1, k)
            return(valor)
}

Por exemplo, sua função deve retornar 6, para n = 4 e k = 2, e deve retornar 10, para n = 5 e k = 2, conforme os resultados a seguir:

coefBin(n = 4, k = 2); coefBin(n = 5, k = 2)
## [1] 6
## [1] 10

Entretanto, note-se que a função acima calcula repetidamente os mesmos subproblemas. levando a um aumento de seu custo computacional [sua complexidade de tempo é da ordem de \(O(n\times \max(k,n-k)\))]. Observe na árvore de recursão para \(n = 5\) e \(k = 2\), a seguir.

Árvore de recursão para coeficiente binomial (Fonte: Geeks for Geeks, 2023.)
Árvore de recursão para coeficiente binomial (Fonte: Geeks for Geeks, 2023.)

Percebe-se que a função \(C(3, 1)\) é chamada duas vezes. Para grandes valores de \(n\), haverá muitos subproblemas comuns. Uma vez que os mesmos subproblemas são chamados novamente, este problema tem a propriedade de sobreposição de subproblemas. Assim, o problema do Coeficiente Binomial tem duas propriedades (Overlapping subproblems e Optinal Structure Property) de um problema de Programação Dinâmica, que, no presente caso, buscará evitar a re-computação dos mesmos subproblemas. Não discutiremos nessa altura do curso o uso de Programação Dinâmica.

A utilização da fórmula (*) no cálculo dos valores de todos os coeficientes binomiais gera um arranjo denominado Triângulo de Pascal.

Triângulo de Pascal (Fonte: PRTI1, 2023.)
Triângulo de Pascal (Fonte: PRTI1, 2023.)

Ele oculta um grande número de padrões, muitos deles descobertos pelo próprio Pascal e até conhecidos antes do seu tempo.

Alguns dos padrões mais importantes do Triângulo de Pascal são:

  • Cada elemento é a soma dos dois números acima dele.

  • Todos os elementos exteriores são iguais a 1.

  • O triângulo é simétrico. Em termos dos coeficientes binomiais \(\binom{n}{k} = \binom{n}{n-k}\).

  • A soma dos elementos da \(n\)-ésima linha do triângulo é \(2^n\). Por outro lado, a soma dos elementos da linha \(n + 1\) é o dobro da soma dos elemntos da linha \(n\). Esse fato pode ser verificado por meio da expansão binomial quando \(x=y=1\), em que:

\[\begin{align} (x + y)^n & = \binom{n}{0}x^n y^0 + \binom{n}{1}x^{n-1} y^1 + \binom{n}{2}x^{n-1} y^2 + \dots + \binom{n}{n-1}x^1 y^{n - 1} + \binom{n}{n}x^0 y^n \\ & = \sum_{k=0}^n \binom{n}{k}x^{n - k} y^k = \sum_{k=0}^n \binom{n}{k}x^k y^{n-k}. \end{align}\]

  • Se o primeiro elemento de uma linha for um número primo (lembre-se, o \(0\)-ésimo elemento de cada fila é 1), todos os números dessa linhas (excluindo-se os 1’s) são divisíveis por ele. Por exemplo, na 7ª linha \((1, 7, 21, 35, 35, 21, 7, 1)\), todos os números diferentes de 1 são divisíveis por 7.

  • A sequência de Fibonacci pode ser obtida somando dos números na diagonal do Triângulo de Pascal.

Sequência de Fibonacci no triângulo de Pascal (Fonte: PRTI1, 2023.)
Sequência de Fibonacci no triângulo de Pascal (Fonte: PRTI1, 2023.)

Você encontrará, nos posts indicados a seguir, maiores detalhes sobre o Triângulo de Pascal, suas propriedades e padrões: Patterns in Pascal’s Triangle e Pascal’s Triangle and Its Patterns.

Salienta-se que, além da sobreposição de subproblemas descrita acima, há cálculos redundantes em decorrência da propriedade de simetria dos coeficientes binomiais. De fato \(\binom{n}{k} = \binom{n}{n - k}\), podendo-se provar por indução finita a partir da definição dada em (*), ou seja considerando que a hipótese de indução é verdadeira verdadeira para \(n\), ou seja, \(\binom{n}{k} = \binom{n}{n - k}\). Precisa-se provar que \(\binom{n + 1}{k} = \binom{n + 1}{n + 1 - k}\). Por definição (*) \(\binom{n + 1}{k} = \binom{n}{k - 1} + \binom{n}{k}\), com a hipótese de indução de que \(\binom{n}{k} = \binom{n}{n - k}\), temos que:

\[\begin{align} \binom{n + 1}{k} & = \binom{n}{k} + \binom{n}{k - 1} = \binom{n}{n - k} + \binom{n}{n - (k-1)} \\ & = \binom{n + 1}{n + 1 - k} \end{align}\]

Uma alternativa para o cálculo dos coeficientes binomiais é por meio de outra formulação, que também leva a uma recorrência unidimensional. Ela é obtida a partir da formulação por fatoriais do coeficiente binomial, ou seja:

\[\begin{align} \tag{I} \binom {n}{k}& = \frac {n!}{k!(n-k)!} \\ & = \frac{n}{k} \binom{n-1}{k-1} \text{, com } \binom{1}{0} = 1 \end{align}\]

Outra possibilidade de cálculo dos coeficientes binomiais utiliza sua formulação multiplicativa, levando também a uma recorrência unidimensional:

\[\begin{align} \tag{II} \binom{n}{k} & =\frac {n\times (n-1)\times \cdots \times (n-k+1)}{k\times (k-1)\times \cdots \times 1} \\ & = \binom{n}{k - 1} (n - k + 1)/k \text{, com } \binom{n}{0} = 1. \end{align}\]

As recorrências oriundas das expressões (I) e (II) são bastante usadas computacionalmente, porque elaa não exigem a construção de uma tabela, como no caso da recorrência bidimensional e nem envolvem inteiros muito grandes, como na aplicação pura e simples da fórmula fatorial, embora ainda possa ocorrer sobreposição de subproblemas.

Entretanto, uma maneira de evitar a sobreposição de chamadas de funções, situação oriunda da abordagem recursivo, é usar um método não-recursivo, como quando efetuamos o cálculo de maneira iterativa. Para simplificar a iteração considere a seguinte decomposição da expressão do coeficiente binomial:

\[\begin{align} \binom{n}{k} & = \frac {n!}{k!(n-k)!} = \frac{(k+1) \times (k+2) \times \dots \times n}{ 1 \times 2 \dots \times (n-k)} \\ & = \frac{k+1}{1} \times \frac{k+2}{2} \times \dots \times \frac{k + (n-k)}{n-k} \\ & = \prod_{i=1}^{n-k} \frac{k + i}{i} = \prod_{i=1}^{n-k} \left( 1 + \frac{k}{i} \right) \end{align}\]

Abaixo encontra-se um esboço de código para cálculo do coeficiente binomial usando um método iterativo usual (loop em for). Para evitar loops longos, usa-se a propriedade de que \(\binom{n}{k} = \binom{n}{n-k}\), ou seja, nesse caso, troca-se, na expressão acima, \(k\) por \(n-k\) e vice-versa:

\[ \binom{n}{n - k} = \prod_{i=1}^{k} \frac{(n-k) + i}{i} = \prod_{i=1}^{k} \left( 1 + \frac{n-k}{i} \right) \]

        # ------- Abordagem não-recursiva -------

coefBin.iter <- function(n, k){
              # numerador da iteração
            if(n - k < k) num <- k else num <- (n - k)
              # iteração para cálculo do produtório
            res <- 1
            for(i in 1:(n - num)){ res <- res*(1 + num/i)}
            return(res)
}

Como na aplicação da função anterior, a função coefBin.iter() deve retornar 6, para n = 4 e k = 2, e, 10, para n = 5 e k = 2, conforme os resultados a seguir:

coefBin.iter(n = 4, k = 2); coefBin.iter(n = 5, k = 2)
## [1] 6
## [1] 10

As funções do R Base choose() e lchoose() tem como saída os coeficientes binomiais e seus logaritmos, respectivamente. Note que choose(n, k) é a função relevante em R para calcular coeficientes binomiais, sendo definida para \(n \in \mathbb{R}\) e \(k \in \mathbb{Z}\). Para \(k \geq 1\) é definido como \(n(n - 1) \dots (n - k + 1)/k!\), como \(1\), para \(k = 0\) e como \(0\), para \(k < 0\). Os valores não-inteiros de k são arredondados para um número inteiro, com um aviso. Por razões de velocidade e precisão, para pequenos valores de k, a função choose(n, k) usa aritmética direta (em vez de chamar a função [l]gamma()). Note-se que a função combn{utils} enumera todas as combinações possíveis. As soluções vetorizadas para o coeficiente binomial e para seu logaritmo são bastante simples e são desenvolvidas a partir da formulação multiplicativa de coeficiente binomial, ou seja:

\[ \binom {n}{k} = \frac{n \times (n-1) \times \cdots \times (n-k+1)}{k \times (k-1) \times \cdots \times 1} \]

Perceba que tanto o denominador quanto o numerador têm \(k\) termos, assim, dados n e k, construímos o vetor vet.k contendo os inteiros entre 1 e k, calculando o coeficiente binomial com a função prod(). O código abaixo é uma possibilidade inicial, lembrando que, por definição, \(\binom{n}{0} = 1\).:

  # valores do coeficiente binomial
n <- 5; k <- 2

  # vetor com os k elementos
vet.k <- 1:k

  # cálculo do coeficiente binomial
prod((n - vet.k + 1)/vet.k)
## [1] 10

Por essa formulação, pode-se calcular, por exemplo, \(\binom{1000}{500}\), que não conseguimos obter se tentamos calcular por meio da fórmula com fatoriais, pois ocorre overflow (calcule factorial(500)).

  # valores do coeficiente binomial
n <- 1000; k <- 500

  # vetor com os k elementos
vet.k <- 1:k

  # cálculo do coeficiente binomial
prod((n - vet.k + 1)/vet.k)
## [1] 2.702882e+299

Esse valor é o mesmo que aquele obtido por choose(1000, 500). Esse é um motivo para termos uma transformação por logaritmos, que pode ser muito eficiente computacionalmente quando trabalhamos com valores muito grandes:

\[\begin{align} \binom {n}{k} & = \frac {n!}{k!(n-k)!} = \frac{ \prod_{i=1}^n i}{\prod_{i = 1}^k i \prod_{i = 1}^{n - k} i} \\ \log \left(\binom {n}{k}\right) & = \sum_{i=1}^n \log(i) - \sum_{i = 1}^k \log(i) - \sum_{i = 1}^{n - k} \log(i) \end{align}\]

O esboço de código abaixo é uma aplicação da expressão acima, criando o vetor log.vet <- log(1:n). assim

        # ------- Cálculo vetorizado com transformação log ---------

  # valores do coeficiente binomial
n <- 1000; k <- 500

  # vetor com os k elementos
vet.log <- log(1:n)

  # cálculo do log do coeficiente binomial
log.Coef <- sum(vet.log) - sum(vet.log[1:k]) - sum(vet.log[1:(n-k)])
log.Coef
## [1] 689.4673
  # Cálculo do coeficiente binomial
exp(log.Coef)
## [1] 2.702882e+299

Verifiquem as funções do R Base factorial(n), lchoose(n, k), gamma(x) e lgamma(x) que também estão relacionadas com o cálculo de fatoriais e coeficientes binomiais.

Enunciado:

Baseando-se na equação de recorrência advinda das expressões (I) OU (II), elabore uma função em R para cálculo de coeficiente binomial definida para \(n \in \mathbb{R}\) e \(k \in \mathbb{Z}\). No código deve-se utilizar obrigatoriamente a abordagem recursiva. Compare (e comente) o desempenho de sua função com aquela baseada em recursão multinomial e com a função choose(), do R Base.

Referências:

Youtube fibonacci