QuickTip: Adding Search Paths in Octave

Suppose you want a utility (such as matrixTeX.m) to be visible to every octave session.

1. Select a suitable folder to keep your local octave files in (e.g. /home/sachin/octave/ on my Linux box). If such a folder doesn't exist, make one.

2.  Fire up an octave session

3. In Octave, add the path to the path variable:

octave:> addpath(genpath("/home/sachin/octave")) 

4. List directories in default path (your output will likely look different)

octave:> path 

/home/sachin/octave
/home/sachin/octave/io-2.4.11
/home/sachin/octave/io-2.4.11/doc
/home/sachin/octave/io-2.4.11/packinfo
/home/sachin/octave/io-2.4.11/templates

etc.

4. So far, we haven't permanently written the modified path to disk. To do this:

octave:> savepath

Nicer Octave Plots for Presentations/Documents

GNU Octave has a very capable plotting system (which is transitioning away from gnuplot as its backend, even as we "speak"). Quite often, you generate a figure and make it look like you want:

You then want to incorporate it into some other document - say a presentation or a document.

So you say:

print -dpng test0.png

Note that you could also use -depsc2 or -dpdf to generate EPS or PDF formats among a whole host of other possibilities.

You get a figure that looks like:

Clearly imperfect.

There are a number of ways to improve the image.

These include writing to tikz, TeX, which gives you exquisite control over all features, but this can take too much effort.

Or you could spend time fiddling with adjustable parameters to make the plot look just right. If you are trying to work on "the" central plot of a paper, or a PhD thesis, then such care is certainly warranted.

But if you are merely dissatisfied (as I am) with the size and font of the text, there is an easy fix.

print -dpng -FHelvetica:18 test1.png

Conditional Column Switching in Octave or Matlab

Consider a simple Matlab/Octave exercise.

I have an N by 2 vector x. Thus, x(i,1) and x(i,2) are the entries of the i-th row. Say I want to rearrange x so that the first entry in any row is larger than the second. That is, I want to switch the entries of a row, only if a certain condition is met.

It is easy to write this in a element-by-element fashion

for i = 1:N
  if x(i,2) > x(i,1)
    c = x(i,1);
    x(i,1) = x(i,2);
    x(i,2) = c;
  endif
endfor

But we all know that element-by-element operations in Matlab or Octave suck! Here is another vectorized method:


SwitchCond = x(:,2) < x(:,1);
x(SwitchCond,[1,2]) = x(SwitchCond,[2,1]);

Besides being faster, it is also more succinct.

Linear Least Squares in GNU Octave: Part Deux

I recently blogged about doing LLS in GNU Octave, getting not only the best fit parameters, but also the standard error associated with those estimates (assuming the error in the data is normally distributed).

Octave as an built-in function to do ordinary least squares: ols. While it does not directly report the standard errors associated with the regressed parameters, it is relatively straightforward to get them from the variables returned.

Here is how you would solve the same example as before:

Nobs = 10;
x = linspace(0,1,Nobs)';
y = 1 + 2 * x + 0.05 * randn(size(x));
X = [ones(Nobs,1) x];

[beta sigma] = ols(y,X)
yest         = X * beta;

