Octave

In this section I show you the Octave session. On our systems Octave lives on

/afs/ovpit.indiana.edu/@sys/gnu/bin/octave

and, if you have the ../gnu/bin directory in front of your command search path, all you need to do to run Octave is to type

$ octave

Octave may take a while to load, because it is a very large application. But once you get it going you should be greeted by a message and a prompt, which looks like that:

gustav@blanc:../gustav 12:22:45 !530 $ octave
Octave, version 2.0.13 (sparc-sun-solaris2.6).
Copyright (C) 1996, 1997, 1998 John W. Eaton.
This is free software with ABSOLUTELY NO WARRANTY.
For details, type `warranty'.

octave:1>

We begin our computation by defining and instantiating Octave arrays  x_i, y_i and $\sigma _i$:

octave:1> x = [ 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 ]
x =

   1   2   3   4   5   6   7   8   9  10

octave:2> y = [ 2.2 3.3 3.7 3.8 4.6 4.9 5.3 6.2 6.4 7.1 ]
y =

 Columns 1 through 9:

  2.2000  3.3000  3.7000  3.8000  4.6000  4.9000  5.3000  6.2000  6.4000

 Column 10:

  7.1000

octave:3> sigma_y = [ 0.2 0.4 0.4 0.2 0.4 0.2 0.2 0.4 0.2 0.4 ]
sigma_y =

 Columns 1 through 8:

  0.20000  0.40000  0.40000  0.20000  0.40000  0.20000  0.20000  0.40000

 Columns 9 and 10:

  0.20000  0.40000

octave:4>

Now we calculate $1/\sigma_i^2$: 

octave:5> one_by_sigma_2 = 1.0 ./ sigma_y .^ 2
one_by_sigma_2 =

 Columns 1 through 8:

  25.0000   6.2500   6.2500  25.0000   6.2500  25.0000  25.0000   6.2500

 Columns 9 and 10:

  25.0000   6.2500

octave:6>

Observe that every time you perform an operation, Octave returns the result immediately. You can ask it not to print it, by terminating your commands with semicolons. It works the same way as in Mathematica: the semicolon is the output silencer. 

The notation ./ and .^ means element by element division and exponentiation . In Fortran all normal array operations are element by element, so no special notation is needed. In Matlab and Octave certain operations, such as *, may have a meaning of a scalar product, for example, or of other contraction operations, such as multiplication of matrices or multiplication of a matrix by a vector.

Having evaluated $1/\sigma_i^2$ we can now finish our computation very quickly, proceeding roughly in the same way we did in the Fortran program:

octave:6> S_sigma = sum(one_by_sigma_2)
S_sigma = 156.25
octave:7> S_x = sum( x .* one_by_sigma_2 )
S_x = 850.00
octave:8> S_y = sum( y .* one_by_sigma_2 )
S_y = 720.62
octave:9> S_xx = sum( x .^ 2 .* one_by_sigma_2 )
S_xx = 5837.5
octave:10> S_xy = sum( x .* y .* one_by_sigma_2 )
S_xy = 4545.6
octave:11> one_by_delta = 1.0 / (S_sigma * S_xx - S_x ^ 2)
one_by_delta =  5.2740e-06
octave:12> a = one_by_delta * ( S_xx * S_y - S_x * S_xy )
a = 1.8083
octave:13> b = one_by_delta * ( S_sigma * S_xy - S_x * S_y )
b = 0.51539
octave:14> sigma_a = sqrt( one_by_delta * S_xx )
sigma_a = 0.17546
octave:15> sigma_b = sqrt( one_by_delta * S_sigma )
sigma_b = 0.028706
octave:16> chi_2 = sum(one_by_sigma_2 .* ( y .- a .- b .* x ) .^ 2)
chi_2 = 3.8276
octave:17> Q = 1.0 - gammai ( ( 10 - 2 ) / 2, chi_2 / 2.0 )
Q = 0.87233
octave:18>

The incomplete gamma function  is predefined in Octave. You will find the operational definition of gammai on

/afs/ovpit.indiana.edu/@sys/gnu/share/octave/2.0.13/m/specfun/gammai.m

and when you look inside that file, you'll see that all it does is to call gammainc with a changed order of arguments. In turn, gammainc is a built-in function that has been compiled into Octave binary.

To view the results of our $\chi ^2$ fit, we begin by evaluating y_prime as we did in our Fortran program, i.e., y'_i = a + b x_i:

octave:17> y_prime = a + b * x
y_prime =

 Columns 1 through 9:

  2.3237  2.8391  3.3544  3.8698  4.3852  4.9006  5.4160  5.9314  6.4468

 Column 10:

  6.9622

octave:18>

Now we have to construct a matrix, which looks very much the same as our plot data file written by program chi:

octave:18> data = [ x' y' sigma_y' y_prime' ]
data =

   1.00000   2.20000   0.20000   2.32367
   2.00000   3.30000   0.40000   2.83906
   3.00000   3.70000   0.40000   3.35445
   4.00000   3.80000   0.20000   3.86984
   5.00000   4.60000   0.40000   4.38523
   6.00000   4.90000   0.20000   4.90062
   7.00000   5.30000   0.20000   5.41601
   8.00000   6.20000   0.40000   5.93140
   9.00000   6.40000   0.20000   6.44679
  10.00000   7.10000   0.40000   6.96218

octave:19>

The primes adoring x, y, sigma_y and y_prime transpose  row vectors into column vectors. Otherwise instead of obtaining a matrix we would get one row vector 40 elements long.

The last step is to issue a plot command, which looks almost exactly the same as the Gnuplot  command we've used in our Gnuplot script:

octave:19> gplot [0.5:10.5] [2 : 8] data using 1:2:3 with errorbars, \
                                    data using 1:4 with lines

This should not be surprising, because Octave uses Gnuplot for all its graphics.

You can use other Gnuplot commands, slightly renamed in Octave (gset instead of set for example), in order to generate an encapsulated PostScript plot for inclusion in a LaTeX    document:

octave:23> gset term postscript eps 22
octave:24> gset output "octave-chi.eps"
octave:25> gset title "chi^2 fitting"
octave:26> replot 
octave:27>

Compare these to Gnuplot commands shown in Figure 2.7 on page .