Programação para as
Ciências Experimentais
2008/9
Teórica 12
Ludwig Krippahl, 2009
Na aula de hoje...



Minimização multidimensional
Exemplo: estimativa do efeito de um
antibiótico no crescimento bacteriano.
Trabalho Prático 2
Ludwig Krippahl, 2009
2
Crescimento bacteriano

Equação de Verhulst:
dB/dt = cB – mB2
Ludwig Krippahl, 2009
3
Crescimento bacteriano

Equação de Verhulst:
dB/dt = cB – mB2
Variação do número de organismos
Ludwig Krippahl, 2009
4
Crescimento bacteriano

Equação de Verhulst:
dB/dt = cB – mB2
É o ritmo de crescimento vezes o número
de organismos
Ludwig Krippahl, 2009
5
Crescimento bacteriano

Equação de Verhulst:
dB/dt = cB – mB2
Menos a taxa de mortalidade vezes o
quadrado desse número. A mortalidade
resulta da competição por recursos.
Ludwig Krippahl, 2009
6
Crescimento bacteriano
Ludwig Krippahl, 2009
7
Crescimento bacteriano

Problema:
• Dado um conjunto de medições, ajustar os
parâmetros da equação
Ludwig Krippahl, 2009
8
Crescimento bacteriano

Problema:
• Dado um conjunto de medições, ajustar os
•
parâmetros da equação
Mas são dois parâmetros: crescimento e
mortalidade. Precisamos de uma minimização
a duas dimensões.
Ludwig Krippahl, 2009
9
Minimização multidimensional

Método:
• Minimizar uma variável de cada vez até
chegar a um ponto fixo, a menos da precisão
desejada
Ludwig Krippahl, 2009
10
Minimização multidimensional

Método :
• Partir de um ponto inicial, um valor para cada
•
•
•
•
variável.
Encontrar o mínimo de uma variável.
Alterar o vector das variáveis
Encontrar o mínimo da próxima.
Repetir para todas, as vezes que for
necessário.
Ludwig Krippahl, 2009
11
Minimização multidimensional
Ponto inicial
Ludwig Krippahl, 2009
12
Minimização multidimensional
Ponto inicial
Mínimo de X
Ludwig Krippahl, 2009
13
Minimização multidimensional
Ponto inicial
Mínimo de X
Mínimo de Y
Ludwig Krippahl, 2009
14
Minimização multidimensional
Ponto inicial
Mínimo de X
Novo mínimo de X
Mínimo de Y
Ludwig Krippahl, 2009
15
Minimização multidimensional

