Root-Finding: O método de Newton-Raphson

Acho que a primeira coisa vem a cabeça, depois de ler o título desse post e lembrando que a gente está num blog de um biólogo é…

“Porque diabos eu quero saber isso?”

É uma boa pergunta, mas veja, quanto mais a gente vai estudando alguma coisa, algum sistema, a gente vai refinando nosso conhecimento, a gente vai fazendo predições melhores, mais específicas, e predições melhores implicam em melhores modelos, e algumas vezes ainda a gente vai cair em modelos não lineares, infelizmente, nem tudo é uma reta. Eu mesmo ja cai aqui nesse caso.

Mas beleza, e qual a relação com esse negócio de achar raiz?

Se você for usar a função nls do R por exemplo, ela tem um argumento chamado algorithm, de um ?nls para ver o help.

Usage nls(formula, data, start, control, algorithm, trace, subset, weights, na.action, model, lower, upper, ...)

Se a gente olhar esse argumento , a gente ve que

algorithm - character string specifying the algorithm to use. The default algorithm is a Gauss-Newton algorithm. Other possible values are "plinear" for the Golub-Pereyra algorithm for partially linear least-squares models and "port" for the ‘nl2sol’ algorithm from the Port library

Temos esse Gauss-Newton como o método default, e se você olhar sobre ele no wikipedia aqui, vai ver a primeira frase.

“It can be seen as a modification of Newton’s method for finding a minimum of a function.”

Na minha tradução nórdica livre, “ele pode ser visto como uma modificação do método de Newton para achar o mínimo de uma função.”

Esse método de achar por sua vez é uma extensão do método de Newton-Raphson.

Newton acho que não precisa de introduções, simplesmente o cara do da teoria da gravidade, e temos esse tal de Raphson, que publicou um ano antes do Newton o método, mas o Newton desenvolveu o método independentemente, e para uso mais genérico.

Então estamos no ponto de partida aqui, entender aqui vai nos permitir entender melhor os métodos de optimização mais complexos, o nosso objetivo, e a hora que a gente pegar o jeito aqui, usar o optmin do R pode ficar mais fácil, e plausível. E usar o optmin, assim como métodos de otimização em geral abrem várias oportunidades de fazer coisas legais. E aqui você pode ler, testar as hipóteses que realmente nos interessam em biologia, e não testar o que da com os dados! Não podemos prender as idéias das pessoas a modelos pre-definidos meramente porque a gente não sabe expandir ou melhorar estes.

Mas vamos começar a embrenhar nesses algorítimos, que no fim das contas são idéias bem legais e úteis.

Para começar, vamos tentar entender o método de Newton-Raphson para achar a raiz de uma função.

Antes de tudo vamos definir melhor nosso objetivo.

Nosso objetivo é achar a raiz de uma função, mas o que é a raiz de uma função?
A raiz de uma determinada função f(x) é o valor de x para o qual f(x)=0, simples assim.

Então existe uma função f qualquer a qual entra um valor real e sai um valor real, f:\mathbb{R} \to \mathbb{R}, e para essa função pode existir um a \in \mathbb{R} tal que f(a)=0, ou sempre existe, não tenho certeza se posso afirmar isso, mas não vamos parar aqui.

Então quando a gente faz aquele velho gráfico, colocando a resposta no eixo y e o preditor no eixo x, a raiz é o lugar do preditor onde a resposta fica 0.

Agora eu concordo que até agora falamos somente blablabla. Porque o interesse nisso?
Suponha que temos uma dívida que tem um valor P, um taxa de juros mensal de r e uma duração de n meses pagar, e um pagamento mensal no valor de A. O valor da dívida após n pagamentos, ou n meses, será P_n.

Então, pensando em tudo isso como uma formula, temos o seguinte:

P_0=P

Primeiro, a gente tem uma dívida inicial P dai a recorrência é a seguinte.

P_{n+1}=P_n \cdot (1+r)-A

