Estou tentando ajustar alguns dados a uma lei de energia usando python. O problema é que alguns dos meus pontos são limites superiores, que não sei como incluir na rotina de adaptação.

Nos dados, eu coloquei os limites superiores como erros em y iguais a 1, quando o restante é muito menor. Você pode colocar esses erros em 0 e alterar o gerador de lista de uplims, mas o ajuste é terrível.

O código é o seguinte:

import numpy as np
import matplotlib.pyplot as plt
from scipy.odr import *

# Initiate some data
x = [1.73e-04, 5.21e-04, 1.57e-03, 4.71e-03, 1.41e-02, 4.25e-02, 1.28e-01, 3.84e-01, 1.15e+00]
x_err = [1e-04, 1e-04, 1e-03, 1e-03, 1e-02, 1e-02, 1e-01, 1e-01, 1e-01]
y = [1.26e-05, 8.48e-07, 2.09e-08, 4.11e-09, 8.22e-10, 2.61e-10, 4.46e-11, 1.02e-11, 3.98e-12]
y_err = [1, 1, 2.06e-08, 2.5e-09, 5.21e-10, 1.38e-10, 3.21e-11, 1, 1]

# Define upper limits
uplims = np.ones(len(y_err),dtype='bool')
for i in range(len(y_err)):
    if y_err[i]<1:
        uplims[i]=0
    else:
        uplims[i]=1

# Define a function (power law in our case) to fit the data with.
def function(p, x):
     m, c = p
     return m*x**(-c)

# Create a model for fitting.
model = Model(function)

# Create a RealData object using our initiated data from above.
data = RealData(x, y, sx=x_err, sy=y_err)

# Set up ODR with the model and data.
odr = ODR(data, model, beta0=[1e-09, 2])
odr.set_job(fit_type=0)   # 0 is full ODR and 2 is least squares; AFAIK, it doesn't change within errors
# more details in https://docs.scipy.org/doc/scipy/reference/generated/scipy.odr.ODR.set_job.html

# Run the regression.
out = odr.run()


# Use the in-built pprint method to give us results.
#out.pprint()   #this prints much information, but normally we don't need it, just the parameters and errors; the residual variation is the reduced chi square estimator

print('amplitude = %5.2e +/- %5.2e \nindex = %5.2f +/- %5.2f \nchi square = %12.8f'% (out.beta[0], out.sd_beta[0], out.beta[1], out.sd_beta[1], out.res_var))

# Generate fitted data.
x_fit = np.linspace(x[0], x[-1], 1000)    #to do the fit only within the x interval; we can always extrapolate it, of course
y_fit = function(out.beta, x_fit)


# Generate a plot to show the data, errors, and fit.
fig, ax = plt.subplots()

ax.errorbar(x, y, xerr=x_err, yerr=y_err, uplims=uplims, linestyle='None', marker='x')
ax.loglog(x_fit, y_fit)
ax.set_xlabel(r'$x
O resultado do ajuste é:
amplitude = 3.42e-12 +/- 5.32e-13
index =  1.33 +/-  0.04
chi square =   0.01484021
Como você pode ver na plotagem, os dois primeiros e os dois últimos pontos são limites superiores e o ajuste não os leva em consideração. Além disso, no penúltimo ponto, o ajuste é ultrapassado, mesmo que isso fosse estritamente proibido.
Eu preciso que o ajuste saiba que esses limites são muito rígidos, e não tente ajustar o ponto em si, mas considere-os apenas como limites. Como eu poderia fazer isso com a rotina odr (ou qualquer outro código que me ajusta e me dá um estimador de qui-quadrado)?
Por favor, leve em consideração que preciso alterar a função para outras generalizações facilmente, para que coisas como o módulo powerlaw não sejam desejáveis.
Obrigado!
)
ax.set_ylabel(r'$f(x) = m·x^{-c}
O resultado do ajuste é:
amplitude = 3.42e-12 +/- 5.32e-13
index =  1.33 +/-  0.04
chi square =   0.01484021
Como você pode ver na plotagem, os dois primeiros e os dois últimos pontos são limites superiores e o ajuste não os leva em consideração. Além disso, no penúltimo ponto, o ajuste é ultrapassado, mesmo que isso fosse estritamente proibido.
Eu preciso que o ajuste saiba que esses limites são muito rígidos, e não tente ajustar o ponto em si, mas considere-os apenas como limites. Como eu poderia fazer isso com a rotina odr (ou qualquer outro código que me ajusta e me dá um estimador de qui-quadrado)?
Por favor, leve em consideração que preciso alterar a função para outras generalizações facilmente, para que coisas como o módulo powerlaw não sejam desejáveis.
Obrigado!
)
ax.set_title('Power Law fit')

plt.show()

O resultado do ajuste é:

amplitude = 3.42e-12 +/- 5.32e-13
index =  1.33 +/-  0.04
chi square =   0.01484021

Como você pode ver na plotagem, os dois primeiros e os dois últimos pontos são limites superiores e o ajuste não os leva em consideração. Além disso, no penúltimo ponto, o ajuste é ultrapassado, mesmo que isso fosse estritamente proibido.

Eu preciso que o ajuste saiba que esses limites são muito rígidos, e não tente ajustar o ponto em si, mas considere-os apenas como limites. Como eu poderia fazer isso com a rotina odr (ou qualquer outro código que me ajusta e me dá um estimador de qui-quadrado)?

Por favor, leve em consideração que preciso alterar a função para outras generalizações facilmente, para que coisas como o módulo powerlaw não sejam desejáveis.

Obrigado!