Curve fitting in python

Question

Hey, I have a set of values for frequency and power spectrum and I have to plot Power spectrum Versus frequency on log scale. Once done, I need to pass the best fit straight line through it.. I get the line on a linear scale.. but when I try to superimpose it onto the freq-power spectrum plot, the resultant plot does not show any line, instead the data points of 1st plot are merely shifted in space. Also, the same line, if plotted on log scale using loglog function, does not show up.

Can somebody tell me what I should do in order to get the line on a Log scale?

SO I have a file having three columns; Frequency, Power spec. Power signal.. Here is a piece of what i wrote to plot the data and line..

#initialize all variables to 0

#open the data file

while 1:
  ln = datafile.readline()
  if ln:
    data = ln.split()
    x = float(n)
    y = float(data[0])
    z = float(data[1])
    xval.append(float(n))
    yval.append(y)
    zval.append(z)
    n += 1
    sum_z += z
    sum_y += y
    sum_y_squared += y*y
    sum_yz += y*z
  else:
    break
datafile.close()

# calculate slope and intercept using formulae
for num in xval:
    res = intercept + slope*num
    line.append(res)

#Plot data
pylab.figure(0)
matplotlib.pylab.loglog(yval,zval)

#Plot line
pylab.figure(0)
pylab.plotloglog(line)

Answer 1

As I understand your problem, you want to plot two lines to the same diagram. Here is how it is done in general:

import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(line1_x, line1_y)
ax.plot(line2_x, line2_y)
ax.set_yscale("log")

So, first you put them both in the same Axes, so they appear in the same diagram. TO modify the scaling, you can use set_xscale and set_yscale respectively.

Apart from that, I cannot help but notice that your code for reading the file is horrible. As @Bernhard suggests in his answer, try using numpy.loadtxt. This could look like this:

data = numpy.loadtxt("data.txt")
n = len(data)
x = numpy.arange(n)
sum_z = sum(data.T[1])
sum_y = sum(data.T[0])
sum_y_squared = sum(data.T[0]**2)
sum_yz = sum(data.T[0]*data.T[1])

This should give you the same results as your loop, only it is much more concise. I strongly recommend you read the Tentative NumPy Tutorial, as it explain a lot of the really cool features of numpy arrays.

Answer 2

Despite the fact that the plot line commands are not correct in your example I assume it is similar to what you actually do.

The second plot command plots on a different x range:

loglog(yval,zval) # plot  yval vs zval
loglog(line) #  plots range(0,len(line)) vs line

Also have you look at the values of line, do they make sense are they in the same range as yval, zval?

Additionally you might want to use numpy.loadtxt to load your data file.