 <h1><center>Guia básico para experimental no Python</center></h1>
 

 
 # 1. Bibliotecas
 
Primeiramente é importante verificar corretamente quais bibliotecas importar, em geral importamos: 


```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt 
import scipy.optimize as optimization
import scipy.odr.odrpack as odrpack
```

### Numpy
$*$ Numpy é uma biblioteca usada para administrar arrays, sendo esses "vetores" que podem ser interpretados como as coordenadas no $\mathbb{R}^{n}$. Então, você pode colocar seus dados em arrays para fazer operações matemáticas com eles, por exemplo, podemos criar um array com incertezas de dados e operar com eles.

$*$ Também podemos usar ferramentas matemáticas com muita precisão, por exemplo, Cos(x), Sin(x), Sinh(x) e MUITOS outros.

$*$ Cálculo numérico ( integração, solução de sistemas lineares e não lineares, solução de EDO,etc)

$*$ Todas as informações podem ser encontradas na referência da biblioteca: https://numpy.org/doc/stable/reference/ ,
ou também nos StackOverflow da vida


```python
incertezas = np.array([0.5,0.6,0.4,0.3,0.6])
incertezas = incertezas*2
incertezas
```

Podemos mudar apenas uma componente no array, ou então, podemos mudar conjuntos dentro do array.

(Perceba que os valores estão mudados na célula de baixo devido as operações acima)

(Perceba também que o python não conta a terceira coordenda do vetor na separação, ou seja, o intervalo é aberto no final)


```python
incertezas[1:3] = incertezas[1:3]*np.cos(2)
incertezas
```

Podemos separar o array nas componentes que nos interessam. Isso se chama "Slicing".
###### Quando declarar uma variável como um slicing de um array temos um cópia, ou seja, no nosso caso, mudar o qualquer elemento do "pedaco" mudará o mesmo elemento na variável "incerteza". Para evitar isso, basta fazer uma cópia com np.copy(array)


```python
pedaco = incertezas[:4]
pedaco
```

Vou usar numpy para criar vetores com valores para embarcar no ajuste e no gráfico, então, vou definir os valores de x, y e sigmaX.


```python
x = np.linspace(0,100,10)

```

### Modelo


```python
def f(x, b, a):
    return a*x + b

# Criação dos valores de y
y = f(x,0.5,1)

#criando valores aleatórios entre o range [0.5,0.6]
sigma_y = np.random.uniform(low=2, high=5, size=(10,))
```

## Scipy

$*$ Essa biblioteca é um grupo de várias bibliotecas, nela temos quase todas as ferramentas numéricas para integração, Álgebra Linear, Analise de Fourier, estatística, etc.

$*$ A referência https://docs.scipy.org/doc/scipy/reference/index.html

### Scipy.optimization

$*$ Scipy.optimization é uma sub-biblioteca do scipy para optimização de parâmetros, no caso, usaremos o curve_fit, que serve para ajuste de funções lineares e não lineares usando MMQ, da mesma forma que o webROOT(para apenas incertezas em y), além de fornecer a matriz de covariância.
###### $*$ O método curve_fit só considera o $\sigma_y$, uma vez que utiliza o MMQ ordinário. Para incerteza em ambas as variáveis temos o MMQ total, que é feito com o método ODR do Scipy. Isso será abordado na última parte do notebook! ( parte prática).

$*$ A referência https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html


```python
#xo é o chute dos parâmetros
xo = [0,0]
fit = optimization.curve_fit(f, x, y, p0=xo, sigma=sigma_y,absolute_sigma=True)  
fit
```

Agora perceba que o "fit" devolvou um dois arrays, como se fosse uma lista contendo dois arrays. O primeiro array é o chute do parâmetros e o segundo array é a matriz de covariância( note que esse array tem 2 colunas e 2 linhas). Podemos apresentar tudo isso de maneira mais agradável, como será feito na próxima sessão com a biblioteca Pandas.

(Você pode calcular o $\chi^2$ sem biblioteca, visto que é bem simples). Equação em baixo só para mostrar como usa latex aqui.

$$
\begin{align}
{\chi}^2=\sum_{k=1}^{n} \frac{(y_i - f(x_i))^2}{\sigma_y^2}
\end{align}
$$