Vamos modificar a minfn
• Para partir do ponto dado e não ser preciso
•
especificar os três pontos iniciais (é mais
eficiente começar com 3 pontos juntos
quando próximo do mínimo)
Para procurar o mínimo de uma de várias
variáveis.
Ludwig Krippahl, 2009
16
Os 3 pontos iniciais
X1 é o ponto
dado
Ludwig Krippahl, 2009
17
Os 3 pontos iniciais
X1
Xm próximo de X1
Ludwig Krippahl, 2009
18
Os 3 pontos iniciais
X1
Desce?
Xm próximo de X1
Se não, troca
Ludwig Krippahl, 2009
19
Os 3 pontos iniciais
X1
Xm
X2 a 1.618*(Xm-X1)
Ludwig Krippahl, 2009
20
Os 3 pontos iniciais
Y2>Ym?
Não, continua:
X1=Xm
X1
Xm
X2
Xm=X2
Ludwig Krippahl, 2009
21
Os 3 pontos iniciais
Y2>Ym?
Não, continua:
X1=Xm
X1
Xm
Xm=X2
Ludwig Krippahl, 2009
22
Os 3 pontos iniciais
Y2>Ym?
Não, continua:
X1=Xm
X1
Xm
X2
Xm=X2
Ludwig Krippahl, 2009
23
Os 3 pontos iniciais
Y2>Ym?
Sim.Devolve:
X1 Xm X2 Ym
X2
X1
Xm
(Ym para começar a minimização)
Ludwig Krippahl, 2009
24
Os 3 pontos iniciais
function
[x1,xm,x2,ym]=
mininicial(funcao,params,vars,indice,delta)
Ludwig Krippahl, 2009
25
Os 3 pontos iniciais
function
[x1,xm,x2,ym]=
mininicial(funcao,params,vars,indice,delta)
Os valores a devolver, os 3 pontos de x e o y
do meio que precisamos para começar a
minimização.
Ludwig Krippahl, 2009
26
Os 3 pontos iniciais
function
[x1,xm,x2,ym]=
mininicial(funcao,params,vars,indice,delta)
Nome da função a minimizar e os parâmetros
constantes que precisamos para a avaliar.
E.g.: os coeficientes do polinómio, os dados
experimentais, etc
Ludwig Krippahl, 2009
27
Os 3 pontos iniciais
function
[x1,xm,x2,ym]=
mininicial(funcao,params,vars,indice,delta)
Um vector com os valores das N variáveis da
função (a função tem várias dimensões)
Ludwig Krippahl, 2009
28
Os 3 pontos iniciais
function
[x1,xm,x2,ym]=
mininicial(funcao,params,vars,indice,delta)
O índice no vector da variável que estamos a
minimizar agora (a função tem várias, mas só
conseguimos lidar com uma de cada vez)
Ludwig Krippahl, 2009
29
Os 3 pontos iniciais
function
[x1,xm,x2,ym]=
mininicial(funcao,params,vars,indice,delta)
Tamanho do passo inicial para calcular o
primeiro Xm a partir do X inicial, que será o
valor que vem em vars(indice).
Ludwig Krippahl, 2009
30
Os 3 pontos iniciais
razaod=1.618;
Razão dourada
x1=vars(indice);
y1=feval(funcao,params,vars);
xm=x1+delta;
vars(indice)=xm;
ym=feval(funcao,params,vars);
Ludwig Krippahl, 2009
31
Os 3 pontos iniciais
razaod=1.618;
Calcular x1 e y1.
x1=vars(indice);
y1=feval(funcao,params,vars);
Nota: a função cujo nome
xm=x1+delta;
foi dado em funcao precisa
vars(indice)=xm;
de todo o vars, mas só
ym=feval(funcao,params,vars);
alteramos o elemento
indice
Ludwig Krippahl, 2009
32
Os 3 pontos iniciais
razaod=1.618;
x1=vars(indice);
y1=feval(funcao,params,vars);
xm=x1+delta;
Mesma coisa para xm e ym
vars(indice)=xm;
ym=feval(funcao,params,vars);
Ludwig Krippahl, 2009
33
Os 3 pontos iniciais
if ym>y1
Se for “a subir” troca o x1
t=ym;
com o xm, e o y1 com o ym.
ym=y1;
y1=t;
t=xm;
xm=x1;
x1=t;
endif
x2=xm+razaod*(xm-x1);
vars(indice)=x2;
y2=feval(funcao,params,vars);
Ludwig Krippahl, 2009
34
Os 3 pontos iniciais
if ym>y1
Calcula o x2 a uma distância
t=ym;
de xm igual a 1.618.. vezes o
ym=y1;
y1=t;
intervalo (xm-x1).
t=xm;
xm=x1;
x1=t;
endif
x2=xm+razaod*(xm-x1);
vars(indice)=x2;
y2=feval(funcao,params,vars);
Ludwig Krippahl, 2009
35
Os 3 pontos iniciais
if ym>y1
Se x1>xm, continua para
t=ym;
valores maiores.
ym=y1;
y1=t;
Se foi trocado, x1-xm é
t=xm;
xm=x1;
negativo e segue para
x1=t;
valores menores.
endif
x2=xm+razaod*(xm-x1);
vars(indice)=x2;
y2=feval(funcao,params,vars);
Ludwig Krippahl, 2009
36
Os 3 pontos iniciais
Enquanto continua “a
descer”, avança com os
pontos,
while y2<ym
x1=xm;
xm=x2;
ym=y2;
x2=xm+razaod*(xm-x1);
vars(indice)=x2;
y2=feval(funcao,params,vars);
endwhile
Ludwig Krippahl, 2009
37
Minimização
Antes era:
function
xm=minfn(func,params,x1,xm,x2,prec)
 Agora é:
