Итак, здесь происходит несколько вещей, которые, к сожалению, мешают тому, как вы делаете вещи.Взгляните на этот код:
for i,j in zip(x,y):
b += i
c += i**2
d += i**3
e += i**4
m += j
n += j*i
p += j*i**2
Вы строите объекты так, что значения x не только возводятся в квадрат, но и в кубы, а четвертое питание.
Если вы распечатываете каждый изэти значения перед тем, как поместить их в матрицу 3 x 3 для решения:
In [35]: a = b = c = d = e = m = n = p = 0
...: a = len(x)
...: for i,j in zip(xx,y):
...: b += i
...: c += i**2
...: d += i**3
...: e += i**4
...: m += j
...: n += j*i
...: p += j*i**2
...: print(a, b, c, d, e, m, n, p)
...:
...:
3 18.744836 117.12356813829001 731.8283056811686 4572.738547313946 0.9294744720000001 5.807505391292503 36.28641270376207
При работе с арифметикой с плавающей точкой и особенно для небольших значений порядок операций имеет значение.Здесь происходит то, что, по счастливой случайности, сочетание как малых значений, так и больших значений, которые были вычислены, приводит к значению, которое очень мало.Поэтому, когда вы вычисляете детерминант с использованием факторизованной формы и расширенной формы, обратите внимание на то, как вы получите немного отличающиеся результаты, но также посмотрите на точность значений:
In [36]: det = determinant(a,b,c,b,c,d,c,d,e)
In [37]: det
Out[37]: 1.0913403514223319e-10
In [38]: det = determinantAlt(a,b,c,b,c,d,c,d,e)
In [39]: det
Out[39]: 2.3283064365386963e-10
Определитель имеет порядок 10 -10 !Причина расхождения заключается в том, что с помощью арифметики с плавающей точкой теоретически оба метода детерминантов должны давать один и тот же результат, но, к сожалению, в действительности они дают немного разные результаты, и это происходит из-за того, что называется распространение ошибки .Поскольку существует конечное число битов, которые могут представлять число с плавающей запятой, порядок операций изменяет способ распространения ошибки, поэтому, даже если вы удаляете скобки и формулы по существу совпадают, порядок операций, чтобы добраться доРезультат сейчас другой.Эта статья очень полезна для любого разработчика программного обеспечения, который регулярно занимается арифметикой с плавающей запятой: Что должен знать каждый компьютерный специалист об арифметике с плавающей запятой .
Поэтому, когда вы пытаетесьчтобы решить систему с помощью правила Крамера, неизбежно, когда вы делите на основной детерминант в вашем коде, даже если изменение порядка 10 -10 , изменение между этими двумя методами незначительно, но вы будетеПолучите очень разные результаты, потому что вы делите на это число при решении для коэффициентов.
Причина, по которой у NumPy нет этой проблемы, состоит в том, что они решают систему методом наименьших квадратов и псевдо -обратный и без использования правила Крамера.Я бы не рекомендовал использовать правило Крамера для нахождения коэффициентов регрессии, в основном из-за опыта и более надежных способов сделать это.
Однако, чтобы решить вашу конкретную проблему, хорошо бы нормализовать данные таким образом, что динамический диапазон теперь центрирован на 0. Следовательно, функции, которые вы используете для построения матрицы коэффициентов, являются более разумными, и, таким образом, вычислительному процессу легче обрабатывать данные.В вашем случае должно работать что-то простое: вычитать данные со средним значением x.Таким образом, если у вас есть новые точки данных, которые вы хотите предсказать, вы должны вычесть среднее значение данных x, прежде чем делать прогноз.
Поэтому в началеваш код, выполнить среднее вычитание и регрессировать на этих данных.Я показал вам, где я изменил код, указанный выше:
import numpy as np
x = [6.230825,6.248279,6.265732]
y = [0.312949,0.309886,0.306639472]
# Calculate mean
me = sum(x) / len(x)
# Make new dataset that is mean subtracted
xx = [pt - me for pt in x]
#toCheck = x[2]
# Data point to check is now mean subtracted
toCheck = x[2] - me
def evaluateValue(coeff,x):
c,b,a = coeff
val = np.around( a+b*x+c*x**2,9)
act = 0.306639472
error= np.abs(act-val)*100/act
print("Value = {:.9f} Error = {:.2f}%".format(val,error))
###### USing numpy######################
coeff = np.polyfit(xx,y,2) # Change
evaluateValue(coeff, toCheck)
################# Using explicit formula
def determinant(a,b,c,d,e,f,g,h,i):
# the matrix is [[a,b,c],[d,e,f],[g,h,i]]
return a*(e*i - f*h) - b*(d*i - g*f) + c*(d*h - e*g)
a = b = c = d = e = m = n = p = 0
a = len(x)
for i,j in zip(xx,y): # Change
b += i
c += i**2
d += i**3
e += i**4
m += j
n += j*i
p += j*i**2
det = determinant(a,b,c,b,c,d,c,d,e)
c0 = determinant(m,b,c,n,c,d,p,d,e)/det
c1 = determinant(a,m,c,b,n,d,c,p,e)/det
c2 = determinant(a,b,m,b,c,n,c,d,p)/det
evaluateValue([c2,c1,c0], toCheck)
######Using another explicit alternative
def determinantAlt(a,b,c,d,e,f,g,h,i):
return a*e*i - a*f*h - b*d*i +b*g*f + c*d*h - c*e*g # <- barckets removed
a = b = c = d = e = m = n = p = 0
a = len(x)
for i,j in zip(xx,y): # Change
b += i
c += i**2
d += i**3
e += i**4
m += j
n += j*i
p += j*i**2
det = determinantAlt(a,b,c,b,c,d,c,d,e)
c0 = determinantAlt(m,b,c,n,c,d,p,d,e)/det
c1 = determinantAlt(a,m,c,b,n,d,c,p,e)/det
c2 = determinantAlt(a,b,m,b,c,n,c,d,p)/det
evaluateValue([c2,c1,c0], toCheck)
Когда я запустил это, мы теперь получим:
In [41]: run xor
Value = 0.306639472 Error = 0.00%
Value = 0.306639472 Error = 0.00%
Value = 0.306639472 Error = 0.00%
В качестве окончательного прочтения для вас, это аналогичная проблема, с которой столкнулся кто-то другой, с которой я обратился в своем вопросе: Подгонка квадратичной функции в python без nimpy polyfit .В итоге я посоветовал им не использовать правило Крамера и использовать наименьшие квадраты через псевдообратное.Я показал им, как получить точно такие же результаты, не используя numpy.polyfit.Кроме того, использование метода наименьших квадратов позволяет определить, где, если у вас более 3 точек, вы все равно можете уместить квадратик через свои точки, чтобы у модели была наименьшая возможная ошибка.