0

I have a function as the following

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

I have experimental answers of above equation and once I must to put some data into above function and solve equation below

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 and q_exp are in the data file(2.txt). Now I have to choose a range for o_m (0.2 to 0.4) and find in what o_m, the chi function will be minimized.

my code is:

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] are z1, z2, q_exp and err, respectively in the data file. and z1 and z2 are the integration range. I need your help and I appreciate your time and your attention. Please do not rate a negative value. I need your answers.

mikuszefski
  • 3,943
  • 1
  • 25
  • 38
Ethan
  • 180
  • 2
  • 15
  • 1: What is your problem? 2: Please share some minimal data set. 3: why `ant()` has `o_m` and `o_d` while `q` above has only `o_m`. – mikuszefski Aug 04 '17 at 09:52
  • Why not use `scipy.optimize.leastsq`? BTW, if the data is not too large, just load it once in the beginning, maybe with `numpy.loadtxt()` – mikuszefski Aug 04 '17 at 09:58
  • Moreover, there are some typos; in the last `print()` you open with `'` but close with `"` and you probably want to print `c` not `q`, which you should name, e.g. `c2` as it is actually the square. This avoids confusion. Some indents seem wrong. Comment on typing pythonic: `min == None` works but `min is None` looks better. Maybe even `if not min`. And do you really want to compare to `h` or to `c`? – mikuszefski Aug 04 '17 at 10:08
  • it is chi square test. we have an integration in which a data file will be used to reach final answer. we must minimizing the chi square and find what parameter in what value minimized the chi square. om and od are omega d and omega m in the main function in print we only need om or od. – Ethan Aug 04 '17 at 10:55
  • I corrected some of parts that you mentioned – Ethan Aug 04 '17 at 10:58
  • OK, but 1 and 2 from my first comment remain. What are you actually asking. That should be somewhat improved to avoid downvotes. – mikuszefski Aug 04 '17 at 12:36
  • Another important question, is the interval `[0.2,0.4]` a known boundary for the fit or do you just assume it lies inside this interval. Are values outside the interval allowed? – mikuszefski Aug 04 '17 at 12:43

2 Answers2

1

Here is my understanding of the problem. First I generate some data by the following code

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)

which I would then fit in the following way

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()

giving output

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

and the graph Fit Note that for leastsq the covariance matrix is in reduced form and needs to be rescaled.

mikuszefski
  • 3,943
  • 1
  • 25
  • 38
  • Thnk you very much, it is perfect. But fitting data is not what i expect. We change om and we check the out puts that if there is any om makes the chi function minimum. My code is correct. But needs some revisions – Ethan Aug 05 '17 at 13:27
  • @Ehsan The original version I wouldn't call "correct"; anyhow, that's why I asked above what it is that you actually want. You show a one-dimensional fit in an inefficient way. As you have `o_m` and `o_d` I presented the general 2D case, which you can easily modify to 1D. If you want a version that is working, you should really post a minimum data set, maybe 20 to 30 points. Apart from that my code is doing exactly what your code is doing except for the inefficient `for` loop you are using. From that point of view, it is a revision of your code, like revision 10 leaving out steps 2 to 9. – mikuszefski Aug 07 '17 at 06:11
  • @Ehsan I really recommend to really answer the questions I put in the comments above. – mikuszefski Aug 07 '17 at 06:14
  • thank you so much. I read it over and over from ###fitting til end, but I could not quite understand what them are. and what is the relation between them and the last section of my question. In fact i revised my code again but the problem is: I must find a "o_m" that results in minimum "chi" or minimum "c". I can find it but in wrong way. in each range that i choose for o_m in for loop, the minimum value of that loop results in minimum value of "c". for example in my code we have (20,40,1) for o_m and 20results in minimum. it is wrong and i do not know how to fix it – Ethan Aug 10 '17 at 08:25
  • because minimum must be between 25 to 31 or a value which is not the minimum number or maximum number or our arbitrary range in loop. – Ethan Aug 10 '17 at 08:27
  • My friend @mikuszefski, maybe this question help you to have another sight of view of this question. I appreciate your help and your attention my friend https://stackoverflow.com/q/45631732/8339406 – Ethan Aug 11 '17 at 09:46
0

Unconcsiosely, this question overlap my other question. The correct answer is:

 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)
Ethan
  • 180
  • 2
  • 15