Блог

У меня есть одномерный массив данных, который я пытаюсь смоделировать как гиперболу, используя три параметра. Я пытаюсь реализовать алгоритм Левенберга-Марквардта, используя функцию leastsq из библиотеки scipy.optimize. Однако моя программа застревает на итерации, где число делится на ноль, и я не понимаю почему.

Некоторая справочная информация: данные одномерного массива в основном представляют собой значения лакунарности для блоков разного размера. Я сгенерировал данные о лакунарности из некоторых звуковых файлов, контекст к которым можно найти здесь.

В алгоритме функция наименьших квадратов принимает три входных данных:
(а) первоначальное предположение для трех параметров
(б) координата x для задачи наименьших квадратов — это в основном одномерный массив целых чисел от 1 до 100 в моей задаче
(в) координата y для задачи наименьших квадратов — это одномерный массив, в котором хранятся значения лакунарности. Итак, значения лакунарности являются функцией x, где x варьируется от 1 до 100.

Гипербола моделируется с использованием трех параметров a, b и c как

Код выдает следующую ошибку:
«OverflowError: не удается преобразовать бесконечность с плавающей точкой в целое число»

Код:

 #import
from scipy import *
from scipy.optimize import leastsq
import matplotlib.pylab as plt
import numpy as np
import codecs, json
from math import *

# Define your function to calculate the residuals. 
#The fitting function holds your parameter values.  
def residuals(p, y, x):
        err = y-pval(x,p)
        return err

def pval(x, p):
        z = x
        for i in range(100):
                print(x)
                print(x[i]**p[1])
                z[i] = p[0]/(x[i]**p[1]) p[2]
        return z

#read in your data
obj_text = codecs.open('textfilesCC1.json', 'r', encoding='utf-8').read()
b_new = json.loads(obj_text)
data = np.array(b_new)
x = np.arange(1,101)
y = data[1:101]

#guess at initial parameters
A1_0=1.0
A2_0=1.0
A3_0=0.5

#leastsq package calls the Levenberg-Marquardt algorithm
pname = (['A1','A2','A3'])
p0 = array([A1_0 , A2_0, A3_0])
plsq = leastsq(residuals, p0, args=(y, x), maxfev=2000)

# Now, plot your data
plt.plot(x,y,'xo',x,pval(x,plsq[0]),'x')
title('Least-squares fit to data')
xlabel('x')
ylabel('y')
legend(['Data', 'Fit'],loc=4)

# Your best-fit paramters are kept within plsq[0].
print(plsq[0])
  

Согласно ошибке, значение x изменяется на 0 в какой-то момент итерации, а первый параметр a в конечном итоге делится на ноль, что приводит к ошибке.

Для устранения неполадок я напечатал значения x [i] ^ b и массив x во время выполнения кода, и вы можете увидеть значения здесь. Я вижу, что массив x модифицируется, чего не должно происходить. x должен оставаться одномерным массивом натуральных чисел от 1 до 100 и не изменяться на итерации. Я не смог определить, где именно находится код, модифицирующий массив x.

Я ожидаю, что массив x останется неизменным, а код выведет последние три значения параметров a, b и c.

РЕДАКТИРОВАТЬ: Я внес некоторые изменения в свой код, после чего он успешно заработал. Ниже приведены те изменения, которые могут кому-либо быть интересны:

1. Не определил z как z = x, а просто определил его как z = np.диапазон (1,101). Результатом стало то, что массив x больше не изменился, чего и ожидалось.

2. Изменен тип данных массивов x и y на float с использованием

 x = np.array(x, dtype=np.float64)
  

3. Я снова застрял на фрагменте кода, который выводит данные на график. Я получил сообщение об ошибках «‘title’ не определен. Аналогичные ошибки для xlabel, ylabel. Итак, я просто удалил эти строки и просто застрял с

 plt.plot(x,y,'red',x,pval(x,plsq[0]),'blue')
plt.show()
  

Это не прямой ответ на ваш вопрос, но поскольку вы используете возведение в степень ( ** ), я настоятельно рекомендую вам заранее преобразовать все ваши числа в Decimal , чтобы избежать потери точности, присущей арифметике с плавающей запятой при больших значениях.

Например:

 import decimal
decimal.getcontext().prec = 100

A1_0=Decimal("1.0")
A2_0=Decimal("1.0")
A3_0=Decimal("0.5")

x = [Decimal(f) for f in x]
y = [Decimal(f) for f in y]
  

Возможно, ваш ноль «окажется» небольшим значением, близким к нулю…