Isso é, a cada mês, você paga o juros no referente a quanto você devia mês passado, e então reduz da quantidade da dívida a quantidade A que você sempre paga, essa é uma equação de recorrência de primeira ordem, e tem a seguinte solução:

{P_{n}}={{P \cdot(1+r)^n} -  {{A\cdot((1+r)^n-1 )}\over{r}}}

Colocando P_N=0, ou seja, a divida acabou, temos o seguinte.

{A\over{P}}={r \cdot (1+r)^n \over{(1+r)^n-1}}

Dessa equação, podemos isolar parâmetros, como a taxa de juros (r) por exemplo, suponha que sabemos P, n e A, podemos achar r achando a raiz da seguinte função.

f(x) ={{A\over{P}}-{x \cdot (1+x)^n \over{(1+x)^n-1}}}

Agora acho que ja da para ter uma idéia, que achar raízes tem a ver com otimizar parâmetros certo? Tem a ver não, na verdade, mas pode servir de caminho.
Primeiro a gente vai ver como o método de Newton-Raphson funciona achando raízes, e depois podemos modificar ele um pouco para começar estimar parâmetros em modelos.

Mas vamos la. Agora vamos ver como é que o Newton-Raphson funciona.

Suponha que temos uma função f que é diferenciável com uma derivativa continua f\prime e uma raiz a. Agora se temos um x_0 \in \mathbb{R} como nosso chute inicial de a, a linha entre o ponto (x_0,f(x_0)) com inclinação f\prime(x_0) é a melhor linha aproximação de linha reta da função f(x) para o ponto x_0, esse é o significado da derivativa, a equação para essa linha reta é dada por

{f\prime(x_0)}={{f(x_0)-y}\over{x_0-x}}

Agora, essa linha reta cruza o eixo-x até o ponto x_1, que deve ser uma melhor aproximação que x_0 para a, para achar x_1 nos fazemos

{f\prime(x_0)}={{f(x_0)-0}\over{x_0-x_1}}

logo, a gente tem para x_1 que

 {x_1} = {x_0 - {f(x_0) \over{f\prime(x_0)}}}

De forma mais simples, ja da para sacar o que o algorítimo vai fazer, a gente vai ficar subtraindo  f(x) \over{f\prime(x)} do nosso chute para gerar outro chute, até não aguentar mais.

Ou seja, apos achar x_1, a gente vai atrás do x_2, que vai ser:

 {x_2} = {x_1 - {f(x_1) \over{f\prime(x_1)}}}

Ou de forma geral, podemos escrever esse método da seguinte forma.

 {x_{n+1}} = {x_n - {f(x_n) \over{f\prime(x_n)}}}

Uma relação de recorrência.

Pode ser demonstrado (não que eu consiga) que se f é “bem comportado” em a (de uma olhada aqui para ver o caos, mas se  f\prime(a) = 0 e  f\prime\prime(a) é finito e contínuo em a, vai dar certo) e você começar com x_0 perto o suficiente de a, então x_n vai convergir rápido. Infelizmente, quase sempre a gente não conhece o comportamento de f, a gente pensa nos modelos primeiro, depois vai estudar eles, então a gente não sabe a princípio se esse método vai funcionar ou não.

O método de Newton-Raphson não é garantido de te dar uma resposta, você não tem garantias de uma solução, entretanto, se x_n \to a, então f e f\prime são contínuos e nos temos

{a} = {\lim \limits_{n\to\inf} x_{n+1}} = {\lim \limits_{n\to\inf} ({x_n - {f(x_n) \over{f\prime(x_n)}}})} = {\lim \limits_{n\to\inf} x_n - {f(\lim \limits_{n\to\inf} x_{n}) \over{f\prime(\lim \limits_{n\to\inf} x_{n})}}} = {{a} - {f(a)\over{f\prime(a)}}}