function
xm=minfnvec(func,params,vars,indice,prec)

Ludwig Krippahl, 2009
38
Minimização multidimensional
function
xm=minfnvec(func,params,vars,indice,prec)
Valor da variável que está a minimizar no
mínimo da função considerando apenas
esta variável.
Ludwig Krippahl, 2009
39
Minimização multidimensional
function
xm=minfnvec(func,params,vars,indice,prec)
Função (o nome, em string) e parâmetros
constantes, como de costume.
Ludwig Krippahl, 2009
40
Minimização multidimensional
function
xm=minfnvec(func,params,vars,indice,prec)
Vector com os valores das várias variáveis
no ponto inicial, de onde parte à procura
do mínimo.
Ludwig Krippahl, 2009
41
Minimização multidimensional
function
xm=minfnvec(func,params,vars,indice,prec)
Índice da variável onde procurar o mínimo.
Ludwig Krippahl, 2009
42
Minimização multidimensional
function
xm=minfnvec(func,params,vars,indice,prec)
Precisão (tamanho do intervalo abaixo do
qual consideramos ter encontrado o
mínimo)
Ludwig Krippahl, 2009
43
Minimização multidimensional
Esta função é igual à minfn, excepto:
• Usa o mininicial para determinar os 3 pontos e
o valor de ym
[x1,xm,x2,ym]=mininicial(func,params,vars,indice,prec);
Nota: um bom valor para o delta é a precisão:
começamos do intervalo mais pequeno.
Ludwig Krippahl, 2009
44
Minimização multidimensional
Esta função é igual à minfn, excepto:
• Usa o mininicial para determinar os 3 pontos e
•
o valor de ym
Tem de atribuir o valor correcto a vars(indice)
antes de chamar a função fornecida em func
xn=c1*xm+c2*x1;
vars(indice)=xn;
yn=feval(func,params,vars);
Ludwig Krippahl, 2009
45
Minimização multidimensional
Esta função é igual à minfn, excepto:
• Usa o mininicial para determinar os 3 pontos e
•

o valor de ym
Tem que atribuir o valor correcto a vars(indice)
antes de chamar a função fornecida em func
Mas minfnvec ainda só minimiza numa
dimensão (a dimensão indicada em
indice).
Ludwig Krippahl, 2009
46
Minimização multidimensional
Precisamos de:
function
xs=multimin(funcao,params,xs,prec)

