Crescimento continuo em populações.

“A solução é igual, só não tão retinha quanto simplesmente ligar os pontos, porque a cada momento alguém esta nascendo, a população esta mudando continuamente, ou seja aqui não estamos dando passos discretos, mas passos contínuos.”

E aqui uma frase do ultimo post. Na verdade isso é algo que sempre me trouxe dificuldade, porque aquela linha não passa exatamente sobre os pontos?

Mas primeiro vamos pensar somente no crescimento não dependente de densidade.

Então temos que o crescimento populaticonal é:

 N_{t+1}= N_t\cdot(1+r)\cdot (1-{N_t\over K})

Essa equação é a que gente montou no R no ultimo post, única diferença é que trocamos o k por alpha, que seria isso.

 N_{t+1}= N_t\cdot(1+r)\cdot (1-{N_t\cdot \alpha})

Ao invés de pensar em o quanto cabe de indivíduos, k, a gente pode pensar o quanto cada individuo incomoda o outro, 1/k, que é alpha. Bem não tem muito milagre nessa mudança.

Mas então como queremos somente o crescimento independente da densidade, queremos somente essa primeira metade:

 N_{t+1}= N_t\cdot(1+r)

Isso é a formula do crescimento populacional, so que ai estou colocando a mudança de um tempo para o outro, é primeira formula do crescimento populacional que esta no wikipedia, so que se você quiser ver bastantes tempo, é so fazer isso varias vezes, mais ou menos isso:

 N_{t+1}= N_t\cdot(1+r)\cdot (1+r)\cdot

Mas isso é o mesmo que eu elevar tudo ao quadrado

 N_{t+1}= N_t\cdot(1+r)^2

Então se quiser-mos ver três tempo, pensando em tempo como anos por exemplo é so elevar por 3, quatro ano 4, assim podemos generalizar para:

 N_{t+1}= N_t\cdot(1+r)^t

Pronto, agora é exatamente a formula do wikipedia. Mas vamos voltar para somente uma tempo, t=1, que é:

 N_{t+1}= N_t\cdot(1+r)

Então podemos pensar o seguinte, se tempo é medido de 1 em 1 ano, como dissemos, nesse caso estamos pensando na contribuição de cada individuo uma vez por por ano, como se todo mundo fosse como tartarugas, que põe os ovos uma vez no ano na areia, e nasce todo mundo ao mesmo tempo. Nesse caso seria perfeito. Mas mamíferos, a gente por exemplo, tem filhos ao longo do ano, todo dia esta nascendo gente, sem parar.

Então poderíamos pensar nessa reprodução de 6 em 6 meses, ou seja partir o r em dois momentos. Dois eventos, isso aqui:

 N_{t+1}= N_t\cdot(1+r/2)^2

Mas poderíamos ainda pensar em como seria o crescimento de 4 em 4 meses, ou seja três eventos no ano…

 N_{t+1}= N_t\cdot(1+r/3)^3

Então podemos picotar em 2 eventos, mas podemos também picotar em 3 eventos, mas podemos querer ser mais e mais criteriosos e fazer isso mais em mais até o infinito.
Generalizando isso temos o seguinte:

 N_{t+1}= N_t\cdot(1+r/n)^n

Bem simples a generalização não?

Podemos mexer nessa generalização ainda onde a razão, o tanto que a população cresceu, o N1 dividido pelo N0 fica:

 {N_{t+1}\over N_t}= (1+{r \over n})^n

O que nos leva a questão o que acontece com  (1+{r over n})^n aumentando n mais e mais até o infinito.

 lim_{n\to +\infty}(1+{r \over n})^n

Isso a gente pode tentar resolver a partir de força bruta, que é, a gente vai aumentando o n para algum valor de r e ve o que acontece.
No r a gente faz o seguinte:

n <- 0:100
N0 <- 1
r <- 1
N1 <- N0 * (1 + r/n)^n
 
plot(n, N1/N0, type = "l")
text(50, 2, "For n = 100,")
text(50, 1.6, bquote((1 + frac("r", "n"))^"n" ==.(round(N1[101]/N0, 3))))

E o que vemos é o seguinte:

01