Veja que pode parecer um monte de formula complexa, mas a gente ta substituindo a formula dali de cima dentro dela, dentro dela denovo e denovo, recorrentemente.
Então, desde que f\prime(a) \ne \pm \inf, nos devemos ter f(a) = 0.
Como nos estamos esperando f(x_n) \to 0, um boa condição de parada para o algorítimo é quando \vert f(x_n) \vert \leq Q para uma tolerância Q que nos determinarmos.
Se a sequência \{x_n\}_{n}^{\inf} = {0} está convergindo para a raiz a então para x próximo de a temos que f(x) \approx f(a)\cdot(x-a), então se \vert f(x_n) \vert \leq \epsilon então {\vert (x-a) \vert } \leq {\epsilon\over{f\prime(a)}} aproximadamente.

Agora, depois de todas essa formulas, vamos ver o algorítimo funcionando, e talvez vendo ele funcionando, a gente melhore a nossa crença que as coisas ali em cima funcionam, ou melhor que isso, vendo funcionando a gente melhora nosso entendimento.

newtonraphson <- function(ftn, x0, tol = 1e-9, max.iter = 100) {
    x <- x0
    fx <- ftn(x)
    iter <- 0
 
    while ((abs(fx[1]) > tol) && (iter < max.iter)) {
        x <- x - fx[1]/fx[2]
        fx <- ftn(x)
        iter <- iter + 1
        cat("Na iteração", iter, "o valor de x foi de", x, "\n")
    }
    # A saida depende do que acontece
    if (abs(fx[1]) > tol) {
        cat("O algoritimo falhou em convergir\n")
        return(NULL)
    } else {
        cat("O algoritimo convergiu\n")
        return(x)
    }
}

E é isso ai, pequenininho. O que é nossa entrada? A função, oras bolas, um valor inicial pra começar a recorrência, uma tolerância, pra gente saber onde parar e um número máximo de iterações, isso porque como a gente disse, a gente não tem certeza se toda função tem raiz. Então a gente vai tentar esse tanto de vezes ai, se num deu, a gente para o algorítimo e abandona, vamos pra outro método.

Ai olha a nossa função, a gente começa com

x <- x0
fx <- ftn(x)
iter <- 0

A gente começa o valor de x como x_0
Salvamos em fx, a saída da nossa função, que vai ser uma função que vai calcular o valor de f(x) e o valor de f\prime(x) e devolver o resultado num vetor de dois valores. Lembrando que o R pode facilitar encontrar a derivada com a função deriv, de uma olhada aqui para ver como faz.
E a entrada é a iteração zero.

Dai pra frente é exatamente a mesma coisa.

while ((abs(fx[1]) > tol) && (iter < max.iter)) {
        x <- x - fx[1]/fx[2]
        fx <- ftn(x)
        iter <- iter + 1
        cat("Na iteração", iter, "o valor de x foi de", x, "\n")
    }

Olha so a gente calcula o próximo x com a nossa formulinha de recorrência, dai usando esse valor novo de x e a gente usa de novo nossa função para pegar novos valores de f(x) e f\prime(x) e vamos para uma nova iteração.
A gente vai fazer isso até que ou tenhamos ultrapassado o valor máximo de iterações (para não correr o risco de fazer um loop que não para nunca) ou até não ultrapassarmos nossa tolerância entre as iterações.

Agora vamos fazer um exemplo aqui pra ver se vai dar certo.

Primeiro pegamos uma equação.

f(x)=log(x) - exp(-x)

Bem, depois de uma aulinha de calculo 1 na faculdade ou a matéria de calculo do coursera ou olhar muitos videos do khan academy, a gente sabe que

{f\prime(x)}={{1\over{x}} + {exp(-x)}}

E passando para a linguagem do R a gente faz o seguinte.

funcao <- function(x) {
    fx <- log(x) - exp(-x)
    dfx <- 1/x + exp(-x)
    return(c(fx, dfx))
}

Veja que a entrada é o valor de x, e a saída é um vetor de dois valores, o f(x) e f\prime(x) respectivamente.

Usando nossa o método do Newton-Raphson temos o seguinte.

