Полиномиальный фитинг, проходящий через 1 точку с производной силы =0

Я ищу что-то похожее на numpy.polyfit У меня есть кривая с фиксированной точкой, которая выглядит как полином степени 2

Что я хочу сделать, так это:

  • прохождение точно через первую точку (в следующем примере (0.05 , 1.0))
  • По первому пункту иметь derivative =0

пример :

TabX :

0,050 ; 0,055 ; 0,060 ; 0,065 ; 0,070 ; 0,075 ; 0,080 ; 0,085 ; 0,090 ; 0,095 ;     0,100 ; 0,110 ; 0,120 ; 0,130 ; 0,140 ; 0,150 ; 0,160 ; 0,170 ; 0,180 ; 0,190 ; 0,200 ; 0,210 ; 0,220 ; 0,230 ; 0,243

TabY :

1,000000000 ; 1,000446069 ; 1,000395689 ; 1,000359466 ; 1,000313867 ; 0,999937484 ; 0,999760969 ; 0,999811533 ; 0,999966352 ; 0,999767956 ; 1,000148567 ; 1,000634904 ; 1,000735849 ; 1,001199937 ; 1,001510678 ; 1,001722496 ; 1,001992602 ; 1,002487029 ; 1,003492247 ; 1,004006533 ; 1,004832845 ; 1,005730132 ; 1,006327527 ; 1,007109894 ; 1,008266254

Я уже нашел «простое, но варварское» решение для прохождения первой точки: я добавляю к этой точке много веса либо с помощью функции веса, либо добавляя много 0.05 в TabX и много 1.0 в TabY и используя обычная функция np.polyfit. Это некрасиво, но это работает.

Но я действительно не знаю, как было derivative=0 в точке (0.05;1.0)

также в этой ветке говорится, что Лагранж множитель мог бы помочь, но я не смог использовать функцию, написанную @Jaime, у меня возникла ошибка LinAlgError, 'Singular matrix'

Кроме того, поскольку я должен иметь возможность запускать этот скрипт в базовом Abaqus, решение может использовать только базовые python 2.7 и numpy. В этой реализации нет возможности использовать scipy или matplotlib.

_____________________________________________________________________

Редактировать: используя ответ DanielF, я мог бы сделать что-то наполовину работающее (также спасибо DanielF за исправление в моем исходном сообщении, которое намного проще читать)

Однако использование нового источника - хорошая идея, я не могу сделать его эффективным с моими данными, это работает:

def WorkingDeg4 ():
    x = np.arange(100)
    y0 = 4.0*x**4+0.07 * x ** 3 + 0.3 * x ** 2 + 1.1 * x
    y = y0 + 1000 * np.random.randn(x.shape[0])
    XX = np.vstack((x**4,x ** 3, x ** 2, x, np.ones_like(x))).T
    p_all = np.linalg.lstsq(XX, y)[0]
    pp = np.polyfit(x, y, 3)
    p_no_offset = np.linalg.lstsq(XX[:, :-1], y)[0]
    y_fit = np.dot(p_no_offset, XX[:, :-1].T)
    for i in range(0,len(x)):
        print x[i],y0[i],y[i],y_fit[i]

Но если я хочу сделать int, используя мои данные в основном, я поставлю:

if __name__ == '__main__':
    MyDataX=[0.050 , 0.055 , 0.060 , [...], 0.243]
    MyDataY=[1.000000000 , 1.000446069 , 1.000395689 , [...] , 1.008266254]
    TabX=[0.0]*len(MyDataX)
    TabY=[0.0]*len(MyDataY)
    for i in range(0,len(TabX)):
        TabX[i]=MyDataX[i]-MyDataX[0]
        TabY[i]=MyDataY[i]-MyDataY[0]

Итак, на этом этапе я выполнил фазу «возврата к исходному состоянию».

И я хочу сделать то же самое, что и def Working, но для своих данных, поэтому я сделал копию tge WorkingDeg4 и просто избавился от создания x и y и поместил его в аргумент.

def NOTWorkingDeg4 (x,y): 
    XX = np.vstack((x**4,x ** 3, x ** 2, x, np.ones_like(x))).T
    p_all = np.linalg.lstsq(XX, y)[0]
    pp = np.polyfit(x, y, 3)
    p_no_offset = np.linalg.lstsq(XX[:, :-1], y)[0]
    y_fit = np.dot(p_no_offset, XX[:, :-1].T)

а этот не работает.... у меня мистейк на линии

XX = np.vstack((x**4,x ** 3, x ** 2, x, np.ones_like(x))).T

говорю TypeError: unsupported operand type(s) for ** or pow(): 'list' and 'int'

Итак, насколько я понимаю, он не хочет этого делать, потому что когда он делает x**4 , x не является целым числом. Но я не знаю, как обойти проблему

____________________________________________________________________________ edit 2: нашел это: проблема заключалась в том, чтобы инициировать TabX и TabY не как массив, а как np.array неправильно:

TabX=[0.0]*len(MyDataX)
TabY=[0.0]*len(MyDataY)

Верно :

TabX=np.array([0.0]*len(LongueursFissureGlobale))
TabY=np.array([0.0]*len(CourbeInterpolationGlobale))

person texas    schedule 31.08.2017    source источник
comment
Ааа abaqus. Я знаю эту боль.   -  person Daniel F    schedule 31.08.2017
comment
Вы, вероятно, захотите сдвинуть свои данные так, чтобы (0.5, 1.0) было вашим источником (0, 0), чтобы облегчить ваши вычисления, а затем выполните numpy.linalg.lstsq. См. этот ответ. Затем вы можете установить первую производную равной нулю, отбрасывая не только константы из матрицы a, но и значения x.   -  person Daniel F    schedule 31.08.2017
comment
если вы посмотрите, как вычисляется минимизация такого рода, то будет довольно легко вывести модифицированную формулу. Через несколько часов что-нибудь напишу.   -  person Mad Physicist    schedule 31.08.2017


Ответы (1)


Вы, вероятно, захотите сдвинуть свои данные так, чтобы (0.05, 1.0) было вашим источником (0, 0), чтобы облегчить ваши расчеты, а затем выполните numpy.linalg.lstsq

x = TabX - 0.05
y = TabY - 1.0
X_poly = np.vstack((x ** 4, x ** 3, x ** 2))
poly_coeffs = np.linalg.lstsq(X_poly.T, y)
y_fit = np.dot(poly_coeffs, X_poly)

Если это необходимо, вам придется преобразовать полином обратно в ваши старые координаты, но это преобразование значительно упрощает подгонку.

Подробнее см. в этом ответе.

person Daniel F    schedule 31.08.2017
comment
Да, когда имеешь дело с полиномиальными производными, обнуление как можно большего количества вещей значительно упрощает вычисления. Если это ответ на ваш вопрос, не забудьте поставить галочку. - person Daniel F; 01.09.2017
comment
Это был один шаг, но мне нужно было время, чтобы отредактировать первый пост, так как я еще не мог его решить. Но спасибо за вашу помощь. - person texas; 01.09.2017