Ludwig Krippahl, 2009
47
Minimização multidimensional
function
xs=multimin(funcao,params,xs,prec)
Vector com os valores de todas as variáveis
no mínimo da função
Ludwig Krippahl, 2009
48
Minimização multidimensional
function
xs=multimin(funcao,params,xs,prec)
Nome da função a minimizar.
Ludwig Krippahl, 2009
49
Minimização multidimensional
function
xs=multimin(funcao,params,xs,prec)
Parâmetros constantes...
Ludwig Krippahl, 2009
50
Minimização multidimensional
function
xs=multimin(funcao,params,xs,prec)
Ponto inicial (valores de todas as variáveis
da função de onde partir à procura do
mínimo)
Ludwig Krippahl, 2009
51
Minimização multidimensional
function
xs=multimin(funcao,params,xs,prec)
Precisão
Ludwig Krippahl, 2009
52
Minimização multidimensional
while true
xvs=xs;
for f=1:length(xs)
xs(f)=minfnvec(funcao,params,xs,f,prec);
endfor
disp("Valores até agora:")
disp(xs);
fflush(stdout)
if sum(abs(xs-xvs))<prec
break
endif
Ciclo
só
termina
no
endwhile
Ludwig Krippahl, 2009
break.
53
Minimização multidimensional
while true
Guarda os valores antigos
xvs=xs;
for f=1:length(xs) dos xs
xs(f)=minfnvec(funcao,params,xs,f,prec);
endfor
disp("Valores até agora:")
disp(xs);
fflush(stdout)
if sum(abs(xs-xvs))<prec
break
endif
endwhile
Ludwig Krippahl, 2009
54
Minimização multidimensional
while true
Minimiza em cada dimensão,
xvs=xs;
for f=1:length(xs) actualizando o valor nos xs
xs(f)=minfnvec(funcao,params,xs,f,prec);
endfor
disp("Valores até agora:")
disp(xs);
fflush(stdout)
if sum(abs(xs-xvs))<prec
break
endif
endwhile
Ludwig Krippahl, 2009
55
Minimização multidimensional
while true
Mostra o progresso do
xvs=xs;
cálculo indicando os
for f=1:length(xs)
valores correntes.
xs(f)=minfnvec(funcao,params,xs,f,prec);
endfor
disp("Valores até agora:")
disp(xs);
fflush(stdout);
if sum(abs(xs-xvs))<prec
break
endif
endwhile
Ludwig Krippahl, 2009
56
Minimização multidimensional
while true
xvs=xs;
for f=1:length(xs)
xs(f)=minfnvec(funcao,params,xs,f,prec);
Se o total da variação
endfor
absoluta das variáveis é
disp("Valores até agora:")
inferior à precisão, acabou.
disp(xs);
fflush(stdout)
if sum(abs(xs-xvs))<prec
break
endif
endwhile
Ludwig Krippahl, 2009
57
Crescimento bacteriano

De volta ao problema:
dB/dt = cB – mB2
• Integramos pelo método de Euler, com a
função:
function
mat=crescimento(cresc,mort,dt,qini,tfinal)
Ludwig Krippahl, 2009
58
Crescimento bacteriano
function
mat=crescimento(cresc,mort,dt,qini,tfinal)
Matriz com os valores de tempo e número de
bactérias em duas colunas
Ludwig Krippahl, 2009
59
Crescimento bacteriano
function
mat=crescimento(cresc,mort,dt,qini,tfinal)
Taxa de crescimento
Ludwig Krippahl, 2009
60
Crescimento bacteriano
function
mat=crescimento(cresc,mort,dt,qini,tfinal)
Taxa de mortalidade
Ludwig Krippahl, 2009
61
Crescimento bacteriano
function
mat=crescimento(cresc,mort,dt,qini,tfinal)
Passo de integração
Ludwig Krippahl, 2009
62
Crescimento bacteriano
function
mat=crescimento(cresc,mort,dt,qini,tfinal)
Quantidade inicial de organismos.
Ludwig Krippahl, 2009
63
Crescimento bacteriano
function
mat=crescimento(cresc,mort,dt,qini,tfinal)
Tempo final.
Ludwig Krippahl, 2009
64
Crescimento bacteriano
mat=[0,qini];
B=qini;
for t=dt:dt:tfinal
dB=B*cresc-B^2*mort;
B=B+dB*dt;
mat=[mat;t,B];
endfor
Ludwig Krippahl, 2009
65
Crescimento bacteriano