Ou seja, é so pensar um pouco, isso quer dizer que, quanto mais a gente corta um ano em 6 meses, em 4 meses, em 1 mes, em 1 dia, sempre a gente vai ver a população mudar menos e menos. E o limite disso é mais ou menos 2.705, que não é um número qualquer. Esse é o famoso e da matemática, que corramos o que for, cedo ou tarde em qualquer lugar vamos nos deparar com ele, como aqui na biologia :), mas essa conclusão aqui foi do Jacob Bernoulli tentando melhor compreender o crescimento logístico.

Mas então agora vemos que se tentarmos ver as coisas continuamente, temos que:

 lim_{n\to+\infty}(1+{r \over n})^n=e^r

Agora o tempo que tiramos la em cima, podemos devolver ele na equação, assim como o N0 e ficamos com:

 N_{t+1}= N_t\cdot e^{r\cdot t}

Ou seja, podemos chamar o r aqui de taxa de crescimento instantâneo, ou crescimento intrínseco, e estamos vendo o crescimento, infinitamente precisamente quanto as diferenças de tempo.

Isso é um pouco difícil de começar a assimilar, principalmente que não é algo como uma regressão linear, o tamanho da população em algum momento sempre depende da população anterior aquele momento, mas é importante sempre seguir o conselho do senhor John von Neumann, “Young man, in mathematics you don’t understand things. You just get used to them.”.

E por último, aqui tem um antigo post meu perguntando no fórum do R sobre equações diferencias exatamente isso e algumas considerações do Thomas Petzoldt, que é um dos desenvolvedores do pacote deSolve do R


https://stat.ethz.ch/pipermail/r-sig-dynamic-models/2011q1/000063.html

Espero que não tenho cometido nenhuma atrocidade contra a matemática ou biologia aqui.

Referencia: Stevens, M. H. H. A Primer of Ecology with R. UseR! Series, Springer (2009)

Rosalind – Rabbits and Recurrence Relations – FIB

Que emoção, peguei um problema no comecinho, acho que acabou de ser criado.

Antes de começar eu decidi que vou parar de numerar os problemas, acho que numeração é pouco útil, já que problemas vem e vão.

Então começando, aqui a pergunta envolve o uso da sequência de Fibonnaci, o artigo na Wikipédia é bem legal.
O exemplo clássico é o crescimento de coelhos, onde só casais tem filhos, quem esta sozinho não aumenta a população e isso tem que ser contabilizado.

A diferença aqui é que o problema pede para a gente dizer um número n da sequência, e essa sequência ainda terá outra condição que é o número de filhotes que um casal tem não precisa ser necessariamente 1.

01

Então ela tem uma cara de crescimento exponencial, e quanto mais filhotes um casal gera, mais rápido esse crescimento.

Eu acho que a idéia do problema era introduzir o conceito de programação dinâmica pelas explicações, mas para começar eu fiz uma função com o R bem “força bruta”.

fibonacci<-function(a1=1,a2=1,t=10,k=1) {
n<-vector()
n[1]<-a1
n[2]<-a2
for(i in 3:t) {
n[i]<-n[i-1]+k*n[i-2]
}
return(n)
}

Então entramos com a1 e a2, que são os valores iniciais, como descrito no problema, quanto tempo vamos observar (o número de números na serie) e k, que é quantos filhotes cada casal tem.
Esse k na verdade eu não consegui pensar automaticamente onde ele ia, eu multipliquei o primeiro número e depois que fui ver que so funcionava se multiplicar o segundo número da soma, para nunca cair a multiplicação em alguém que não tem par.

Feito isso a solução é simples, só usar a função.

out <- fibonacci(t=5,k=3) format( out[length(out)] , scientific = F) [1] "19"

Importante desligar o formato científico da notação numérica, já que rapidinho o número fica gigante e a resposta tem que ser o número sem notação cientifica, default no R.

O código para o gráfico no R foi esse:

plot(c(fibonacci(t=10,k=1))~c(1:10),xlab="Tempo",ylab="Casais",frame.plot=F,
pch=19,type="b",main="Sequência de Fibonacci")
points(c(fibonacci(t=10,k=2))~c(1:10),col=2,pch=19,type="b")
points(c(fibonacci(t=10,k=3))~c(1:10),col=3,pch=19,type="b")
points(c(fibonacci(t=10,k=4))~c(1:10),col=4,pch=19,type="b")
legend("topleft",col=1:4,pch=19,lty=1,legend=paste("k = ",1:4))