> newtonraphson(funcao, 2, 1e-06) Na iteração 1 o valor de x foi de 1.12202 Na iteração 2 o valor de x foi de 1.294997 Na iteração 3 o valor de x foi de 1.309709 Na iteração 4 o valor de x foi de 1.3098 O algoritimo convergiu [1] 1.3098

Bang, f(1.3098) \approx 0

A gente pode olhar num gráfico isso pra ficar mais fácil.

01

Olha ai, quando o x é 1.3098, no y a gente ta no zero.

A gente pode colocar esse número na função que fizemos também.

> funcao(1.3098) [1] 4.280091e-07 1.033349e+00

Lembrando que o f(x) é o primeiro número e esta em notação científica, ou seja

> format(4.280091e-07,scientific=F) [1] "0.0000004280091"

E isso é bem próximo de zero, se quiser mais precisão, é só diminuir a tolerância, mais ai precisara de mais iterações, para mais precisão, e isso depende na verdade da aplicação que a gente vai dar para esse número. E em ultima instancia, estamos limitados ao menor número que o computador pode representar, e isso não é infinito.

Agora veja como é a derivativa da função

02

Veja que a derivativa nos da a direção na qual mudar o valor de x, a gente sabe se o valor de x ta mudando, e como está mudando.

Veja que imprimimos os valores que fomos usando, então podemos olhar como foram as iterações.
Ou é so modificar a função e guardar em algum lugar essa informação, transformando a saída numa lista por exemplo.

03

Veja que começamos no 2, ai diminuímos bastante, e fomos voltando, cada vez mais devagarinho, veja que as setas vão diminuindo, até chegar num momento que não da mais nem para desenhar setinhas.

Certo, mas vai que foi sorte, vamos tentar para outra função, a boa e velha x^2, pra ela temos o seguinte

{f(x)}={x^2}

e a derivada é fácil, é

{f\prime(x)}={2 \cdot x}

Passamos essa função pro R, pra poder usar o método de Newton-Raphson

funcao <- function(x) {
    fx <- x^2
    dfx <- 2*x
    return(c(fx, dfx))
}

e vamos usar ela

> newtonraphson(funcao, 10, 1e-06) Na iteração 1 o valor de x foi de 5 Na iteração 2 o valor de x foi de 2.5 Na iteração 3 o valor de x foi de 1.25 . . . Na iteração 12 o valor de x foi de 0.002441406 Na iteração 13 o valor de x foi de 0.001220703 Na iteração 14 o valor de x foi de 0.0006103516 O algoritimo convergiu [1] 0.0006103516

Veja que a diferença da iteração 13 para 14 é menor que a nossa tolerância

> format(1e-06,scientific=F) [1] "0.000001"

Veja o último número que muda durante as iterações.
Agora o lado bom dessa função é que a gente tem outro método pra achar raiz, para testa o que estamos usando, podemos usar a boa e velha bhaskara.

Se a gente usa bhaskara, a gente acha as duas raízes, e ambas são zero, veja que a entrada da função são os coeficientes, a, b e c.

> bhaskara(1,0,0) [1] 0 0

A gente pode olhar a solução também num gráfico como no exemplo anterior.

04

Mas e se tivesse duas raízes? Vamos pegar uma segunda função de segundo grau aqui, que dai a gente pode conferir com bhaskara.

Vamos usar essa função

f(x)=2x^2+4x-4

a derivada dela é bem fácil também

f\prime(x)=4x+4

E no R fica assim

funcao <- function(x) {
    fx <- 2*x^2+4*x-4
    dfx <- 4*x+4
    return(c(fx, dfx))
}

Primeiro vamos ver onde queremos chegar usando bhaskara

> bhaskara(2,4,-4) [1] 0.7320508 -2.7320508

Dai vamos usar o método de Newton-Raphson.