```python
a = fit[0][1]
b = fit[0][0]
#np.diag pega a diagonal da matriz de covariância, outro artifício do numpy
stdevs = np.sqrt(np.diag(fit[1]))
#stdevs é um array contendo duas componentes stdevs[0] é a incerteza do parâmetro (b) e stdevs[1] do parâmetro (a)
Chi2 = chi = sum((y - f(x, b, a))**2/(sigma_y**2))
Chi2
```

## Pandas

$*$ Essa biblioteca se preocupa em manipular dados de forma estrutural, por exemplo, criando e operando "tabelas" que nós chamamos de DATA FRAME. Mas ela é muito mais complexa do que isso, tendo muita aplicação em machine learning( numpy também tem bastante), então, só vou abordar de forma bem superfícial, mas encorajo você a pesquisar e aprender a mexer em DATA FRAMES.

$*$ Vou usa-la para criar uma tabela dos coeficientes, uma aplicação bem porca para a biblioteca rs. Em geral eu monto os dados numa planilha do Google Sheets, baixo em csv e passo para o python em forma de Data Frame, assim fica muito fácil trabalhar. O pandas faz com que as colunas de Data Frames se comportem como arrays, dessa forma, posso fazer operações matemáticas com colunas e linhas. 

$*$ O guia/referência é : https://pandas.pydata.org/docs/user_guide/index.html


```python
df = pd.DataFrame({'Parâmetro(a)' : [a,stdevs[1]], 'Parâmetro(b)' : [b,stdevs[0]]}, index=['Valor', 'Incerteza'])
df
```

### Matplotlib

$*$ O nome da biblioteca fala por ela. Para plotar todo tipo de coisa e ainda dá configurar ao máximo os gráficos.

$*$ Tem muito conteúdo na internet e a referencia dela é bacana e recheada de exemplo https://matplotlib.org/tutorials/index.html#introductory, novamente o StackOverflow é uma ótima fonte de pesquisa.


```python
h = plt.figure(1)
h = plt.plot(x,f(x,b,a), label="model", color='black')
h = plt.errorbar(x,y,sigma_y,fmt='.',color='red')
h = plt.grid()
h = plt.ylabel(r'$f(x)$')
h = plt.xlabel(r'$x$')
#para salvar o gráfico no mesmo diretório
plt.savefig('teste_1.png', format='png', dpi=150,bbox_inches = "tight")
plt.show(h)
```

Agora o gráfico de resíduos. (Lembrando, que estes são os gráficos defaults mais simples possíveis, caso queira sofisticar é só pesquisar, tem muito conteúdo!)


```python
g = plt.figure(2)
g.set_size_inches(6, 1)
g =plt.errorbar(x,y-f(x,b,a),sigma_y,fmt='.',color='red')
g = plt.hlines(0,x.min(),x.max(),color='black')
g = plt.grid()
g = plt.ylabel("f (x)")
g = plt.xlabel('x')
g = plt.savefig('teste_1_resid.png', format='png', dpi=150,bbox_inches = "tight")
g = plt.show(g)
```

## Dicas adicionais

Quero relembrar que todas as bibliotecas são vastas e os conteúdos sobre elas não se acabam, então, sempre que tiver dúvida você provavelmente vai achar com bastante facilidade na internet, principalmente, em referências e no StackOverflow.

$(*)$ Ótimo vídeo sobre o básico de curve fitting no python: 
https://www.youtube.com/watch?v=Jl-Ye38qkRc&t=810s

$(*)$ Página com os códigos e noções básicas sobre o curve fitting no python: 
https://towardsdatascience.com/basic-curve-fitting-of-scientific-data-with-python-9592244a2509


$(*)$ Eu usual faço tabelas no Google Sheets e import para cá com o pandas. Não é difícil fazer e nem mexer, mas exige um pouco para se acostumar e tals.

# 2. Colocando na prática

#### Vamos usar dados com incerteza em x e y ( como geralmente ocorre na disciplina de experimental)

$*$ É importante lembrar que a função curve_fit vista acima é apenas para o caso com incerteza em x

$*$ Para simular uma atividade foram gerados dados que estão em dois formatos, csv e txt.

#### Importando txt com Pandas


```python
data = pd.read_csv('exemplo.txt',sep="\t")

# No caso do separador com vírgulas é preciso fazer
# data_1 = pd.read_csv('exemplo.txt', sep="\t", decimal=",")

data.head()
```

#### Importando csv com Pandas


```python
data = pd.read_csv('exemplo.csv')
data.head()
```

### Análise dos dados

Vamos começar pelo modelo ( os dados criados são feitos para serem ajustados por uma parábola)