Problema solucionado, eu ainda consegui ser o terceiro a postar minha solução.

Mas como a idéia do problema era o uso de programação dinâmica, eu fui olhar o wikipedia do que se trata o tal do Dynamic programming, até porque em outro problema eu ja não tinha conseguido resolver com o R porque não vectorizei o problema, e o R demorava mais que os 5 minutos para a resposta do problema, e acabou que o python me iluminou, junto com a solução no stack overflow.

Mas a idéia é quebrar o problema em problemas menores, mesmo que similares, mas só computar a solução uma vez. Então em python eu tentei aplicar o primeiro algoritmo geral que esta no artigo do wikipedia sobre Dynamic programming.

def fib_rek(n):
    if n == 0:
        return 0
    if n == 1:
        return 1
    else:
        return fib_rek(n-1) + fib_rek(n-2)

É muito louco pensar que a gente ta fazendo uma função que reutiliza ela, que reutiliza ela até o mais simples possível.

Então essa chamada de função ta re-chamando ela mesmo, que como descrito no wikipedia faz uma árvore assim para n=5

1 fib(5) 2 fib(4) + fib(3) 3 (fib(3) + fib(2)) + (fib(2) + fib(1)) 4 ((fib(2) + fib(1)) + (fib(1) + fib(0))) + ((fib(1) + fib(0)) + fib(1)) 5 (((fib(1) + fib(0)) + fib(1)) + (fib(1) + fib(0))) + ((fib(1) + fib(0)) + fib(1))

E após apanha um pouco, eu consegui adicionar o k, que é o número de filhotes para o casal, que aqui a gente tem que usar ele e passar o k para frente na função denovo, esse “passar para frente denovo” custou um tempo para assimilar.

def fib_rek(n,k):
    if n == 0:
        return 0
    if n == 1:
        return 1
    else:
        return fib_rek(n-1,k) + k*fib_rek(n-2,k)

E por último, a outra solução que consiste em memorizar o resultado uma vez computado para não fazer a mesma conta 2 vezes. Como aqui eu ja não estava mais dando conta eu peguei uma solução na internet para olhar.

# dynamic programming
def fib_mem(n):
    m = []
    for i in range(n+1):
        m.append(-1)
    m[0] = 0
    m[1] = 1
    def fib_mem_rek(n):
        if m[n] != -1:
            return m[n]
        else:
            val = fib_mem_rek(n-1) + fib_mem_rek(n-2)
            m[n] = val
            return val
    return fib_mem_rek(n)

Então criamos uma lista vazia, enchemos de -1 (acho que o caso é encher de algo para comparar depois), preenchemos o começo da lista, e ai rodamos dentro outra função que, se o valor é diferente de -1, a gente deixa ele ja que ja esta solucionado, senão a gente soluciona. Eu acho que é mais ou menos isso, eu achei bem complicado. Mas uma vez recebida a função de mão beijada assim não é dificil adicionar a variação da produção de filhotes.

def fib_mem(n,k):
    m = []
    for i in range(n+1):
        m.append(-1)
    m[0] = 0
    m[1] = 1
    def fib_mem_rek(n,k):
        if m[n] != -1:
            return m[n]
        else:
            val = fib_mem_rek(n-1,k) + k*fib_mem_rek(n-2,k)
            m[n] = val
            return val
    return fib_mem_rek(n,k)

E para fechar, se alguém chegou nesse parágrafo, o coursera abriu uma disciplina de Bioinformática, que o conhecimento básico requisito é, tharararam, os primeiros problemas do Rosalind, e mais foda que isso, ela é dada pelos caras que fundaram o Rosalind. Não sei se vai ser difícil, mas pode ser uma oportunidade legal, eu estarei por la apanhando para a bioinformática :).

Crescimento discreto vs crescimento contínuo

Vamos olhar mais um pouco o crescimento logístico.

Podemos fazer uma função no R para descrever ele assim:

logistic <- function(alpha=0.01, rd=1, N0=2, t=10) {
N <- c(N0, numeric(t))
for(i in 1:t) N[i+1] <- {N[i] + rd * N[i] * (1 – alpha * N[i])}
return(N) }

Olha o que essa função vai fazer, vai pegar os 4 parâmetros, capacidade suporte (alpha, K é o mesmo que 1/alpha aqui, depois veremos porque fazer isso), taxa de crescimento, o estado inicial da população, e quanto temos queremos acompanhar ela.
Para cada tempo o processo é simples, pegamos a população anterior, nesse caso no tempo 0 é 2 indivíduos e aplicamos o que vimos em crescimento logístico.