> newtonraphson(funcao, 10, 1e-06) Na iteração 1 o valor de x foi de 4.636364 Na iteração 2 o valor de x foi de 2.084311 Na iteração 3 o valor de x foi de 1.028488 Na iteração 4 o valor de x foi de 0.753711 Na iteração 5 o valor de x foi de 0.7321846 Na iteração 6 o valor de x foi de 0.7320508 O algoritimo convergiu [1] 0.7320508

E veja que conseguimos a primeira raiz, mas veja que começamos chutando o valor de 10, dai encontramos a raiz mais próxima desse chute inicial, isso quer dizer que se chutarmos um número mais próximo a outra raiz, acharemos ela

> newtonraphson(funcao,-10, 1e-06) Na iteração 1 o valor de x foi de -5.666667 Na iteração 2 o valor de x foi de -3.654762 Na iteração 3 o valor de x foi de -2.892403 Na iteração 4 o valor de x foi de -2.738845 Na iteração 5 o valor de x foi de -2.732064 Na iteração 6 o valor de x foi de -2.732051 O algoritimo convergiu [1] -2.732051

E pimba, temos as duas raízes, e com a certeza que o algorítimo funciona beleza, já que já sabíamos o que esperar.
Podemos também ver o método no gráfico.

05

Mas então porque não usar sempre bhaskara? Ai a gente cai na mesma pergunta, de porque não usar sempre funções lineares nos modelos, porque as vezes não da.
Vamos a um último exemplo, para as funções esquisitas.

Pense na função

f(x)=exp(x)

espantosamente, a derivada dela é

f\prime(x)=exp(x)

Não vamos debater isso, mas é legal né, a derivada é a própria função, aqui bhaskara não funciona, mas Newton-Raphson pode dar certo.

Estabelecemos nossa função no R

funcao <- function(x) {
    fx <- exp(x)
    dfx <- exp(x)
    return(c(fx, dfx))
}
> newtonraphson(funcao, 10, 1e-06,100) Na iteração 1 o valor de x foi de 9 Na iteração 2 o valor de x foi de 8 Na iteração 3 o valor de x foi de 7 . . . Na iteração 23 o valor de x foi de -13 Na iteração 24 o valor de x foi de -14 O algoritimo convergiu [1] -14

E, pimba, deu certo também.

Podemos olhar no gráfico para ver se esse número faz sentido

06

E parece que está tudo certinho.

E esse é o método Newton-Raphson. Como a gente já viu no algorítimo do gradient descent, daqui, com algumas modificações, a gente vai parametrizar modelos também.

Agora imagine que o Newton, um cara que olhava planeta, tinha dados de posição de planeta e parametrizava as coisas fazendo contas na mão, ainda bem que esse método não precisa de sei la, 10 mil iterações, senão ele tava ferrado hehehe.

Mas é isso ai, o primeiro algorítimo de encontrar raízes que vemos, num próximo post nos progredimos ele um pouco para determinar parâmetros.

Bem, eu vou ficar por aqui, que já faz muito tempo que estou escrevendo aqui. Scripts no repositório recologia e aqui em baixo, e até o próximo post, este aqui espero dar sequência e usar o método de Newton para achar o mínimo de funções, ou seja, estimar parâmetros pela máxima verossimilhança.

Referencia:
Owen Jones, Robert Maillardet & Andrew Robinson 2009 – Introduction to Scientific Programming and Simulation Using R. CRC Press 499pp

bhaskara<-function(a,b,c) {
    delta = b*b - 4*a*c
    if(delta < 0) {
        print("A equacão não possui raizes reais.\n");
        saida<-c(NA,NA)
        return(saida)
    }else{
        x1 = (-b + sqrt(delta)) / (2*a)
        x2 = (-b - sqrt(delta)) / (2*a)
        saida<-c(x1,x2)
        return(saida)
    }
}
 