#### Modelo


```python
# Importante que notar que agora nossos parâmetros são representados por array
# Outro fator a se notar é a ordem na a f, ou seja, f(p,x) "p" vem primeiro

def f(p,x):
    return p[0]*x**2 + p[1]*x + p[2]
```

##### Usando os valores de $x$, $y$, $\sigma_x$ e $\sigma_y$


```python
x = data['x (cm)']
y = data['y (cm)']
sigma_x = data['sigma_x (cm)']
sigma_y = data['sigma_y (cm)']
# O x é uma "Série", seria algo como um array com indexes, mas se comporta de maneira
# bastante similar

x.head()
```

### O ajuste

Como agora nossa incerteza está nas duas variáveis, precisamos utiliar um outro método da mesma biblioteca( Scipy ). O método é o ODR.


```python
# objeto modelo ( isso é um termo técnico para dizer que o modelo está guardado em um "objeto")
model_object = odrpack.Model(f)
# Objeto data ( termo técnico para dizer que dizer que o modelo está guardado em um "objeto")
data_object = odrpack.RealData(x,y,sx=sigma_x, sy=sigma_y)
# Estabelece a relação entre data e modelo, e em geral necessita de das condições iniciais
# No caso do nosso modelo temos 3 parâmetros, ou seja beta0 = lista de 3 elementos
odr = odrpack.ODR(data_object, model_object, beta0=[1,1,1])
# Roda a regressão
odr.set_job(fit_type=0)  # fit_type = 0 executa a regressão por ODR e fit_type = 2 por MMQ
# O tópico do fit_type é um complicado e vai além do proposito dessa introdução
out = odr.run()
out.pprint()
```

Agora só falta implementar o $\chi^2$

$$
\begin{align}
{\chi}^2=\sum_{k=1}^{n} \frac{(y_i - f(x_i))^2}{\sigma_i^2 + \big(\frac{\partial f}{\partial x}\big)^2 \sigma_{x_i}^2 }
\end{align}
$$

Não impletarei, mas você pode fazer essa linha de código fácilmente sabendo a derivada analítica

### Gráficos

Vale ressaltar que o método plt.savefig() salvará os gráficos no diretório do .ipynb ou do .py


```python

ax2 = plt.figure(1)
ax2 = plt.errorbar(x, y, xerr=sigma_x, yerr=sigma_y, linestyle='None', marker='+',color='black')
#ax2 = plt.ylim(-1,2)
#ax2 = plt.xlim(0, 20000)
ax2 = plt.plot(x, f([2.01584576e-01, -3.87102908e+01, 1.85196706e+03],x), label='model',color='red')
# nomeando os eixos
ax2 = plt.xlabel(r'$x (cm)$')
ax2 = plt.ylabel(r'$y (cm)$')
ax2 = plt.grid()
#Salvando a figura
ax2 = plt.savefig('exemplo.png',dpi=150,bbox_inches = "tight")

ax3 = plt.figure(2)
ax3 = plt.errorbar(x, abs(y - f([2.01584576e-01, -3.87102908e+01, 1.85196706e+03],x))
                   ,xerr=sigma_x, yerr=sigma_y, linestyle='None', marker='+',color='black')

ax3 = plt.hlines(0,x.min(),x.max(),color='red')
ax3 = plt.grid()
# Nomeando os eixos
ax2 = plt.xlabel(r'$x (cm)$')
ax2 = plt.ylabel(r'$y(cm) - modelo$')
ax3 = plt.gcf()
ax3.set_size_inches(6, 1)
#Salvando a figura
ax3 = plt.savefig('exemplo_resid.png',dpi=150,bbox_inches = "tight")
```

4e447a2b449a034f9ccfdeac4e831c901e5b7f69
ipynb
Jupyter Notebook
Notebook-pt.ipynb
Rodrigo-Motta/Intro-Data-Analysis-Exp-Physis-pt-en
d97b50a9739ca3825ff149d607e64a7342a511da
MIT
Notebook-pt.ipynb
Rodrigo-Motta/Intro-Data-Analysis-Exp-Physis-pt-en
d97b50a9739ca3825ff149d607e64a7342a511da
MIT
Notebook-pt.ipynb
Rodrigo-Motta/Intro-Data-Analysis-Exp-Physis-pt-en
d97b50a9739ca3825ff149d607e64a7342a511da
MIT
Qwen/Qwen-72B
1. YES
2. YES