Agora precisamos de calcular o erro do
modelo aos dados experimentais.
Análogo ao que fizemos para as
reacções químicas, mas desta vez com
duas variáveis.
Ludwig Krippahl, 2009
66
Crescimento bacteriano
function err=errocres(dados,vars)
mat=crescimento(vars(1),vars(2),10,0.1,400);
y=interpol(mat,dados(:,1));
err=sum((y-dados(:,2)).^2);
endfunction
Matriz com a simulação, 400 minutos, passo
de 10 minutos.
Ludwig Krippahl, 2009
67
Crescimento bacteriano
function err=errocres(dados,vars)
mat=crescimento(vars(1),vars(2),10,0.1,400);
y=interpol(mat,dados(:,1));
err=sum((y-dados(:,2)).^2);
endfunction
Nota: Quantidade em “kilobactérias”.
Explicação adiante...
Ludwig Krippahl, 2009
68
Crescimento bacteriano
function err=errocres(dados,vars)
mat=crescimento(vars(1),vars(2),10,0.1,400);
y=interpol(mat,dados(:,1));
err=sum((y-dados(:,2)).^2);
endfunction
Interpolar os valores simulados para os pontos
dos dados experimentais.
Ludwig Krippahl, 2009
69
Crescimento bacteriano
function err=errocres(dados,vars)
mat=crescimento(vars(1),vars(2),10,0.1,400);
y=interpol(mat,dados(:,1));
err=sum((y-dados(:,2)).^2);
endfunction
Erro quadrático....
Ludwig Krippahl, 2009
70
Crescimento bacteriano
function err=errocres(dados,vars)
mat=crescimento(vars(1),vars(2),10,0.1,400);
y=interpol(mat,dados(:,1));
err=sum((y-dados(:,2)).^2);
endfunction
Para o erro quadrático médio usar mean ou
dividir pelo total.
Ludwig Krippahl, 2009
71
Crescimento bacteriano
Para testar, simulamos dados com estes
parâmetros: 10 pontos de 30 em 30 minutos.
Nota: cada linha da matriz são 10 minutos.
vals=[0.040234,0.001877];
mat=crescimento(vals(1),vals(2),10,0.1,400);
dados=mat(3:3:30,:)
Ludwig Krippahl, 2009
72
Crescimento bacteriano
Escolhemos um ponto inicial diferente, e
minimizamos:
xs=multimin("errocres",dados,[0.05,0.005],1e-4)
Ludwig Krippahl, 2009
73
Crescimento bacteriano
Simulamos com os parâmetros calculados e
comparamos:
mat2=crescimento(xs(1),xs(2),10,0.1,400);
hold off
plot(dados(:,1),dados(:,2),"or");
hold on
plot(mat2(:,1),mat2(:,2));
Ludwig Krippahl, 2009
74
Crescimento bacteriano
Ludwig Krippahl, 2009
75
Crescimento bacteriano

Nota sobre as “kilobactérias”:
•
•
Com esta equação, se contarmos em unidades de uma
bactéria o parâmetro da mortalidade tem de ser mil
vezes mais pequeno.
Em geral, é melhor escolher as unidades de forma a
que a função tenha uma escala semelhante nas várias
dimensões. Desta forma a taxa de crescimento e de
mortalidade têm apenas uma ordem de grandeza de
diferença em vez de quatro.
Ludwig Krippahl, 2009
76
Processar dados experimentais



Queremos estudar o efeito da meticilina no
crescimento de uma bactéria.
Duas pessoas, Ana e Carlos, cresceram lotes
da bactéria em meios com e sem meticilina, e
contaram as colónias de amostras retiradas de
30 em 30 minutos.
A concentração inicial era de 100 bactérias por
ml.
Ludwig Krippahl, 2009
77
Processar dados experimentais

Os dados estão em 20 ficheiros 1.txt a 20.txt
Dados de crescimento
Meio:Normal
Preparador:Ana
25;0
59;1
...
296;40
Ludwig Krippahl, 2009
78
Processar dados experimentais

Objectivo: ajustar o modelo de crescimento às 2
condições e comparar os parâmetros
• Ler os ficheiros para uma lista de estruturas
• Separar as medições por meio e/ou preparador,
•
em matriz
Calcular parâmetros.
Ludwig Krippahl, 2009
79
Processar dados experimentais