Assim obtemos:

01

Mas aqui tratamos o tempo como algo contínuo, por exemplo olhando a população de 1 em 1 ano, mas ai podemos pensar, o que acontece em 6 meses, quer dizer todos os filhos não nascem de uma vez num dia do ano?

Muitas populações não tem seu crescimento ocorrendo de forma assim discreta, e nem sempre temos medidas perfeitamente simétricas(exatos 365 dias de diferença entre medidas), o que torna essa abordagem complicada. Mas então fica a dúvida, o que acontece no intervalo entre um tempo e outro?

Poderíamos pensar que basta passar uma linha em cada ponto e teríamos uma estimativa entre esses pontos assim.

02

Mas outra abordagem é pensar continuamente, descrevendo o crescimento populacional como uma equação diferencial ordinária.

A única mudança é que ao invés de falar na mudança da população de ano baseado no ano anterior, vemos como a população muda ao longo do tempo, dN/dt.

clogistic <- function(times, y, parms) {
n <- y[1]
r <- parms[1]
alpha <- parms[2]
dN.dt <- r * n * (1 – alpha * n)
return( list( c(dN.dt) ) )
}

Podemos usar a função ode do pacote desolve, que soluciona equações diferenciais ordinárias. E temos o seguinte:

03

A solução é igual, só não tão retinha quanto simplesmente ligar os pontos, porque a cada momento alguém esta nascendo, a população esta mudando continuamente, ou seja aqui não estamos dando passos discretos, mas passos contínuos.

Bem vamos deixar aqui essas duas funções guardadas para olhar mais alguns pontos depois.

Referencia: Stevens, M. H. H. A Primer of Ecology with R. UseR! Series, Springer (2009)

logistic <- function(alpha=0.01, rd=1, N0=2, t=10) {
N <- c(N0, numeric(t))
for(i in 1:t) N[i+1] <- {N[i] + rd * N[i] * (1 – alpha * N[i])}
return(N) }
 
pop<-logistic()
 
plot(pop~c(0:10),xlab="Tempo",ylab="População",frame.plot=F,pch=19)
legend("topleft",legend="Crescimento Discreto",pch=19,bty="n")
 
plot(pop~c(0:10),xlab="Tempo",ylab="População",frame.plot=F,pch=19,type="b")
legend("topleft",legend="Crescimento Discreto",pch=19,lty=1,bty="n")
 
clogistic <- function(times, y, parms) {
n <- y[1]
r <- parms[1]
alpha <- parms[2]
dN.dt <- r * n * (1 – alpha * n)
return( list( c(dN.dt) ) )
}
 
prms <- c(r=1, alpha=.01)
init.N <- 1
t.s <- seq(0.1, 10, by=0.1)
 
library(deSolve)
out <- ode(y=init.N, times=t.s, clogistic, parms=prms)
 
