Программа для подгонки гиперболы к линейным данным с использованием метода наименьших квадратов (алгоритм Левенберга-Марквардта) не работает должным образом - PullRequest
0 голосов
/ 07 апреля 2019

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

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

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

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

enter image description here

Код выдает следующую ошибку:
«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('textfiles\CC1.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.arange (1,101).В результате массив x больше не изменился, что и ожидалось.

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

x = np.array(x, dtype=np.float64)
Я снова застрял в куске кода, который отображает данные.Я получил ошибки "title" не определен. Подобные ошибки для xlabel, ylabel. Поэтому я просто удалил эти строки и просто застрял с
plt.plot(x,y,'red',x,pval(x,plsq[0]),'blue')
plt.show()

1 Ответ

0 голосов
/ 07 апреля 2019

Не является прямым ответом на ваш вопрос, но, поскольку вы используете возведение в степень (**), я настоятельно рекомендую заранее преобразовать все ваши числа в 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]

Возможно, ваш ноль "превратится" в небольшое значение, близкое к нулю ...

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...