p    = length(beta);             % #parameters
varb = sigma * ((X'*X)\eye(p));  % variance of "b"
se   = sqrt(diag(varb))          % standard errors

Linear Least Squares in GNU Octave with Standard Errors

Consider a linear model with \(N\) observations of the quantity \(Y\), as a function of \( p\) regressors, \(Y = \sum \beta_i X_i\).

\[\begin{bmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_N \end{bmatrix} = \begin{bmatrix} X_{11} & X_{12}  &  ... & X_{1p} \\ X_{21} & X_{22}  & ... & X_{2p} \\  & & \ddots & \\ X_{N1} & X_{N2}  & ... & X_{Np} \end{bmatrix} \begin{bmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_p \end{bmatrix} + \epsilon,\]

where \(\epsilon\) is normally distributed error. Each row corresponds to an observation, and each column and associated \(\beta\) corresponds to a parameter to be regressed

In general, this can be written as \(Y = X \beta\).

As a simple illustrative case consider fitting the model \(Y = \beta_0 + \beta_1 X\). The above equation becomes: \[\begin{bmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_N \end{bmatrix} = \begin{bmatrix} 1 & X_1 \\ 1 & X_{2}  \\ \vdots & \vdots \\   1 & X_{N} \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix}.\]

Given the expressions in the wikipedia article on LLS, we can easily write an Octave program that takes in y and X as above, and spits out the best estimate, the standard error on those estimates, and the set of residuals.

%

% For y = X b + gaussian noise:

% compute b, standard error on b, and residuals

%

function [b se r] = lls(y, X)

        

    [Nobs, p] = size(X);          % size of data

    b   = (X'*X)\X'*y;            % b = estimate beta

    df  = Nobs - p;               % degrees of freedom

    r   = y - X * b;              % residuals

    s2  = (r' * r)/df;            % SSER

    varb = s2 * ((X'*X)\eye(p));  % variance of "b"

    se   = sqrt(diag(varb));      % standard errors

endfunction

To test the model I "create" the data y = 1 + 2x + white noise.

> x = linspace(0,1,10)';

> y = 1 + 2 * x + 0.05 * randn(size(x));

> X = [ones(t,1) x];

> [beta se] = lls(y,X)

beta =

   0.98210

   2.03113

se =

   0.039329

   0.066302

A plot of the data and the regression looks like this:

Non-Uniform Quadrature Points in Octave/Matlab

Suppose you are given a bunch of points in the form of a vector x (between "a" and "b") and an accompanying density function p(x)>0. For concreteness, consider

x = [0:0.1:1.0]'; px = 0.3 - (x-0.5).^2;

The top subfigure depicts the density function p(x), and the cumulative density function C(x). The bottom subfigure depicts the  quadrature points "z" (circles), and the lines mark the end of the intervals, using N=11.

This density function is shown in the picture above.

We want to create a N*1 vector "z" of quadrature points that are distributed according to the density p(x). We also want a vector  "dz" corresponding to the width of strip belonging to a particular quadrature point. As an additional constraint, we assume that we want z(1) = a, and z(N)  = b to match the end-points.

The program attached at the end is able to do this relatively efficiently. In fact, I used the following command to generate the picture above:

[z dz] = GridDensity(x,px,11,1);

GNU Octave Code:

%

%  PROGRAM:

% Takes in a PDF or density function, and spits out a bunch of points in

%       accordance with the PDF

%

%  INPUT:

%       x  = vector of points. It need *not* be equispaced

%       px = vector of same size as x: probability distribution or

%            density function. It need not be normalized but has to be positive.

%   N  = Number of points >= 3. The end points of "x" are included

%       necessarily

%       Pt = Optional argument. If present, then some plotting.

% 

%  OUTPUT:

%       z  = Points distributed according to the density

%       hz = width of the "intervals" - useful to apportion domain to points

%            if you are doing quadrature with the results, for example.

%

%  (*) Sachin Shanbhag, March 5, 2012

%  GNU Public License

%

function [z h] = GridDensity(x,px,N,Pt)

  npts = 100;                              % can potentially change

  xi   = linspace(min(x),max(x),npts)';    % reinterpolate on equi-spaced axis

  pint = interp1(x,px,xi,'spline');        % smoothen using splines

  ci   = cumtrapz(xi,pint);                

  pint = pint/ci(npts);

  ci   = ci/ci(npts);                      % normalize ci

  alfa = 1/(N-1);                          % alfa/2 + (N-1)*alfa + alfa/2

  zij  = zeros(N,1);                       % quadrature interval end marker

  z    = zeros(N,1);                       % quadrature point

  z(1)  = min(x);  

  z(N)  = max(x); 

%

% ci(Z_j,j+1) = (j - 0.5) * alfa

%

  beta  = [0.5:1:N-1.5]'*alfa;

  zij   = [z(1); interp1(ci,xi,beta,'spline'); z(N)];

  h        = diff(zij);

  clear beta;

%

% Quadrature points are not the centroids, but rather the center of masses

% of the quadrature intervals

%

  beta     = [1:1:N-2]'*alfa;

  z(2:N-1) = interp1(ci,xi,beta,'spline');

%

% Some plotting if required

%

  if(nargin>3)

    subplot(2,1,1)

    plot(xi,ci,'b-','LineWidth',2,xi,pint,'r-','LineWidth',2);

    axis([min(xi),max(xi)]);

    legend('CDF','PDF')

    title('Visualization')

    xlabel('x');

    ylabel('CDF(x)/PDF(x)')

    subplot(2,1,2)

    plot(z,0.5*ones(size(z)),'ro');

    axis([min(xi),max(xi)]);

    xlabel('z');

    hold on;

    for i = 2:N

      X = [zij(i); zij(i)];

      Y = [0; 1];

      plot(X,Y,'b')

    endfor

    hold off;

  endif

endfunction

Subfunctions in Octave/Matlab

Consider the contents of the following Matlab program "mainFunction.m".

%

% A function that integrates the function

% given by the routine subFunction(x)

% over the range a and b

%

function I = mainFunction(a,b)
   a = -5;

   b = 5;

   I = quad('subFunction', a, b)

%  keyboard  % uncomment if needed

endfunction
 

%

% some nontrivial function

%

function y=subFunction(x)

   if( x > 0)

      y = x^2;

   else

      y = x;

   endif

endfunction

The function to be integrated, subFunction(x), can be saved as a separate file called subFunction.m, or, as in this case, it can be included in the same file as the main function. There are advantages to both these techniques.

If I store functions in their own individual files, I can call them from any other function. However, I can very easily end up with a lot of files. This may not necessarily be a bad thing, but it can be a challenge beyond a point.

If I use the subfunctions (as above), then the subfunctions themselves remain invisible outside the file in which they are contained. They may be accessed by the main function and other subfunctions within the same file.

Unfortunately, you cannot use "subfunctions" with script files. A useful workaround is to wrap a function declaration around your script, and use the command keyboard (as shown in the commented line in the example above) to transfer control out.

From here you can examine all the variables local to the "script".

When you are done, say dbquit, or return.