12 - Métodos Estocásticos12 - Métodos Estocásticos12 ...ICE-BICE-BICE-BICE-BICE-BICE-BICE-BICE-BICE-B 12 - Métodos Estocásticos12 - Métodos Estocásticos12 - Métodos Estocásticos12

ICE-BICE-BICE-BICE-BICE-BICE-BICE-BICE-BICE-B

12 - Métodos Estocásticos12 - Métodos Estocásticos12 - Métodos Estocásticos12 - Métodos Estocásticos12 - Métodos Estocásticos12 - Métodos Estocásticos12 - Métodos Estocásticos12 - Métodos Estocásticos12 - Métodos Estocásticos

Ludwig Krippahl

1

Métodos EstocásticosMétodos EstocásticosMétodos EstocásticosMétodos EstocásticosMétodos EstocásticosMétodos EstocásticosMétodos EstocásticosMétodos EstocásticosMétodos Estocásticos

Resumo■ Calculo de integrais por Monte Carlo■ Ajuste de modelos: avaliação de parâmetros com bootstrapping

2


Monte CarloMonte CarloMonte CarloMonte CarloMonte CarloMonte CarloMonte CarloMonte CarloMonte Carlo

3

ImagensImagensImagensImagensImagensImagensImagensImagensImagens

Métodos de Monte Carlo■ Por alusão aos casinos■ Métodos de amostragem aleatória para estimar estatísticas■ Exemplo: estimar uma área• Sabemos distinguir entre pontos dentro e fora da região de interesse

• Lançamos pontos ao acaso numa "caixa" de área conhecida que inclui aquela quequeremos determinar

• Contamos a proporção que cai na região de interesse

• A área dessa região é dada pela proporção de pontos lá dentro

4


Calcular um integral

■ É fácil determinar se umponto está abaixo ouacima

■ Mas não sabemos aexpressão da primitiva

∫x1

x0

esin(x)

5



■ Vamos "disparar" pontosdentro da caixa

■ Contar quantos calhamabaixo da função

∫x1

x0

esin(x)

6



■ Exemplo:• Calharam 3677 de 10000• Caixa tem 8x4• Integral:

∫8

0

esin(x)

8 × 4 × 3677/10000 = 11.77

7


Calcular um integral■ Implementação genérica, vamos separar tarefas:• Uma função faz a estimativa por Monte Carlo

• Recebe como argumento a função a integrar, a "caixa" e o número de pontos

import numpy as np def mc_integral_fn(fn,box,num_points): "Estimate integral of fn by Monte Carlo" in_count = 0 x0,x1,y0,y1=box for _ in range(num_points): x=np.random.uniform(x0,x1) y=np.random.uniform(y0,y1) y_f = fn(x) if y<=y_f: in_count = in_count+1 return in_count/num_points*(x1-x0)*(y1-y0)

■ Nota: uma função é um objecto como qualquer outro

8


Calcular o integral ∫ x1

x0e

sin(x)

• Para usar a função precisamos de criar a função cujo integral queremos estimar

def mc_integral_fn(fn,box,num_points): "Estimate integral of fn by Monte Carlo" ... def exp_sin(x): return np.e**np.sin(x)

■ Agora fornecemos esta função à mc_integral_fn• Nota: sem os parênteses, porque isso é para mandar executar

In : mc_integral_fn(exp_sin, (0,8,0,4), 10000) Out: 11.3792 In : mc_integral_fn(exp_sin, (0,8,0,4), 10000) Out: 11.5296

9


Para calcular outros basta criar as funções

∫1

0

1 − x4

‾ ‾‾‾‾‾√

def sqrt_x4(x): return (1-x**4)**0.5

In : mc_integral_fn(sqrt_x4, (0,1,0,1.2), 10000) Out: 0.8737199999999999

10


Para calcular outros basta criar as funções

∫4

2

1

ln(x)

def one_ln(x): return 1/np.log(x)

In : mc_integral_fn(one_ln, (2,4,0,2), 10000) Out: 1.914

11


■ Para fazer o gráfico precisamos dos pontos dentro e fora■ Vamos tornar o gráfico opcional, e delegá-lo noutra função■ Primeiro, a função do gráfico, recebendo listas de tuplos (x,y)

def plot_mc(fn,box,ins,outs): "plots MC result" ins = np.array(ins) #better as numpy arrays outs = np.array(outs) plt.figure(figsize=(10,10)) #create a square figure plt.plot(ins[:,0],ins[:,1],'.b') #in is blue plt.plot(outs[:,0],outs[:,1],'.r') #out is red x_axis = np.linspace(box[0],box[1],200) #plot the function plt.plot(x_axis,fn(x_axis),'-k',linewidth=3) plt.axis(box) #set plot to the box size plt.show()

12


■ Agora alteramos a função original para desenhar o gráfico

def mc_integral_fn(fn,box,num_points,do_plot=False): "Estimate integral of fn by Monte Carlo" in_count = 0 ins = [] # listas onde guardar os pontos outs = [] x0,x1,y0,y1=box for _ in range(num_points): x=np.random.uniform(x0,x1) y=np.random.uniform(y0,y1) y_f = fn(x) if y<=y_f: ins.append((x,y)) in_count = in_count+1 # guarda na lista de dentro else: outs.append((x,y)) # guarda na lista de fora if do_plot: # por omissão não faz plot_mc(fn,box,ins,outs) return in_count/num_points*(x1-x0)*(y1-y0)

