next up previous index
Next: Using help in Octave Up: Interactive Fortran: Octave and Previous: Octave

Matlab

Now we repeat the same computation in Matlab. You will notice that nearly everything works almost exactly as in Octave. But there are small and subtle differences in some places, so watch out! 

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

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

    2.2000    3.3000    3.7000    3.8000    4.6000    4.9000    5.3000

  Columns 8 through 10 

    6.2000    6.4000    7.1000

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

    0.2000    0.4000    0.4000    0.2000    0.4000    0.2000    0.2000

  Columns 8 through 10 

    0.4000    0.2000    0.4000

>> one_by_sigma_2 = 1.0 ./ sigma_y .^ 2

one_by_sigma_2 =

  Columns 1 through 7 

   25.0000    6.2500    6.2500   25.0000    6.2500   25.0000   25.0000

  Columns 8 through 10 

    6.2500   25.0000    6.2500

>> S_sigma = sum(one_by_sigma_2)

S_sigma =

  156.2500

>> S_x = sum( x .* one_by_sigma_2 )

S_x =

  850.0000

>> S_y = sum( y .* one_by_sigma_2 )

S_y =

  720.6250

>> S_xx = sum( x .^ 2 .* one_by_sigma_2 )

S_xx =

   5.8375e+03

>> S_xy = sum( x .* y .* one_by_sigma_2 )

S_xy =

   4.5456e+03

>> one_by_delta = 1.0 / (S_sigma * S_xx - S_x ^ 2)

one_by_delta =

   5.2740e-06

>> a = one_by_delta * ( S_xx * S_y - S_x * S_xy )

a =

    1.8083

>> b = one_by_delta * ( S_sigma * S_xy - S_x * S_y )

b =

    0.5154

>> sigma_a = sqrt( one_by_delta * S_xx )

sigma_a =

    0.1755

>> sigma_b = sqrt( one_by_delta * S_sigma )

sigma_b =

    0.0287

>> chi_2 = sum(one_by_sigma_2 .* ( y - a - b .* x ) .^ 2)

chi_2 =

    3.8276

>>
Here is the first difference. The operators - and + are overloaded in Matlab  and there is no .- or .+. They work both on arrays and scalars exactly as you expect them to (and exactly as they do in Fortran).

The other difference is in the incomplete gamma function . There is no gammai in Matlab. Instead we must call gammainc, which swaps the order of its arguments:

>> Q = 1.0 - gammainc ( chi_2 / 2.0, ( 10 - 2 ) / 2 )

Q =

    0.8723

>>
In other words:  

\begin{displaymath}\texttt{gammainc}(x, a) = \texttt{gammai}(a, x) = \frac{1}{\Gamma(a)} \int_0^x e^{-t} t^{a - 1}\,\textrm{d} t
\end{displaymath}

In order to plot the fit we calculate y_prime as in Octave and as in our Fortran program:

>> y_prime = a + b * x

y_prime =

  Columns 1 through 7 

    2.3237    2.8391    3.3544    3.8698    4.3852    4.9006    5.4160

  Columns 8 through 10 

    5.9314    6.4468    6.9622

>>
Matlab's plot command  is not as nice as Gnuplot plot, and, for example, you can't use it to draw a plot with error bars. But drawing plots like that is so important that Matlab has a separate command for that. The command is errorbar , and you should use it in combination with plot, as shown below:
>> plot (x, y_prime)
>> hold on
>> errorbar(x, y, sigma_y, '.')
>> hold off
>>
The way it works is as follows. First we plot y' = a + bx versus xusing the plot command. Then we tell Matlab to hold on to the plot and add other elements to it rather than overwrite it with a new plot. Then we plot errorbars. The last argument in the errorbar command tells Matlab that the points given by (xi, yi) pairs should not be connected. Matlab draws them connected by default. The last command, hold off tells Matlab that the plot is finished, and no new elements will be added to it.

In order to generate a PostScript file  for inclusion in a TeX document, pull down the File menu on the plot window generated by Matlab and select Print. Then when the Print window pops up, tell it to send the output to a File instead of a Printer, and click on the Print button in the lower right corner. The generated PostScript will not be encapsulated though.


next up previous index
Next: Using help in Octave Up: Interactive Fortran: Octave and Previous: Octave
Zdzislaw Meglicki
2001-02-26