Минимизация функции с использованием python для подгонки данных

У меня есть следующая функция

q = 1 / sqrt( ((1+z)**2 * (1+0.01*o_m*z) - z*(2+z)*(1-o_m)) )
h = 5 * log10( (1+z)*q ) + 43.1601

У меня есть экспериментальные ответы на вышеуказанное уравнение, и однажды я должен поместить некоторые данные в вышеуказанную функцию и решить уравнение ниже

chi=(q_exp-q_theo)**2/err**2  # this function is a sigma, sigma chi from z=0 to z=1.4 (in the data file)

z, err и q_exp находятся в файле данных (2.txt). Теперь мне нужно выбрать диапазон для o_m (0.2 to 0.4) и найти, в каком o_m функция chi будет свернута.

мой код:

from math import *
from scipy.integrate import quad

min = None
l = None
a = None
b = None
c = 0


def ant(z,om,od):
    return 1/sqrt( (1+z)**2 * (1+0.01*o_m*z) - z*(2+z)*o_d )


for o_m in range(20,40,1):
    o_d=1-0.01*o_m
    with open('2.txt') as fp:
        for line in fp:
            n = list( map(float, line.split()) )
            q = quad(ant,n[0],n[1],args=(o_m,o_d))[0]
            h = 5.0 * log10( (1+n[1])*q ) + 43.1601
            chi = (n[2]-h)**2 / n[3]**2
            c = c + chi
        if min is None or min>c:
            min = c
            l = o_m                            
        print('chi=',q,'o_m=',0.01*l)

n[1], n[2], n[3], n[4] — это z1, z2, q_exp и err соответственно в файле данных. а z1 и z2 — диапазон интегрирования. Мне нужна ваша помощь, и я ценю ваше время и ваше внимание. Пожалуйста, не оценивайте отрицательное значение. Мне нужны ваши ответы.


person Ethan    schedule 03.08.2017    source источник
comment
1: В чем твоя проблема? 2: Пожалуйста, поделитесь минимальным набором данных. 3: почему у ant() есть o_m и o_d, а у q выше есть только o_m.   -  person mikuszefski    schedule 04.08.2017
comment
Почему бы не использовать scipy.optimize.leastsq? Кстати, если данные не слишком велики, просто загрузите их один раз в начале, возможно, с помощью numpy.loadtxt()   -  person mikuszefski    schedule 04.08.2017
comment
Кроме того, есть некоторые опечатки; в последнем print() вы открываете с помощью ', но закрываете с помощью ", и вы, вероятно, хотите напечатать c, а не q, которое вы должны назвать, например. c2 так как на самом деле это квадрат. Это позволяет избежать путаницы. Некоторые отступы кажутся неправильными. Комментарий к вводу pythonic: min == None работает, но min is None выглядит лучше. Может быть, даже if not min. И вы действительно хотите сравнить с h или с c?   -  person mikuszefski    schedule 04.08.2017
comment
это тест хи-квадрат. у нас есть интеграция, в которой файл данных будет использоваться для получения окончательного ответа. мы должны минимизировать хи-квадрат и найти, какой параметр в каком значении минимизирует хи-квадрат. om и od — это омега d и омега m в функции main, в print нам нужны только om или od.   -  person Ethan    schedule 04.08.2017
comment
Я исправил некоторые части, которые вы упомянули   -  person Ethan    schedule 04.08.2017
comment
Хорошо, но 1 и 2 из моего первого комментария остались. Что вы на самом деле спрашиваете. Это должно быть несколько улучшено, чтобы избежать отрицательных голосов.   -  person mikuszefski    schedule 04.08.2017
comment
Другой важный вопрос: является ли интервал [0.2,0.4] известной границей подгонки или вы просто предполагаете, что он находится внутри этого интервала. Допустимы ли значения вне интервала?   -  person mikuszefski    schedule 04.08.2017


Ответы (2)


Вот мое понимание проблемы. Сначала я генерирую некоторые данные с помощью следующего кода

import numpy as np
from scipy.integrate import quad
from random import random


def boxmuller(x0,sigma):
    u1=random()
    u2=random()
    ll=np.sqrt(-2*np.log(u1))
    z0=ll*np.cos(2*np.pi*u2)
    z1=ll*np.cos(2*np.pi*u2)
    return sigma*z0+x0, sigma*z1+x0


def q_func(z, oM, oD):
    den= np.sqrt( (1.0 + z)**2 * (1+0.01 * oM * z) - z * (2+z) * (1-oD) )
    return 1.0/den


def h_func(z,q): 
    out = 5 * np.log10( (1.0 + z) * q ) + .25#43.1601
    return out


def q_Int(z1,z2,oM,oD):
    out=quad(q_func, z1,z2,args=(oM,oD))
    return out

ooMM=0.3
ooDD=1.0-ooMM

dataList=[]
for z in np.linspace(.3,20,60):
    z1=.1+.1*z*.01*z**2
    z2=z1+3.0+.08+z**2
    q=q_Int(z1,z2,ooMM,ooDD)[0]
    h=h_func(z,q)
    sigma=np.fabs(.01*h)
    h=boxmuller(h,sigma)[0]
    dataList+=[[z,z1,z2,h,sigma]]
dataList=np.array(dataList)

np.savetxt("data.txt",dataList)

который я бы тогда подогнал следующим образом

import matplotlib
matplotlib.use('Qt5Agg')