newtonraphson <- function(ftn, x0, tol = 1e-9, max.iter = 100) {
    x <- x0
    fx <- ftn(x)
    iter <- 0
 
    while ((abs(fx[1]) > tol) && (iter < max.iter)) {
        x <- x - fx[1]/fx[2]
        fx <- ftn(x)
        iter <- iter + 1
        cat("Na iteração", iter, "o valor de x foi de", x, "\n")
    }
    # A saida depende do que acontece
    if (abs(fx[1]) > tol) {
        cat("O algoritimo falhou em convergir\n")
        return(NULL)
    } else {
        cat("O algoritimo convergiu\n")
        return(x)
    }
}
 
##########################
#Exemplo 01
##########################
funcao <- function(x) {
    fx <- log(x) - exp(-x)
    dfx <- 1/x + exp(-x)
    return(c(fx, dfx))
}
 
newtonraphson(funcao, 2, 1e-06)
 
funcao(1.3098)
 
format(4.280091e-07,scientific=F)
 
#figura 1
curve(log(x) - exp(-x),0,3,frame=F)
points(1.3098,funcao(1.3098)[1],cex=2)
abline(v=1.3098,h=funcao(1.3098)[1],lty=3,col=2)
legend("topleft",legend="Raiz",pch=1,bty="n")
 
#figura 2
jpeg("02.jpg")
par(mfrow=c(2,1))
curve(log(x) - exp(-x),0,3,frame=F,main="Função")
abline(v=1.3098,h=funcao(1.3098)[1],lty=2)
curve(1/x + exp(-x),0,3,frame=F,main="Primeira derivada")
abline(v=1.3098,lty=2)
 
#figura 3
pontos<-c(2,1.12202,1.294997,1.309709,1.3098)
plot(pontos,log(pontos)-exp(-pontos),type="p",pch=19,cex=0.5,frame=F,xlab="x",ylab="f(x)")
text(pontos[1]-0.1,(log(pontos)-exp(-pontos))[1],"inicio")
#ps, a ultima seta não tem tamanho suficiente para o gráfico
for(i in 1:5) {
    arrows(pontos[i],(log(pontos)-exp(-pontos))[i],pontos[i+1],(log(pontos)-exp(-pontos))[i+1])
}
 
###########################3
#Exemplo 2
############################
funcao <- function(x) {
    fx <- x^2
    dfx <- 2*x
    return(c(fx, dfx))
}
 
newtonraphson(funcao, 10, 1e-06)
format(1e-06,scientific=F)
 
bhaskara(1,0,0)
 
#figura 4
par(mfrow=c(2,1))
curve(x^2,-5,5,frame=F,main="Função")
abline(v=0,h=0.0006103516,lty=3,col=2)
curve(2*x,-5,5,frame=F,main="Primeira derivada")
abline(v=0.0006103516,lty=3,col=2)
 
 
###################
#Exemplo 3
###################
funcao <- function(x) {
    fx <- 2*x^2+4*x-4
    dfx <- 4*x+4
    return(c(fx, dfx))
}
 
bhaskara(2,4,-4)
 
newtonraphson(funcao, 10, 1e-06)
newtonraphson(funcao,-10, 1e-06)
 
#figura 5
par(mfrow=c(2,1))
curve(2*x^2+4*x-4,-10,10,frame=F,main="Função")
abline(h=0,lty=3,col=2)
abline(v=0.7320508,lty=3,col=2)
abline(v=-2.732051,lty=3,col=2)
curve(4*x+4,-10,10,frame=F,main="Primeira derivada")
abline(v=0.7320508,lty=3,col=2)
abline(v=-2.732051,lty=3,col=2)
 
###############
#exemplo 4
###############
 
funcao <- function(x) {
    fx <- exp(x)
    dfx <- exp(x)
    return(c(fx, dfx))
}
 
newtonraphson(funcao, 10, 1e-06,100)
 
#figura 6
par(mfrow=c(2,1))
curve(exp(x),-15,-5,main="Função exp(x)",frame=F)
abline(h=0,v=-14,lty=3,col=2)
curve(exp(x),-15,-5,main="Primeira derivada de exp(x)",frame=F)
abline(v=-14,lty=3,col=2)

Rosalind e o primeiro problema em R e Python

