Skip to content

Commit

Permalink
seccao de propagação de erros
Browse files Browse the repository at this point in the history
  • Loading branch information
viniciusdutra314 committed Jan 21, 2025
1 parent 72f98f2 commit a9076ae
Show file tree
Hide file tree
Showing 22 changed files with 217 additions and 80 deletions.
4 changes: 4 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
#temporario
*.jpg
*.png

# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
Expand Down
2 changes: 1 addition & 1 deletion LabIFSC2/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
'''


from ._arrays import (arrayM, curva_max, curva_min, incertezas, linspace,
from ._arrays import (arrayM, curva_max, curva_min, incertezas, linspaceM,
nominais)
from ._medida import Comparacao, Medida, comparar_medidas, montecarlo
from ._regressões import (MExponencial, MPolinomio, regressao_exponencial,
Expand Down
2 changes: 1 addition & 1 deletion LabIFSC2/_arrays.py
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ def curva_max(arrayMedidas : np.ndarray,unidade:str,sigma:float | int=2)-> np.nd


@obrigar_tipos
def linspace(a:Real,b:Real,n : int,
def linspaceM(a:Real,b:Real,n : int,
incertezas : Real ,unidade : str) -> np.ndarray:
"""Gera um array com N Medidas de valor nominal [a,b]
A incerteza será constante caso 'incertezas' for um número,
Expand Down
5 changes: 2 additions & 3 deletions LabIFSC2/_medida.py
Original file line number Diff line number Diff line change
Expand Up @@ -184,9 +184,8 @@ def __getattr__(self, func_name:str) -> Any:
raise AttributeError
else:
func=getattr(np,func_name)
def FuncaoComMedida() -> Any: return montecarlo(func, self)
return FuncaoComMedida

def funcao_recebe_medida() -> Any: return montecarlo(func, self)
return funcao_recebe_medida
def _adicao_subtracao(self,outro: 'Medida',positivo:bool) -> 'Medida':
if not (isinstance(outro,Medida) or isinstance(outro,Real)):
return NotImplemented
Expand Down
2 changes: 0 additions & 2 deletions docs/Arrays/linspace.md
Original file line number Diff line number Diff line change
@@ -1,2 +0,0 @@
:::LabIFSC2._arrays.linspace

7 changes: 4 additions & 3 deletions docs/arrays.md
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,13 @@ entre arrays, o resultado disso são operações elemento a elemento
```

## Operações matemáticas
Como discutido na secção de [Funçõs matemáticas](funcoes_matematicas.md), as funções
matemáticas da biblioteca (lab.exp,lab.cos,lab.sinh), podem atuar diretamente em arrays
Como discutido na secção de [Funçõs matemáticas](funcoes_matematicas.md), as funções do numpy podem atuar diretamente na classe Medida.

```py
--8<-- "tests/test_doc_sqrt_vetorizado.py:1:8"
```

## linspace
## linspaceM
Em muitas medidas experimentais nós fazemos medições igualmente espaçadas, imagine que
você está medindo o campo magnético de um fio em função da sua distância \(\vec{B}(distância)\),
você realiza uma medição a cada 1cm por exemplo, o lab.linspace recebe o valor da primeira medição,
Expand All @@ -33,6 +32,8 @@ No exemplo abaixo nós fizemos 10 medições entre [1cm,10cm], com precisão de

A função é o análogo do [np.linspace](https://numpy.org/doc/stable/reference/generated/numpy.linspace.html#numpy-linspace) que recebe medidas

##

## arrayM
Agora que temos as distâncias nós medimos o campo magnético pra cada distância, como registramos esses campos?
teríamos que criar 10 objetos `Medida` diretamente? Não, a solução é o lab.arrayM, ele recebe uma lista/array
Expand Down
26 changes: 5 additions & 21 deletions docs/funcoes_matematicas.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,31 +6,15 @@ portanto, operações entre Medidas são idênticas a operações entre números
--8<-- "tests/test_doc_operacoes_basicas.py:1:10"
```

## Decorador Aceita Medida
Pela generalidade do método de propagação de erros usando monte carlo, o LabIFSC2 possui uma decorador
que consegue fazer com que uma função matemática **qualquer do numpy** aceite uma `Medida` e faça propagação de erros,
essa é a magia do `aceitamedida`
## Funções Numpy

O LabIFSC2 implementa uma compatibilidade direta com as funções do Numpy (tecnicamente chamadas de [ufunc](https://numpy.org/doc/stable/reference/ufuncs.html)), então você aplicar funções como np.sin, np.sqrt, np.arctanh naturalmente[^1]

```py
--8<-- "tests/test_doc_aceitamedida.py:1:11"
```

Perceba que pegamos a função do numpy (np.sin) e transformamos ela em uma função que aceita medidas,
temos ainda integrações bonitas com o sistema de unidades, conseguimos usar um ângulo diretamente em **graus**
sem converter para radianos.

## Operando em arrays

As funções além de aceitarem medidas também aceitam arrays numpy, podemos fazer então operações
vetorizadas

```py
```py
--8<-- "tests/test_doc_sqrt_vetorizado.py:1:8"
```

## Todas funções suportadas
Caso queira uma lista de todas as funções matemáticas que estão no LabIFSC2, eis o código fonte com todas elas e seus apelidos (sin=seno, tan=tg)
```py {title=_matematica.py}
--8<-- "LabIFSC2/_matematica.py:39"
```
[^1]:
Para os curiosos, isso é feito implementado métodos específicos na classe Medida, por exemplo, a função np.sqrt(x) verifica se o tipo de x tem o método x.sqrt,se tiver, ele chama x.sqrt, a classe `Medida` possui uma simulação monte carlo de sqrt e todas as funções matemáticas (isso é feito criando-se dinamicamente os métodos com [__getattr__](https://docs.python.org/3/reference/datamodel.html#object.__getattr__))
28 changes: 23 additions & 5 deletions docs/graficos.md
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,29 @@ De maneira análogo podemos também pegar as incertezas com `incertezas(array_me
Eis um exemplo simples de como fazer um gráfico de dispersão com erros tanto em x quanto em y,
basta usar a função do matplotlib [`plt.errorbar`](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.errorbar.html), `nominais` e `incertezas`.

```py
--8<-- "tests/test_doc_grafico_scatter.py:1:16"

<img src="./images/graficos_scatter.jpg" width=400>

```py hl_lines="12-16"
--8<-- "tests/test_doc_grafico_scatter.py:1:19"
```

Repare que as unidade são variáveis no código que podem ser modificadas rapidamente.



## Curva Min/Max

As funções `curva_min` e `curva_max` são utilizadas para calcular a curva teórica considerando as incertezas nos parâmetros encontrados durante a regressão dos dados (chamados de [Confidence and prediction bands]([Confidence and prediction bands)](https://en.wikipedia.org/wiki/Confidence_and_prediction_bands))

## Curva teórica com erro
Regressões de dados inevitavelmente apresentam incertezas nos parâmetros encontrados, podemos representa-las usando as funções `curva_min` e `curva_max` que calculam a curva teórica \(\pm \, \, 2\sigma \)

<img src="./images/graficos_fitting.jpg" width=400>

```py hl_lines="20-29"
--8<-- "tests/test_doc_grafico_fitting.py:1:30"
```

Repare que as unidade são variáveis no código que podem ser modificadas rapidamente, o resultado é esse:
<img src="exemplo.jpg" width=400>
Perceba que para lidar com as unidades acaba que a sintaxe fica infelizmente verbosa, isso é algo pessoal mas nesse caso eu recomendo usar um tipo de indentação chamada `Hanging indentation ` em cada argumento ocupa uma linha de código, assim o código fica mais legível (na minha opinião) e menos horizontal.

##
1 change: 0 additions & 1 deletion docs/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,6 @@ Mesmo que a interface seja intencionalmente parecida a implementação é totalm
| Unidades | Implementação autoral | Baseado no famoso [pint](https://pint.readthedocs.io/)
| Constantes da natureza| ❌ | +350 definidas pela [CODATA(2022)](https://codata.org/initiatives/data-science-and-stewardship/fundamental-physical-constants/)
| Operações com arrays| ❌ | Suportadas pelo [Numpy](numpy.org)
| Tabelas em LaTeX | ❌| ✅
| Segurança de tipos (mypy)| ❌ | ✅
| Docstrings em funções | ❌ | ✅
| Suporte || Ativo |
Expand Down
68 changes: 68 additions & 0 deletions docs/propagacao_de_erros.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
Nesta seção, explicarei em mais detalhes como a biblioteca propaga incertezas, o método usado é mais geral, mas ainda assim compatível com o da apostila. Nos testes unitários da biblioteca, comparamos os erros calculados pelo LabIFSC2 com as bibliotecas [uncertainties](https://pythonhosted.org/uncertainties/) e [LabIFSC](https://github.com/gjvnq/LabIFSC), chegando a um acordo geralmente de \(10^{-3}\) para erros pequenos, onde os métodos devem ser equivalentes.

## Apostila
A apostila se baseia principalmente no GUM (Guide to the Expression of Uncertainty in Measurement)[^1]. O método é uma propagação linear baseada em uma série de Taylor.

Começamos com uma série de Taylor de uma função de \(N\) variáveis, ou seja, uma boa aproximação para **pequenas variações da função**:

$$f(\mathbf{x}+\Delta \mathbf{x})\approx f(\mathbf{x})+\sum_{i} \frac{\partial f}{\partial x_{i}}\Delta x_{i}$$

$$(f(\mathbf{x}+\Delta \mathbf{x})-f(\mathbf{x}))²\approx (\sum_{i} \frac{\partial f}{\partial x_{i}}\Delta x_{i})²$$

Note que o lado esquerdo nada mais é que a variação da função \(\Delta f\). Se pegarmos o valor esperado, teremos a variância de \(f\) (\(\sigma²_{f}\)):

$$\mathbb{E}[(f(\mathbf{x}+\Delta \mathbf{x})-f(\mathbf{x}))²]\approx \mathbb{E}[(\sum_{i} \frac{\partial f}{\partial x_{i}}\Delta x_{i})²]$$

Supondo que \(\Delta x_i\) tenha uma distribuição simétrica \(\mathbb{E}(\Delta x_{i})=0\), os termos cruzados \(\mathbb{E}(\Delta x_{i}\Delta x_{j})=\mathbb{E}(\Delta x_{i})\mathbb{E}(\Delta x_{j})=0\) desaparecem (se supormos que são totalmente independentes).

Logo:

$$\sigma²_{f}=\sum_{i} \left(\frac{\partial f}{\partial x_{i}}\sigma_{x_{i}}\right)²$$

Se quisermos em termos do desvio padrão e não da variância, temos:

$$\sigma_{f}=\sqrt{\sum_{i} \left(\frac{\partial f}{\partial x_{i}}\sigma_{x_{i}}\right)²}$$

Essa é essencialmente a fórmula usada na apostila. Para o caso de uma variável, se reduz a \(\sigma_{f}=|\frac{\partial f}{\partial x}\sigma_{x}|\). A diferença é que a apostila ignora a raiz quadrada na expressão de incertezas com mais de uma variável, superestimando assim o erro.

Pensando intuitivamente, erros não podem simplesmente se somar, visto que, pela sua natureza aleatória, é esperado que existam erros que acabem "compensando" outros.

## Monte Carlo
Imagine que temos uma medida indireta \(y\) que depende de um conjunto de \(N\) medidas:

$$y=f(x_1,x_2,\dots,x_n)$$

Cada variável \(x_i\) tem sua PDF ([probability density function](https://en.wikipedia.org/wiki/Probability_density_function)), que é uma forma matemática de dizer que não temos certeza sobre seus valores. Por simplicidade, assumimos que medidas diretas têm distribuições gaussianas (centradas em uma média \(\mu\) e uma variância \(\sigma²\)).

O Método Monte Carlo consiste em simular diversas medidas experimentais no computador (usando um gerador de números aleatórios com as respectivas PDFs). Dessa forma, geramos um histograma de possíveis valores de \(y\); esse histograma é a PDF de \(y\).

O interessante desse método é que o histograma de \(y\) não necessariamente é analítico (geralmente com formatos bem estranhos para incertezas grandes). Esse histograma é utilizado para diversas coisas na biblioteca:

- Ser usado como PDF para outra propagação de incerteza
- Cálculo do intervalo de confiança
- A média e o desvio padrão da PDF são usados no `print(medida)`

### Exemplo com a gravidade
Retornando ao exemplo da estimativa da gravidade usando um pêndulo, mas agora com incertezas maiores em \(T\) e \(L\) (para que os efeitos fiquem mais visíveis).

A classe `Medida` possui um atributo chamado `histograma`, onde estão guardados os histogramas. No dia a dia, esse atributo deve ser raramente acessado, mas para fins didáticos ele é interessante.

```py
--8<-- "tests/test_doc_gravidade_histograma.py:8:15"
```

Repare como \(T\) e \(L\) são gaussianas (\(\mu_L=15cm\), \(\sigma_L=1cm\)) e (\(\mu_T=780ms\), \(\sigma_T=80ms\)).

<img src="./images/gravidade_histograma.jpg" width=600>

Já o histograma de \(g\) é centralizado em \(10m/s²\), mas observe que ele possui uma cauda para a direita. A distribuição não é simétrica, logo, não é gaussiana. Se usássemos \(g\) para outros cálculos, esse desvio de uma gaussiana provavelmente iria se amplificando. Esse fato não é observado no método GUM, que assume linearidade e basicamente tudo é uma gaussiana.



[^1]: O método GUM é amplamente utilizado em metrologia e calibragem de equipamentos. Existem diversas referências para quem quiser aprender mais. Eu, pessoalmente, achei um material introdutório e interessante em:

Kirkup, L., Frenkel, R. B. (2006). An Introduction to Uncertainty in Measurement: Using the GUM (Guide to the Expression of Uncertainty in Measurement). Reino Unido: Cambridge University Press.

[^2]: O principal material usado na implementação do Monte Carlo foi o próprio material suplementar do GUM sobre Monte Carlo. É interessante notar que esse material explicitamente considera o método Monte Carlo como uma forma mais precisa de calcular incertezas:

“Evaluation of Measurement Data — Supplement 1 to the ‘Guide to the Expression of Uncertainty in Measurement’ — Propagation of Distributions Using a Monte Carlo Method,” 2008. https://doi.org/10.59161/JCGM101-2008.
6 changes: 4 additions & 2 deletions mkdocs.yml
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,10 @@ nav:
- Constantes: Constantes/constantes da natureza.md
- Funções Matemáticas: funcoes_matematicas.md
- Arrays: arrays.md
- Gráficos: graficos.md
- Regressões (Fitting): Regressões (fitting)/linear.md
- Gráficos: graficos.md
- Formataçõs/LaTeX: formatacoes_latex.md
- Como erros são propagados?: propagacao_de_erros.md
- API Completa: api.md
- Como contribuir: contribuir.md

Expand Down Expand Up @@ -49,4 +50,5 @@ markdown_extensions:
- pymdownx.inlinehilite
- pymdownx.snippets
- pymdownx.superfences
- mdx_math
- mdx_math
- footnotes
8 changes: 4 additions & 4 deletions tests/test_arrayM.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
import pytest

from LabIFSC2 import (Medida, arrayM, curva_max, curva_min, incertezas,
linspace, nominais)
linspaceM, nominais)


def test_nominal():
Expand Down Expand Up @@ -44,19 +44,19 @@ def test_curvamax():

def test_linspace():
a=1 ; b=20 ; N=20
x=linspace(a,b,N,0.1,'')
x=linspaceM(a,b,N,0.1,'')
assert np.all(nominais(x,"")==np.linspace(a,b,N))
assert np.all(incertezas(x,"")==0.1)


def test_converter_array():
x=linspace(5,10,10,0.01,'mm')
x=linspaceM(5,10,10,0.01,'mm')
esperado=np.linspace(0.5,1,10) #cm
for index in range(len(x)):
assert np.isclose(x[index].nominal('cm'),esperado[index])

def test_converter_array_si():
x=linspace(10,50,10,0.01,'cm')
x=linspaceM(10,50,10,0.01,'cm')
esperado=np.linspace(0.1,0.5,10) #si
for index in range(len(x)):
assert np.isclose(x[index].nominal('si'),esperado[index])
Expand Down
8 changes: 4 additions & 4 deletions tests/test_classe_para_regressoes.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@


def test_initialization():
coeficientes = lab.linspace(1,3,3,0,'')
coeficientes = lab.linspaceM(1,3,3,0,'')
with pytest.raises(TypeError):
lab.MPolinomio(np.array([1, 2, 3]))
polinomio=lab.MPolinomio(coeficientes)
Expand All @@ -16,11 +16,11 @@ def test_initialization():
assert polinomio.c.nominal("") == 3

def test_call():
coeficientes = lab.linspace(1,3,3,0,'')
coeficientes = lab.linspaceM(1,3,3,0,'')
polinomio=lab.MPolinomio(coeficientes)
with pytest.raises(TypeError):
polinomio(0)
pontos=lab.linspace(0,2,3,0,'')
pontos=lab.linspaceM(0,2,3,0,'')
assert np.isclose(polinomio(pontos[0]).nominal(""), 3)
assert np.isclose(polinomio(pontos[1]).nominal(""),6)
assert np.isclose(polinomio(pontos[2]).nominal("") , 11)
Expand All @@ -32,7 +32,7 @@ def test_call():


def test_unpacking():
coeficientes = lab.linspace(1,3,3,0,'')
coeficientes = lab.linspaceM(1,3,3,0,'')
polinomio = lab.MPolinomio(coeficientes)
a, b, c = polinomio
assert a.nominal("") == 1
Expand Down
28 changes: 18 additions & 10 deletions tests/test_doc_grafico_fitting.py
Original file line number Diff line number Diff line change
@@ -1,25 +1,33 @@
import matplotlib.pyplot as plt
import numpy as np

from LabIFSC2 import *

campo_magnético=arrayM([210,90,70,54,39,32,33,27,22,20],1,'muT')
distancias=linspace(1,10,10,0.01,'cm')
distancias=linspaceM(1,10,10,0.01,'cm')

unidade_x='cm'
unidade_y='muT'

regressao=regressao_potencia(distancias,campo_magnético)
print(regressao)
print(regressao)
#MLeiDePotencia(a=(2,0 ± 0,1)x10² cm⁰⋅⁹⁸⁶³⁵⁵·µT, b=(-9,9 ± 0,3)x10⁻¹ )
plt.style.use('ggplot')
plt.errorbar(nominais(distancias,unidade_x),nominais(campo_magnético,unidade_y),
xerr=incertezas(distancias,unidade_x),yerr=incertezas(campo_magnético,unidade_y),
plt.errorbar(x=nominais(distancias,unidade_x),
y=nominais(campo_magnético,unidade_y),
xerr=incertezas(distancias,unidade_x),
yerr=incertezas(campo_magnético,unidade_y),
fmt='o',label='Dados experimentais',color='red')

x=linspace(1,10,100,0,unidade_x)
plt.plot(nominais(x,unidade_x),nominais(regressao(x),unidade_y),color='blue',
x=linspaceM(1,10,100,0,unidade_x)
plt.plot(nominais(x,unidade_x),
nominais(regressao(x),unidade_y),
color='blue',
label="Curva teórica")
plt.fill_between(nominais(x,unidade_x),curva_min(regressao(x),unidade_y),
curva_max(regressao(x),unidade_y),color='blue',alpha=0.3)

plt.fill_between(x=nominais(x,unidade_x),
y1=curva_min(regressao(x),unidade_y),
y2=curva_max(regressao(x),unidade_y),
color='blue',alpha=0.3)
plt.legend()
#plt.savefig('teste.jpg')
plt.savefig('docs/images/graficos_fitting.jpg',dpi=300)
plt.cla()
12 changes: 8 additions & 4 deletions tests/test_doc_grafico_scatter.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,15 +3,19 @@
from LabIFSC2 import *

campo_magnético=arrayM([250,150,110,90,70,60,55,40,25,20],1,'muT')
distancias=linspace(1,10,10,0.5,'cm')
distancias=linspaceM(1,10,10,0.5,'cm')

unidade_x='cm'
unidade_y='muT'

plt.style.use('ggplot')
plt.errorbar(nominais(distancias,unidade_x),nominais(campo_magnético,unidade_y),
xerr=incertezas(distancias,unidade_x),yerr=incertezas(campo_magnético,unidade_y),fmt='o')
plt.errorbar(x=nominais(distancias,unidade_x),
y=nominais(campo_magnético,unidade_y),
xerr=incertezas(distancias,unidade_x),
yerr=incertezas(campo_magnético,unidade_y),
fmt='o')

plt.xlabel(f"Distancia ({unidade_x})")
plt.ylabel(f"Campo magnético ({unidade_y})")
#plt.savefig("exemplo.jpg")
plt.savefig('docs/images/graficos_scatter.jpg',dpi=300)
plt.cla()
Loading

0 comments on commit a9076ae

Please sign in to comment.