from matplotlib import pyplot as plt
import numpy as np
from scipy.integrate import quad
from scipy.optimize import leastsq 

def q_func(z, oM, oD):
    den= np.sqrt( (1.0 + z)**2 * (1+0.01 * oM * z) - z * (2+z) * (1-oD) )
    return 1.0/den


def h_func(z,q): 
    out = 5 * np.log10( (1.0 + z) * q ) + .25#43.1601
    return out


def q_Int(z1,z2,oM,oD):
    out=quad(q_func, z1,z2,args=(oM,oD))
    return out


def residuals(parameters,data):
    om,od=parameters
    zList=data[:,0]
    yList=data[:,3]
    errList=data[:,4]
    qList=np.fromiter( (q_Int(z1,z2, om,od)[0] for  z1,z2 in data[ :,[1,2] ]), np.float)
    hList=np.fromiter( (h_func(z,q) for z,q in zip(zList,qList)), np.float)
    diffList=np.fromiter( ( (y-h)/e for y,h,e in zip(yList,hList,errList) ), np.float)
    return diffList

dataList=np.loadtxt("data.txt")

###fitting 
startGuess=[.4,.8]
bestFitValues, cov,info,mesg, ier = leastsq(residuals, startGuess , args=( dataList,),full_output=1)
print bestFitValues,cov

fig=plt.figure()
ax=fig.add_subplot(1,1,1)
ax.plot(dataList[:,0],dataList[:,3],marker='x')


###fitresult
fqList=[q_Int(z1,z2,bestFitValues[0], bestFitValues[1])[0] for z1,z2 in zip(dataList[:,1],dataList[:,2])]
fhList=[h_func(z,q) for z,q in zip(dataList[:,0],fqList)]
ax.plot(dataList[:,0],fhList,marker='+')


plt.show()

выход

>>[ 0.31703574  0.69572673] 
>>[[  1.38135263e-03  -2.06088258e-04]
>> [ -2.06088258e-04   7.33485166e-05]]

и график Fit Обратите внимание, что для leastsq ковариационная матрица находится в сокращенной форме и требует масштабирования.

person mikuszefski    schedule 04.08.2017
comment
Большое спасибо, это идеально. Но подходящие данные - это не то, что я ожидаю. Мы меняем ом и проверяем вывод, что если ом есть, то функция ци минимальна. Мой код правильный. Но требует некоторых доработок - person Ethan; 05.08.2017
comment
@Ehsan Оригинальную версию я бы не назвал правильной; во всяком случае, именно поэтому я спросил выше, чего вы на самом деле хотите. Вы показываете одномерную подгонку неэффективным способом. Поскольку у вас есть o_m и o_d, я представил общий случай 2D, который вы можете легко преобразовать в 1D. Если вам нужна работающая версия, вам следует опубликовать минимальный набор данных, возможно, от 20 до 30 пунктов. Кроме того, мой код делает именно то, что делает ваш код, за исключением неэффективного цикла for, который вы используете. С этой точки зрения, это ревизия вашего кода, например, редакция 10, в которой отсутствуют шаги со 2 по 9. - person mikuszefski; 07.08.2017
comment
@Ehsan Я действительно рекомендую действительно ответить на вопросы, которые я задал в комментариях выше. - person mikuszefski; 07.08.2017
comment
Спасибо большое. Я перечитывал это снова и снова от ###fitting до конца, но так и не смог понять, что это такое. и какова связь между ними и последним разделом моего вопроса. На самом деле я снова пересмотрел свой код, но проблема в том, что я должен найти o_m, который приводит к минимальному chi или минимальному c. Я могу найти его, но неправильно. в каждом диапазоне, который я выбираю для o_m в цикле for, минимальное значение этого цикла приводит к минимальному значению c. например, в моем коде у нас есть (20,40,1) для o_m и минимум 20 результатов. это неправильно, и я не знаю, как это исправить - person Ethan; 10.08.2017
comment
потому что минимум должен быть между 25 и 31 или значением, которое не является минимальным или максимальным числом или нашим произвольным диапазоном в цикле. - person Ethan; 10.08.2017
comment
Мой друг @mikuszefski, может быть, этот вопрос поможет вам по-другому взглянуть на этот вопрос. Я ценю вашу помощь и ваше внимание, мой друг stackoverflow.com/q/45631732/8339406 - person Ethan; 11.08.2017

Бессознательно, этот вопрос перекрывает мой другой вопрос. Правильный ответ:

 from math import *
 import numpy as np
 from scipy.integrate import quad
 min=l=a=b=chi=None
 c=0
 z,mo,err=np.genfromtxt('Union2.1_z_dm_err.txt',unpack=True)
 def ant(z,o_m):            #0.01*o_m  is steps of o_m
     return 1/sqrt(((1+z)**2*(1+0.01*o_m*z)-z*(2+z)*(1-0.01*o_m)))
 for o_m in range(20,40):
     c=0
     for i in range(len(z)):
         q=quad(ant,0,z[i],args=(o_m,))[0]     #Integration o to z
         h=5*log10((1+z[i])*(299000/70)*q)+25     #function of dL
         chi=(mo[i]-h)**2/err[i]**2               #chi^2 test function
        c=c+chi
        l=o_m
        print('chi^2=',c,'Om=',0.01*l,'OD=',1-0.01*l)
person Community    schedule 14.08.2017