Eu sempre tive um interesse especial por evolução, desde os livros do Richard Dawkins e do Stephen Jay Gould.

E a algum tempo eu comecei a olhar o Project Euler, que é uma proposta muito legal, o projeto é um site, que continuamente vai colocando novas perguntas, ou desafios de programação se preferir, para as pessoas solucionarem.

Você pode solucionar como quiser, até a mão, mas quase sempre é necessário programar algo, uma função ou algoritmo para conseguir chegar a uma resposta. Mas você pode usar a linguagem de programação que desejar. E eu fui la tentar minha sorte com R. O primeiro impacto foi ver como eu tenho uma deficiência enorme em matemática, mas devagarzinho da para ir resolvendo as coisas. O legal que após solucionar um desafio, você tem acesso a um fórum sobre esse desafio, onde pode ver como as pessoas solucionaram e como a sua solução é inferior e menos elegante hehehe. E tem acesso também ao que seria o método esperado para solucionar o problema pelos organizadores do projeto.

Mas a partir do Project Euler eu conheci o Project Rosalind, que segue a mesma idéia, mas é voltado para problemas de bioinformática. Então juntando R com evolução, é algo muito legal.

Tentando resolver os problemas com R, que alias já tem soluções, no pacote ape por exemplo, mas pensando em entender melhor como as coisas funcionam, fui fazer minha função. Mas conforme fui solucionando os problemas vi que tanto no Project Euler como no Project Rosalind, as pessoas gostam da linguagem Python.

O John Cook mesmo (que sigo os mil twitters dele o excelente blog) esta sempre postando coisas com código de Python (Ele também usa R, e já foi o editor do R Journal, nos seus primórdios), então resolvi conciliar o Project Rosalind com tentar aprender outra linguagem que é Python.

Mas só para dar um exemplo vamos ao primeiro desafio.

Ele é fazer uma função, programa, algoritmo, ou sei la como chamar, mas algo que pegue uma sequência de bases e conte o número de ocorrências de “a”, “c”, “t” e “g”.

No R eu peguei fiz uma função que só entra um argumento, a sequência como character, ai a função apenas pica essa sequência em um vetor onde cada letra é um item e depois usa table() para fazer a contagem.

A função ficou assim:

cont.seq<-function(seq) {
seq.split<-strsplit(seq,"")
print(table(seq.split))
}

Então ela funcionou assim:

> cont.seq("AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC") seq.split A C G T 20 12 17 21

Bem tem alguns defeitos, como poderíamos mudar o nome da tabela para não sair “seq.split”, que eu defini isso dentro da função.
Mas no final das contas elas resolveu o serviço.
Minha outra solução ainda foi usando Regular Expressions, mas isso fica para outro post.

Mas o legal que para python a ideia vai ser a mesma, fazer uma função.
Então eu defino minha função:

[sourcecode language=”python”]
def contseq(s):
return s.count("A"), s.count("G"), s.count("C"), s.count("T")
[/sourcecode]

Como eu ainda não peguei o jeito de criar botões, telinhas para entrar dados ou sair o resultado, vai tudo pelo interpretador mesmo. Igual o R
Então o programa inteiro fica assim:

def contseq(s):
    return s.count("A"), s.count("G"), s.count("C"), s.count("T")
 
sequencia = "AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC";
print contseq(sequencia)

E temos como resultado ao rodar o programa:

>>> (20, 17, 12, 21)

Ohhhh, o mesmo resultado do R 🙂

Bem é isso, meu primeiro programinha em Python.
Quem quiser se aventurar em Python comigo ainda esta rolando essa disciplina no coursera sobre python:
https://www.coursera.org/course/interactivepython

Da uma boa ideia inicial de como funciona as coisas em Python. Ele tem alguma similaridades com R, e se não me engano, funções que você fizer em python, existe um pacote no R para rodar ou adaptar elas, não entendo muito bem.

E como ultimo estimulo, um quadrinho muito legal do xkcd.

python