Ler os ficheiros para uma lista de
estruturas
function dados=leficheiros(num)
Ludwig Krippahl, 2009
80
Ler os ficheiros
dados=[];
for f=1:num
fid=fopen([num2str(f),".txt"],"r");
reg.valores=[];
while !feof(fid)
Percorre o número indicado
...
de ficheiros numero.txt
endwhile
dados(f)=reg;
fclose(fid);
endfor
Ludwig Krippahl, 2009
81
Ler os ficheiros
dados=[];
for f=1:num
fid=fopen([num2str(f),".txt"],"r");
reg.valores=[];
while !feof(fid)
...
endwhile
Matriz para os valores neste
dados(f)=reg;
ficheiro
fclose(fid);
endfor
Ludwig Krippahl, 2009
82
Ler os ficheiros
dados=[];
for f=1:num
fid=fopen([num2str(f),".txt"],"r");
reg.valores=[];
while !feof(fid)
...
Ler o ficheiro
endwhile
dados(f)=reg;
fclose(fid);
endfor
Ludwig Krippahl, 2009
83
Ler os ficheiros
dados=[];
for f=1:num
fid=fopen([num2str(f),".txt"],"r");
Acrescenta registo à lista e
reg.valores=[]; fecha o ficheiro
while !feof(fid)
...
endwhile
dados(f)=reg;
fclose(fid);
endfor
Ludwig Krippahl, 2009
84
Ler os ficheiros
s=fgetl(fid);
if !isstr(s)
break
Lê uma linha e testa se o
endif
if findstr(s,"Meio:")!=[] resultado é string. Se não for
reg.meio=s(6:length(s));
é por ser -1, o que indica que
elseif findstr(s,"Preparador:")!=[]
reg.prep=s(12:length(s));
não há linha para ler. Nesse
elseif findstr(s,";")!=[]
caso termina o ciclo (pode
m=split(s,";");
haver linhas vazias no final
reg.valores=[reg.valores;str2num(m(1,:)),str2num(m(2,:))];
endif;
do texto).
endwhile
Ludwig Krippahl, 2009
85
Ler os ficheiros
s=fgetl(fid);
“Meio:” indica que se segue
if !isstr(s)
meio (Normal ou Meticilina)
break
endif
if !isempty(findstr(s,"Meio:"))
reg.meio=deblank(s(6:length(s)));
elseif !isemoty(findstr(s,"Preparador:"))
reg.prep=deblank(s(12:length(s)));
elseif !isempty(findstr(s,";"))
m=split(s,";");
reg.valores=[reg.valores;str2num(m(1,:)),str2num(m(2,:))];
endif;
endwhile
Ludwig Krippahl, 2009
86
o
Ler os ficheiros
s=fgetl(fid);
if !isstr(s)
O preparador (Ana ou
break
Carlos)
endif
if !isempty(findstr(s,"Meio:"))
reg.meio=deblank(s(6:length(s)));
elseif !isemoty(findstr(s,"Preparador:"))
reg.prep=deblank(s(12:length(s)));
elseif !isempty(findstr(s,";"))
m=split(s,";");
reg.valores=[reg.valores;str2num(m(1,:)),str2num(m(2,:))];
endif;
endwhile
Ludwig Krippahl, 2009
87
Ler os ficheiros
s=fgetl(fid);
Um ; indica que é uma linha
if !isstr(s)
com os valores. Split, depois
break
endif
acrescenta à matriz.
if !isempty(findstr(s,"Meio:"))
reg.meio=deblank(s(6:length(s)));
elseif !isemoty(findstr(s,"Preparador:"))
reg.prep=deblank(s(12:length(s)));
elseif !isempty(findstr(s,";"))
m=split(s,";");
reg.valores=[reg.valores;str2num(m(1,:)),str2num(m(2,:))];
endif;
endwhile
Ludwig Krippahl, 2009
88
Ler os ficheiros
Exemplo:
octave:25> l=leficheiros(20);
octave:26> l
l=
(
[1] =
{
meio = Normal
prep = Ana
valores =
34 0
...
304 42
}
[2] = ...

Ludwig Krippahl, 2009
89
Organizar os dados
Queremos receber uma matriz tempo,
contagens para cada meio e/ou
preparador:
function
mat=compiladados(lista,prep,meio)

Recebe a lista e dois strings com o
preparador e meio (“” para qualquer um)
Ludwig Krippahl, 2009
90
Organizar os dados
function mat=compiladados(lista,prep,meio)
mat=[];
for f=1:length(lista) Percorre a lista elemento
a elemento.
reg=lista(f);
if (isempty(meio) || strcmp(meio,reg.meio))
&& (isempty(prep)|| strcmp(prep,reg.prep))
mat=[mat;reg.valores];
endif
endfor
Ludwig Krippahl, 2009
91
Organizar os dados
function mat=compiladados(lista,prep,meio)
mat=[];
Se é este o meio ou
for f=1:length(lista) preparador, ou se “”,
reg=lista(f);
acrescenta.
if (isempty(meio) || strcmp(meio,reg.meio))
&& (isempty(prep)|| strcmp(prep,reg.prep))
mat=[mat;reg.valores];
endif
endfor
Ludwig Krippahl, 2009
92
Organizar os dados

Exemplo: todos os dados
l=leficheiros(20);
dados=compiladados(l,"","");
clearplot;
plot(dados(:,1),dados(:,2),"o");
Ludwig Krippahl, 2009
93
Organizar os dados

Exemplo: todos os dados
Ludwig Krippahl, 2009
94
Organizar os dados

Exemplos: separados por meio
dmet=compiladados(l,"","Meticilina");
dsem=compiladados(l,"","Normal");
clearplot
hold on
plot(dmet(:,1),dmet(:,2),"or;Meticilina;");
plot(dsem(:,1),dsem(:,2),"xg;Normal;");
Ludwig Krippahl, 2009
95
Organizar os dados
Ludwig Krippahl, 2009
96
Organizar os dados

Exemplos: separados por preparador
ana=compiladados(l,"Ana","");
carlos=compiladados(l,"Carlos","");
clearplot
hold on
plot(ana(:,1),ana(:,2),"or;Ana;");
plot(carlos(:,1),carlos(:,2),"xg;Carlos;");
Ludwig Krippahl, 2009
97
Organizar os dados

Exemplos: separados por preparador
Ludwig Krippahl, 2009
98
Ajustar o modelo

Separamos por meio:
dmet=compiladados(l,"","Meticilina");
dsem=compiladados(l,"","Normal");

E minimizamos, a partir de uma
estimativa inicial:
xs=multimin("errocres",dsem,[0.05,0.002],1e-4);
xm=multimin("errocres",dmet,[0.05,0.002],1e-4);
Ludwig Krippahl, 2009
99
Ajustar o modelo

Simulamos com os parâmetros
calculados:
sims=crescimento(xs(1),xs(2),10,0.1,400);
simm=crescimento(xm(1),xm(2),10,0.1,400);
Ludwig Krippahl, 2009
100
Ajustar o modelo
E comparamos os dados com a
simulação:
clearplot
hold on
plot(dmet(:,1),dmet(:,2),"or");
plot(dsem(:,1),dsem(:,2),"xg");
plot(simm(:,1),simm(:,2),"-r;Meticilina;");
plot(sims(:,1),sims(:,2),"-g;Normal;");

Ludwig Krippahl, 2009
101
Ajustar o modelo

Compara-se no gráfico:
Ludwig Krippahl, 2009
102
Exportar o resultado
Escrita formatada: fprintf
fprintf(id, formato, dados)
 Exemplo: escrever tabela em duas
colunas separadas por tab.
fid=fopen("relatorio.txt","w");
mat=compiladados(l,"","");
fprintf(fid,"%i\t%i\r\n",mat’)
fclose(fid)

Ludwig Krippahl, 2009
103
Exportar o resultado
fprintf(fid,"%i\t%i\r\n",mat’)
%i indica que é um número inteiro.
Ludwig Krippahl, 2009
104
Exportar o resultado
fprintf(fid,"%i\t%i\r\n",mat’)
\t é o caracter tab.
Ludwig Krippahl, 2009
105
Exportar o resultado
fprintf(fid,"%i\t%i\r\n",mat’)
\r\n indica uma nova linha.
Ludwig Krippahl, 2009
106
Exportar o resultado



% é um caracter especial na string de
formatação, indica que o que se segue
especifica o formato (ver fprint no
manual)
\ é um caracter especial em qualquer
string, usado para caracteres que não
são visíveis (mudar de linha, tab, etc.)
Para mostrar escrever dois: \\, %%
Ludwig Krippahl, 2009
107
Exportar o resultado


Quando o fprintf ou o printf (para
escrever no ecrã) recebem uma matriz
percorrem todos os elementos da matriz
aplicando a formatação na ordem
indicada;
Atenção: percorre a primeira coluna
toda, depois a segunda, etc.. Como
temos dados em colunas, usar
transposta.
Ludwig Krippahl, 2009
108
Trabalho 2

O trabalho 2 é parecido.
• Ler os ficheiros com dados e modelos.
• Ajustar as constantes.
• Exportar resultados
Ludwig Krippahl, 2009
109
Trabalho 2

Problema
• Medicamento no estômago (ME)
• Medicamento no sangue (MS)
ME  MS  Excretado
Duas constantes
Ludwig Krippahl, 2009
110
Trabalho 2

Problema
• Medicamento no estômago (ME)
• Medicamento no sangue (MS)
d[ME]/dt = -[ME] ka
d[MS]/dt = [ME] ka – [MS] ke
[ ] mg/Kg
Ludwig Krippahl, 2009
111
Trabalho 2

Problema
• Medicamento no estômago (ME)
• Medicamento no sangue (MS)
• As constantes podem depender
• Da idade (<30 anos, >=30 anos)
• De tomar em jejum ou depois da refeição
Ludwig Krippahl, 2009
112
Trabalho 2

Objectivos
• Determinar as constantes para os 4 casos
•
(idade e quando toma)
Determinar a dose necessária (mg/Kg) para
manter uma concentração média de 50mg/Kg
durante 6 dias tomando de 12 em 12 horas,
conforme o caso.
Ludwig Krippahl, 2009
113
Trabalho 2

O que fazer:
• Ler os ficheiros (estruturas?)
Ludwig Krippahl, 2009
114
Trabalho 2
Idade:49
Administrado em jejum
Horas [MS]
2.36 69
Horas: tempo da amostra
4.06 74
[MS]: mg/Kg
6.35 56
...
Ludwig Krippahl, 2009
115
Trabalho 2

O que fazer:
• Ler os ficheiros (estruturas?)
• Calcular as constantes que ajustam o modelo
aos dados de cada ficheiro
• Minimizar para as duas constantes: 18 valores
• Assumir ka dado, minimizar só uma: 15 valores
Ludwig Krippahl, 2009
116
Trabalho 2

O que fazer:
• Ler os ficheiros (estruturas?)
• Calcular as constantes que ajustam o modelo
•
aos dados de cada ficheiro
Calcular as doses necessárias para uma
média de 50mg/Kg, 12/12h, 6 dias
• Opcional, 2 valores
• Não explico como...
Ludwig Krippahl, 2009
117
Trabalho 2

O que fazer:
• Ler os ficheiros (estruturas?)
• Calcular as constantes que ajustam o modelo
•
•
aos dados de cada ficheiro
Calcular as doses necessárias para uma
média de 50mg/Kg, 12/12h, 6 dias
Gravar tudo para 4 ficheiros com os relatórios
• Constantes, dosagem, dados experimentais
Ludwig Krippahl, 2009
118
Próximas


Prática
•
•
Ajustar as curvas de crescimento
2 aulas para o trabalho
Teóricas
•
•
•
20: Excel, revisões
27: Introdução à bioinformática
3: Estrutura do exame, revisões, dúvidas.
Ludwig Krippahl, 2009
119
Dúvidas
Ludwig Krippahl, 2009
120