plot(pop~c(0:10),xlab="Tempo",ylab="População",frame.plot=F,pch=19,type="b")
points(out[,1], out[,2], type=’l',col="red",lty=2)
legend("topleft",legend=c("Crescimento Discreto","Crescimento Continuo"),
pch=c(19,NA),lty=c(1,2),col=c("black","red"),bty="n")

Aproximação por grid em estatística bayesiana.

A gente já viu nesse post como usar estatística bayesiana sem o mcmc. Onde definimos um prior usando uma distribuição (a beta), calculávamos o likelihood e aplicávamos o teorema de bayes.

Acontece que existe outra forma ainda de fazer as coisas, naquele exemplo trabalhamos com uma distribuição continua, ou seja que que acomodava qualquer valor, seja 0.5 ou 0.55. Mas podemos definir as coisas sem usar uma distribuição, definindo um “grid” de possibilidades.

Ou seja eu defino um grid que represente todas as possibilidades, mas eu posso fechar ele de forma a não existir o 0.55, ou a chance de colonizar um arbusto é 50% ou é 60%, sem meio termo como na distribuição beta.

Isso é legal que agora o divisor do teorema de bayes vai ser uma soma, a soma de todas as possibilidades, e não uma integral como era no outro caso.

Então definimos nossa probabilidades do prior:

01

Uma coisa, a gente usa essas barras e não um histograma para ser claro na precisão. Ou seja eu to falando de 10 em 10%, eu poderia ser mais preciso, mas ai seriam mais barras e se quisesse ser infinitamente preciso e cobrir todos os valores possíveis eu usaria a distribuição beta horas.

Outra coisa legal é que agora vemos como é possível estabelecer um prior que não respeita as forma possíveis da distribuição beta como vimos anteriormente.
Nesse caso podemos pensar no prior como duas vertentes, uma que acha que a ocupação é baixa, está em torno dos 20% e outra que é alta, esta em torno dos 60%.

02

E vamos la, coletamos nossos dados, olhamos vários arbustos novamente, calculamos o likelehood, que de 10 arbustos, 7 tinham aranhas e começamos a achar que realmente é algo em torno de 60% a ocupação desses arbustos. Ou seja a parte do prior que estava ocupando uma baixa taxa de ocupação é não la muito verossímil com o que observamos, afinal quase todos os arbustos estão ocupado, não há como a ocupação não ser alta.

Bem acho que assim como o método exato falado no post anterior, a descrição das probabilidades por grid não vai ser muito útil, mas é legal para fazer alguns experimentos e entender melhor a aplicação do teorema de bayes, ainda mais que nem integral precisamos nesse caso, além criar alguns insights de qual o efeito da precisão usada no grid (usar poucos pontos ou muitos para descrever o mundo). Criar prior esquisitos e ver o que acontece ou ainda tentar com quantidades diferentes de dados para entender melhor o que vai acontecer. Daqui é so mudar o script, alterando a precisão do grid e mudar a coleta de dados (a forma como se gera ela com a função rbinom).

####################################################################
#Aproximação por grid
#Exemplo Baseado no livro Doing Bayesian Data Analysis
#http://www.indiana.edu/~kruschke/DoingBayesianDataAnalysis/
#################################################################
#Prior
Theta<-seq(0.01,0.99,0.05)
pTheta<-c(0,0.025,0.05,0.075,0.2,0.075,0.05,0.025,0,
0.025,0.05,0.075,0.2,0.075,0.05,0.025,0,0,0,0)
sum(pTheta)
 
#grafico do prior
plot(pTheta,type="n",frame.plot=F,xaxt="n",xlim=c(0,21))
segments(1:length(pTheta),0,1:length(pTheta),pTheta,lwd=3)
axis(1,at=1:20,labels =Theta,cex.axis=0.5)
 
#Dados
set.seed(123)
Data<-rbinom(10,1,0.8)
 
#Como fazer os gráficos
credib<-0.95
nToPlot<-length(Theta)
 
z = sum( Data==1 ) # Número de arbustos Ocupados
N = length( Data ) # Número de arbustos Observados
 
# likelihood dos dados para cada valor de Theta que usamos
pDataGivenTheta = Theta^z * (1-Theta)^(N-z)
 
# Computando a evidencia para a distribuição posterios
pData = sum( pDataGivenTheta * pTheta )
 
#Multiplicando o Prior vezes o Likelihood dividido pela Evidencia
#(Teorema Bayes)
pThetaGivenData = pDataGivenTheta * pTheta / pData
 
# Gráfico dos Resultados
layout( matrix( c( 1,2,3 ) ,nrow=3 ,ncol=1 ,byrow=FALSE ) ) # Painel 3×1
par( mar=c(3,3,1,0) , mgp=c(2,1,0) , mai=c(0.5,0.5,0.3,0.1) ) # Margens
dotsize = 1 #Tamanho das bolinhas
 
#Se escolhermos 1 milhão de valores para Theta, num vai dar certo o plot
#pois vão ser muitos pontos, mas da pra fazer um plot menor com media entre
#thetas, como definimos la em cima quantos ponto plotar
 
nteeth = length(Theta)
if ( nteeth > nToPlot ) {
thinIdx = seq( 1, nteeth , round( nteeth / nToPlot ) )
if ( length(thinIdx) < length(Theta) ) {
thinIdx = c( thinIdx , nteeth ) # makes sure last tooth is included
}
} else { thinIdx = 1:nteeth }
 
# Grafico do prior.
 
meanTheta = sum( Theta * pTheta ) # mean of prior, for plotting
plot( Theta[thinIdx] , pTheta[thinIdx] , type="p" , pch=19 , cex=dotsize ,
xlim=c(0,1) , ylim=c(0,1.1*max(pThetaGivenData)) , cex.axis=1.2 ,
xlab=bquote(theta) , ylab=bquote(p(theta)) , cex.lab=1.5 ,
main="Prior" , cex.main=1.5 )
if ( meanTheta > .5 ) {
textx = 0 ; textadj = c(0,1)
} else {
textx = 1 ; textadj = c(1,1)
}
text( textx , 1.0*max(pThetaGivenData) ,
bquote( "média(" * theta * ")=" * .(signif(meanTheta,3)) ) ,
cex=2.0 , adj=textadj )
# Grafico do likelihood: p(Data|Theta)
 
plot(Theta[thinIdx],pDataGivenTheta[thinIdx],type="p",pch=19,cex=dotsize
,xlim=c(0,1) ,cex.axis=1.2 ,xlab=bquote(theta)
,ylim=c(0,1.1*max(pDataGivenTheta))
,ylab=bquote( "p(D|" * theta * ")" )
,cex.lab=1.5 ,main="Likelihood" ,cex.main=1.5 )
if ( z > .5*N ) {
textx = 0 ; textadj = c(0,1)
} else {
textx = 1 ; textadj = c(1,1)
}
 
text(textx,1.0*max(pDataGivenTheta),cex=2.0,
bquote( "Data: z=" * .(z) * ",N=" * .(N) ) ,adj=textadj )
 
# Grafico do posterior.
meanThetaGivenData = sum( Theta * pThetaGivenData )
plot(Theta[thinIdx],pThetaGivenData[thinIdx],type="p" ,pch=19,cex=dotsize
,xlim=c(0,1) ,ylim=c(0,1.1*max(pThetaGivenData)) ,cex.axis=1.2
,xlab=bquote(theta) ,ylab=bquote( "p(" * theta * "|D)" )
,cex.lab=1.5 ,main="Posterior" ,cex.main=1.5 )
 
if ( meanThetaGivenData > .5 ) {
textx = 0 ; textadj = c(0,1)
} else {
textx = 1 ; textadj = c(1,1)
}
 
text(textx ,1.00*max(pThetaGivenData),cex=2.0,
bquote("média(" * theta * "|D)=" * .(signif(meanThetaGivenData,3)) ),adj=textadj )
text(textx,0.75*max(pThetaGivenData),cex=2.0,
bquote( "p(D)=" * .(signif(pData,3)) ) ,adj=textadj )
 
# Marcando o intervalo de maxima densidade, decidimos o intervalo la em cima.
#função para calcular o HDI
HDIofGrid = function( probMassVec , credMass=0.95 ) {
sortedProbMass = sort( probMassVec , decreasing=T )
HDIheightIdx = min( which( cumsum( sortedProbMass ) >= credMass ) )
HDIheight = sortedProbMass[ HDIheightIdx ]
HDImass = sum( probMassVec[ probMassVec >= HDIheight ] )
return( list( indices = which( probMassVec >= HDIheight ) ,
mass = HDImass , height = HDIheight ) )
}
 
HDIinfo = HDIofGrid(pThetaGivenData,credMass=credib)
points(Theta[ HDIinfo$indices ],rep(HDIinfo$height,length( HDIinfo$indices ))
,pch=19,cex=0.5,col="red" )
text(mean(Theta[HDIinfo$indices]), HDIinfo$height,
bquote( .(100*signif(HDIinfo$mass,3)) * "% HDI" ),
adj=c(0.5,-1.5),cex=1.5 )
 
#Marcando uma linha no Intervado de 95%
lowLim = Theta[min(HDIinfo$indices)]
highLim = Theta[max(HDIinfo$indices)]
lines(c(lowLim,lowLim),c(-0.5,HDIinfo$height),type="l",lty=2,lwd=1.5)
lines(c(highLim,highLim),c(-0.5,HDIinfo$height),type="l",lty=2,lwd=1.5)
segments(lowLim,HDIinfo$height,highLim,HDIinfo$height,lty=2,lwd=1.5)
text(lowLim,HDIinfo$height,bquote(.(round(lowLim,3))),
adj=c(1.1,-0.1),cex=1.2 )
text(highLim,HDIinfo$height,bquote(.(round(highLim,3))),
adj=c(-0.1,-0.1),cex=1.2 )

Rosalind problema 12 – Enumerating Gene Orders – PERM

Permutações de 3 bolas

Esse é outro problema que tinha tudo para ser fácil, mas não foi.
A ideia é bem simples, dado um número n, no exemplo 3, temos que dizer quantas permutações são possíveis fazer e quais são elas.

Nesse caso 3 quer dizer que vamos permutar 1,2 e 3. So lembrar então que se vamos combinar 3 coisas, então o primeiro item tem 3 possibilidades, o segundo item só tem 2 possibilidades e o último tem uma, logo o total de possibilidades é 6 nesse caso. Como explicado no Wikipédia

Essa parte é bem do segundo grau, fatorial e aulas afim.
Agora a complicação começa em como fazer um gerador para apresentar todas a permutações possíveis. Nesse caso seria:

6 1 2 3 1 3 2 2 1 3 2 3 1 3 1 2 3 2 1

Então temos o 6 que é o número de permutações e quais são elas logo abaixo.

No R minha primeira idéia foi a seguinte

num<-3
n<-1:num
prod(n)
 
for(i in 1:3) {
  for(j in 1:3) {
    for(k in 1:3) {
      if(i != j &amp; j != k &amp; i != k) { print(c(n[i],n[j],n[k])) }
    }
  }
}

Então eu fiz um loop aninhado, precisando de quantos for fossem a quantidade de itens da lista, e imprimia o número sempre que os índices i j e k não fossem iguais.
Mas isso se demonstrou péssimo, porque eu não consegui imaginar um jeito de fazer um função que construísse uma função assim para mim.

Para resolver logo o problema, eu testei usar uma livraria de python, que vem com ele, a itertools. Com ela foi simples, ele tem uma função chamada permutations que faz todas as permutações

import itertools
n=3
objetos=range(1,n+1)
permuta=list(itertools.permutations(objetos))
a=len(permuta)
 
print a
for i in range(a):
    print permuta[i]

Feito isso foi fácil. obtivemos:

>>> 6 (1, 2, 3) (1, 3, 2) (2, 1, 3) (2, 3, 1) (3, 1, 2) (3, 2, 1)

Mas no final das contas eu não aprendi nada sobre como fazer as permutações, mas olhando as soluções da para aprender bastante.

Então que eu vi uma galera que usa python usando a função yield. Foi difícil entender como isso funciona, até que eu vi esse post em outro blog que foi bem instrutivo.

Ai deu para entender o que isso faz, da para fazer uma função assim:

def gerador():
    for i in range(1,n+1):
        yield i
 
teste=gerador()
 
for i in teste:
    print i

Então conforme o i vai rodando, o yield vai soltando um item, dai a gente atribui a um nome essa função, e com um loop como no acima teremos:

>>> 1 2 3

Ou seja a gente cria a lista como um objeto iterável, o que não sabia fazer no R, vale a pena ler o post que citei acima que explica muito bem isso. Mas tendo isso em mente da para entender mais ou menos o primeiro algoritmo proposto no wikipedia, apresentando as permutações na “lexicographic order”, não sei escrever isso em português.

n=3
 
def factorial(x):
    if x == 0:
        return 1
    else:
        return x * factorial(x - 1)
 
def permutations(items):
    def p(items, n):
        if n == 0:
            yield []
        else:
            for i in range(n):
                for j in p(items[:i] + items[i+1:], n - 1):
                    yield [items[i]] + j
    return p(items, len(items))
 
print factorial(n)
 
lst = permutations(range(1,n+1))
 
for i in lst:
    print &quot; &quot;.join(map(str,i))

E temos:

>>>6 1 2 3 1 3 2 2 1 3 2 3 1 3 1 2 3 2 1

Aqui a primeira função é só para calcular o fatorial, a segunda é usando o yield, com dois indexadores, o i e o j apresentando todas as iterações, mas o gerador por si so não imprime nada, então imprimimos com um loop. Particularmente essa foi uma função que copiei das soluções mas que foi a que melhor deu para entender a idéia, mais algum trabalho nela e é possível tirar aquela complicação do último comando print, mas basicamente a saída são varias listas, então ele tira o colchetes, virgulas e tal para a saída ser no estilo do proposto pelo problema, mas como desde o carnaval não posto nada, decidi fazer algo logo 🙂