13


Ajuste de modelos e BootstrappingAjuste de modelos e BootstrappingAjuste de modelos e BootstrappingAjuste de modelos e BootstrappingAjuste de modelos e BootstrappingAjuste de modelos e BootstrappingAjuste de modelos e BootstrappingAjuste de modelos e BootstrappingAjuste de modelos e Bootstrapping

14

Ajuste de modelosAjuste de modelosAjuste de modelosAjuste de modelosAjuste de modelosAjuste de modelosAjuste de modelosAjuste de modelosAjuste de modelos

Cinética, Michaelis-Menten, v =[S]Vmax

+[S]KM

■ Temos um ficheiro com os dados ( e )[S] v

1,2,5,8,12,30,50

11.1,25.4,44.8,54.5,58.2,72.0,60.1

15



+[S]KM

■ Podemos fazer o ajuste directamente minimizando o erroquadrático

• Não precisamos de linearizar com logaritmos

• Podemos considerar o erro directamente nos dados em vez de no logaritmo

■ Primeiro, definimos a função a minimizar (o "custo")

import numpy as np def mm_cost(pars,S,v): """returns quadratic error of Michaelis Menten model""" pred = pars[0]*S/(pars[1]+S) cost = np.sum((v-pred)**2) return cost

16



+[S]KM

■ Depois usamos a função minimize do scipy

from scipy.optimize import minimize ... def adjust_model(data): """returns adjusted MM parameters from data""" S = data[0,:] v = data[1,:] res=minimize(mm_cost,[1,1],args=(S,v),bounds=[(0,None),(0,None)]) vmax,Km=res.x return vmax,Km

■ Os argumentos de minimize são:• A função a minimizar, valores iniciais, argumentos a passar à função de custo e

limites para os valores

■ Devolve um objecto com vector x com os valores óptimosencontrados

17


■ Lemos os dados e ajustamos o modelo (com função de plot):

In : data = np.loadtxt('mm.txt',delimiter=',') In : vmax, Km = adjust_model(data) Out: 73.261345790320036, 3.437154134510513

18


Residual BootstrappingResidual BootstrappingResidual BootstrappingResidual BootstrappingResidual BootstrappingResidual BootstrappingResidual BootstrappingResidual BootstrappingResidual Bootstrapping

19

BootstrappingBootstrappingBootstrappingBootstrappingBootstrappingBootstrappingBootstrappingBootstrappingBootstrapping

Estimar o intervalo de confiança dos parâmetros■ Bootstrapping é um método de reamostragem aleatória• Permite estimar estatísticas repetindo virtualmente a experiência

■ Normalmente, criamos réplicas amostrando os dados comreposição

■ Neste caso, não pode ser os dados porque não são aleatórios• Os valores de [S] são escolhidos

■ Mas podemos assumir que os erros são aleatórios■ Amostrar os erros é Residual Bootstrapping

20


■ Vamos medir os erros (residuais) de cada ponto em relação aomelhor modelo

21


■ Depois vamos criar novos conjuntos reamostrando os residuais

22


■ E usamos cada conjunto para ajustar novas curvas (e parâmetros)

23


■ E usamos cada conjunto para ajustar novas curvas (e parâmetros)

def residual_bootstrapping(N,vmax,Km,data): "creates N replicas by residual bootstrapping and return parameters" curve = vmax*data[0,:]/(Km+data[0,:]) residuals = data[1,:]-curve params=[] for ix in range(N): rnd = np.random.randint(0,data.shape[1],data.shape[1]) replica=curve+residuals[rnd] res=minimize(mm_cost,[1,1],args=(data[0,:],replica), bounds=[(0,None),(0,None)]) params.append(res.x) return np.array(params)

In : params = residual_bootstrapping(500,vmax,Km,data) In : np.mean(params,axis=0) Out: [ 73.06719906 3.53974319] In : np.std(params,axis=0) Out: [ 4.19434788 0.74514656]

24


■ Assim temos uma ideia da incerteza

25


■ Também podemos ver os histogramas dos parâmetros

plt.figure() plt.hist(params[:,0],color = (1,0.5,0.3),linewidth=3) plt.figure() plt.hist(params[:,1],color = (0.7,0.7,1),linewidth=3)

26


ResumoResumoResumoResumoResumoResumoResumoResumoResumo

27


Objectivos de hoje■ Fazer plots (plot bar hist)■ Ficar com uma ideia melhor da utilidade de ICE• Um semestre de programação não dá para muito

• Mas se usarem ICE como ponto de partida podem aprender coisas úteis

Leitura adicional:■ Recomendada: Capítulo 12 dos apontamentos■ Opcional:• Bootstrapping: https://en.wikipedia.org/wiki/Bootstrapping_(statistics)

• SciPy: https://www.scipy.org/

• Monte Carlo Methods: https://en.wikipedia.org/wiki/Monte_Carlo_method

Documents

12 - Métodos Estocásticos12 - Métodos Estocásticos12 ...ICE-BICE-BICE-BICE-BICE-BICE-BICE-BICE-BICE-B 12 - Métodos Estocásticos12 - Métodos Estocásticos12 - Métodos Estocásticos12