# Plotting and Fitting

Hi, I have to learn to plot and fit functions (with error bars in x and y direction).

My question is: 1) Can I do it with sage? Would you prefer gnuplot?

How can I read and plot data like this: http://rn0.ru/show/1430/ easily in sage? 2) And if yes how?

I already tried to plot some randomly created data: http://rn0.ru/show/1429/ And that works great.

3) How can I fit a function (nonlinear) with considering the errorbars?

edit retag close merge delete

Sort by ยป oldest newest most voted

I'm not sure what you mean by "with considering the errorbars". Do you mean merely that you want the covariance of the fit (which leastsq returns), or that your data has errors in y (the usual model) and you want that to affect the fitting process, or that your data has errors in both x and y and you want that to affect the fit?

Anyway, here are two relatively straightforward ways to fit and get covariances.

sage: import numpy as np
sage: import scipy.optimize
sage: import scipy.odr
sage: import random
sage: random.seed(3)
sage:
sage: #
sage: # generate the data
sage: #
sage:
sage: x = np.linspace(0, 10, 100)
sage: y = np.sin(0.7*x)+2.4*x  # underlying function
sage: y += np.random.randn(len(x))/10 # add some noise
sage: x_err = [0.1]*len(x) # estimated errors
sage: y_err = [0.3]*len(x) # estimated errors
sage:
sage: # first, we can use curve_fit
sage:
sage: def fit1(x, a, b):
....:         return np.sin(x*a)+b*x
....:
sage: # note: can't handle errors in x
sage: popt, pcov = scipy.optimize.curve_fit(fit1, x, y, sigma=y_err)
sage: popt
array([ 0.70124584,  2.39768117])
sage: pcov
array([[  6.97634001e-06,  -2.58153187e-06],
[ -2.58153187e-06,   4.22745005e-06]])
sage:
sage: # or use the more powerful odr approach, which unfortunately has
sage: # some different conventions
sage:
sage: def fit2(B, x):
....:         return np.sin(x*B[0])+B[1]*x
....:
sage: mymodel = scipy.odr.Model(fit2)
sage: mydata = scipy.odr.RealData(x, y, sx=x_err, sy=y_err)
sage: myodr = scipy.odr.ODR(mydata, mymodel, beta0=[1., 1.], maxit=1000)
sage: myoutput = myodr.run()
sage: myoutput.pprint()
Beta: [ 0.70139136  2.39769177]
Beta Std Error: [ 0.00258362  0.00198363]
Beta Covariance: [[  9.08706871e-05  -2.35028948e-05]
[ -2.35028948e-05   5.35658971e-05]]
Residual Variance: 0.0734568446584
Inverse Condition #: 0.679929358499
Reason(s) for Halting:
Numerical error detected

more

Thanks for providing this example :) I think criticizing someone's English is in poor taste though, since sage has such a broad international user base -- probably you just wanted the question to be more clear :)

( 2011-01-30 10:25:39 +0200 )edit

Hmm, I didn't think of it that way at all: noting that something doesn't make sense in English doesn't seem like a criticism to me, but an observation. Still, I'll defer to your judgment.

( 2011-01-30 10:34:05 +0200 )edit

thanks -- this is probably just one of those things that came across a little harsher than it was meant.

( 2011-01-30 12:59:16 +0200 )edit

You could also use the R interface in Sage to produce things, using their standard fitting and error bar tools.

Edit: R is probably the industry standard for statistical/modeling/fitting stuff. Of the statistical variety. http://www.r-project.org is the site. You can do

sage: r?


for basic info, and see the whole documentation at this page for more info (apparently it's not in the reference manual yet. I also have a sample at sagenb, but that site is so overloaded that it might be hard to access...

more

This should work in the notebook, incidentally, if that's where you 'live'.

( 2011-01-28 22:51:26 +0200 )edit

hmmm ok... I don't know R... how is it?

( 2011-01-29 14:10:49 +0200 )edit

Have a look at find_fit on the page http://www.sagemath.org/doc/reference...

more

thanks! hmm.... I found a way to fit without error bars (for now..) with scipy.

( 2011-01-29 14:10:27 +0200 )edit

hmmm I think that is just for very easy fits... And they use the function from scipy: scipy.optimize.leastsq ; that I use currentyl... But I don' t get xhi^2 and the correlation matrix and all that.... and that would be important

( 2011-01-29 14:47:52 +0200 )edit

Here is an old thread from sage-devel about error bars on plots, with some suggestions. One of them suggests using matplotlib, and perhaps this would also work for curve fitting.

more

error bar plotting is as I show in this example: http://rn0.ru/show/1429/ not the problem. But fitting is...I Havn't found a clue how to fit a functionto data in sage...

( 2011-01-28 14:27:40 +0200 )edit