Alternative function: Normal cumulative density function (normcdf)

Download
normcdf2.m

Details
This function implements the basic functionality of the normcdf function in Matlab's stats toolbox. Similar to normcdf it takes three inputs, the x axis (x), mean of the distribution (mu) and the standard deviation (sigma) and generates a normal cumulative density function using the equation:

 y = 0.5*(1+erf((x-mu)/sqrt(2*sigma^2)))

See also normpdf2 for the equivalent probability density distribution.

Example

% Define parameters
mu = -1; % Mean
sigma = 0.5; % Standard deviation
x = -10:0.1:10; % x axis

% Generate curve
y1 = normcdf2(x, mu, sigma);
y2 = normcdf2(x, mu-1, sigma+2);
y3 = normcdf2(x, mu+6, sigma+0.5);

% Plot 
figure
plot(x, [y1', y2', y3'])
ylabel('Prob.')
xlabel('x')
legend({...
    ['mu=', num2str(mu), ', sigma=', num2str(sigma)], ...
    ['mu=', num2str(mu-1), ', sigma=', num2str(sigma+2)], ...
    ['mu=', num2str(mu+6), ', sigma=', num2str(sigma+1)], ...
    })

Alternative function: Normal probability density function (normpdf)

Download
normpdf2.m

Details
This function implements the basic functionality of the normpdf function in Matlab's stats toolbox. Similar to normpdf it takes three inputs, the x axis (x), mean of the distribution (mu) and the standard deviation (sigma) and generates a normal probability density function:

y = 1/sqrt(2*pi*sigma^2) * exp(-(x-mu).^2 / (2*sigma^2))

See also normcdf2 for the equivalent cumulative distribution.

Example

% Define parameters
mu = -1; % Mean
sigma = 0.5; % Standard deviation
x = -10:0.1:10; % x axis

% Generate curve
y1 = normpdf2(x, mu, sigma);
y2 = normpdf2(x, mu-1, sigma+2);
y3 = normpdf2(x, mu+6, sigma+0.5);

% Plot 
figure
plot(x, [y1', y2', y3'])
ylabel('Prob.')
xlabel('x')
legend({...
    ['mu=', num2str(mu), ', sigma=', num2str(sigma)], ...
    ['mu=', num2str(mu-1), ', sigma=', num2str(sigma+2)], ...
    ['mu=', num2str(mu+6), ', sigma=', num2str(sigma+1)], ...
    })

Converting hexadecimal to binary

See also: Base conversion of positional number systems for other conversion functions.

Download
hex2bin2.m

Hexadecimal and binary
Hexidecimal and binary, like "normal" decimal are positional number systems, where individual digits have a value determined by their position in the sequence (see this post for a more detailed explanation). The difference between the three is that decimal has a base of 10 - it uses 10 numbers, 0-9, where hexadecimal has a base of 16 (using 0-9 and A, B, C, D, E and F) and binary has a base of 2 (using 0 and 1).

There are two ways to convert between hexadecimal and binary, mathematically, or by simply using a lookup table. Those interested in the mathematical process should check out the hex2dec2 and dec2bin2 functions (bin = dec2bin2(hex2dec2(hex))), which describe how to do each step. But if you'd just like to get it done and don't care about the maths behind it, just use hex2bin2. hex2bin2 uses the following lookup table to get the binary value of each hex digit.

Dec = Hex = Bin
 0   =   0   =   0000
 1   =   1   =   0001
 2   =   2   =   0010
 3   =   3   =   0011
 4   =   4   =   0100
 5   =   5   =   0101
 6   =   6   =   0110
 7   =   7   =   0111
 8   =   8   =   1000
 9   =   9   =   1001
10   =   A =   1010
11   =   B =   1011
12   =   C =   1100
13   =   D =   1101
14   =   E =   1110
15   =   F =   1111

Converting hexadecimal to decimal

See also: Base conversion of positional number systems for other conversion functions.

Download
hex2dec2.m

Hexadecimal numbers
"Normal" decimal numbers use a base 10 system, most likely because humans have 10 fingers. Hexidecimal is a base 16 number system (presumably because computer programmers have 16 fingers), so whereas the decimal number system has 10 numbers (0-9), hexadecimal has 16 (0-9, and A, B, C, D, E and F).

Like decimal (and binary) hexidecimal is a positional number system, where a numbers position determines its value. This means that calculating the value of a hexidecimal number is exactly the same process as calculating the value of a decimal number. Observe:

Consider the decimal value 2345. The value is obvious. It's 2345. But it can be broken down into the sum of 2000 + 300 + 40 + 5. And each of these values depends on the number of zeros, which is determined by position. Positions start from 0 and count up right to left, so the 2 in 2345 is in position 3, the 3 is in position 2, the 4 is in position 1 and the 5 in position 0. Mathematically this means the value of the 2 is 2*10^3 = 2000, of the 3 is 3*10^2 = 300, and so on. Overall the total value of the number is:

(2*10^3) + (3*10^2) + (4*10^1) + (5*10^0) = 2345.

For hexadecimal numbers, exactly the same process applies, except using a base of 16. The only caveat being that the letters, if there are any, need to be replaced with their equivalent decimal numbers first. So what's the value of the hex number 2345 (notice that for hex values without letters, it's not immediately obvious it's a hex number - so it should be written as either 0x2345 or 234516 to avoid this ambiguity).

0x2345 (hex) = (2*16^3) + (3*16^2) + (4*16^1) + (5*16^0) = 9026 (decimal).

If there are letters in the hex number, for example: 0xF34B, these need to be converted to decimal values first, which is simple: A=10, B=11, C=12, D=13, E=14, F=15. So,

0xF34B = (F*16^3) + (3*16^2) + (4*16^1) + (B*16^0) , or
0xF34B = (15*16^3) + (3*16^2) + (4*16^1) + (11*16^0) = 62283

This process is simple to implement in Matlab, and is very similar to the process to convert decimal values to binary (see dec2bin2). The only complication is the handling of the letters; in Matlab hexadecimal numbers must be strings, because they contain letters. This means that isn't not possible to directly perform mathematical operations on them, such as the above conversion, but Matlab has a number of convenient features for dealing with strings that we'll exploit.

Overall, this function will convert the hexadecimal input to decimal numbers, work out the value of each based on its position, the calculate the sum of all these values to give the final decimal value. It provides an alternative to Matlab's hex2dec function.

Event / spike triggered average

Download
evTrigAvg.m

What caused an event?
Generally, event triggered averages are used to try and determine what caused an event to happen. They look at what preceded a number of events and take the average. In neuroscience events might be a spike recorded from a neuron in response to a certain stimulus. If the stimulus was white noise, for example, taking the time window from before the spike and averaging it would give the average stimulus that caused the neuron to fire. In finance the event could be a price change, and one might like to know if a consistent pattern occurs before similar price changes in order to react to them in the future.

evTrigAvg is designed to be general and agnostic to exactly what events and stimuli its dealing with are. All it needs is the "stimulus" and a logical vector of events (0s indicating no events, 1s indicating an event has occurred), and a window size to average over (in points).

Example
Let's generate some fake neural data as an example. Imagine an experiment where white noise is delivered to a subject and activity is recorded from a single neuron. The raw activity from the neuron will be noisy and will normally require preprocessing first, which will include filtering and event (spike) detection. Event detection is often done using a threshold of some value *the root mean square of the trace, where anything above this is marked as a spike. Or, if multiple units have been recorded, by much more complicated "spike sorting". I mention this mainly because it presents a fascinatingly complex problem, but it's not something we need to worry about here. The goal of preprocessing is to obtain a clean vector of events represented as 0s or 1s, we'll imagine the preprocessing has already been done, and we have this clean vector of events (spikes). We want to know what property of the stimulus caused each event.

nPoints = 10000;
stim = randn(1,nPoints);
resp = stim>2;
resp = [zeros(1,100), resp(1:end-100)];

nPoints determines the number of points in the stimulus (stim) and response (resp). stim is random noise drawn from a normal distribution. resp will be our vector of events (spikes), and for the sake of demonstration spikes are generated whenever the stim is greater than the entirely arbitrary value 2. Obviously, this means we know what stimulus property that is eliciting the spikes already, but pretend this is real data that we don't yet know anything about it. The line stim>2 returns a logical vector of 0s and 1s, which we'll also offset slightly relative to stim to emulate the latency that would exist in a real system (spikes aren't instant) by dumping 100 zeros at the start.

Let's plot out the first 1000 points of data we've just magically generated without having to do any pesky experiments.

subplot(2,1,1), plot(stim(1:1000))
line([1, 1000], [2, 2], 'LineStyle', '--', ...
    'color', 'k')
ylabel('Mag, arb'), title('Stimulus')
subplot(2,1,2), plot(resp(1:1000), 'r')
ylabel('Spikes'), title('Response')
xlabel('Time, pts')

Alternative function: Net future value of a cash flow

See also the time value of money for other related Matlab functions.

Download
nfvFlow.m

Introduction
The future value (FV) of a single lump sum of cash is its present value (PV, its value right now) plus interest earned over a period of time. The future value of a cash flow is the sum of the FVs of the money received (or spent) at each time period. As each time periods go by, the amount of periods over which a subsequent income can potentially earn interest reduces.

For an introduction to future values and simple and compound interest, and a simple Matlab function for calculating the future value of a lump sum, see this post.

nfvFlow.m is an alternative to Matlab's fvvar function in the financial toolbox. nfvFlow includes a flag that allows interest payments to be made in advance or in arrears (intPaid=1 for advance, or =0 for arrears). Payments (payments) and periods (periods) can be explicitly specified, or entered as a single value if regular (see below for examples). The interest rate (rate) should be a decimal value, not a percentage. The outputs give the overall net future value of the cash flow (nfv) and the future value of each payment made at each time point (fvAtPeriod).

[nfv, fvAtPeriod]= nfvFlow(payments, rate, periods, intPaid)

For example, for a regular payment of £2000, over 4 periods with a rate of 5% and interest paid at the end of time periods (see first example below for details):

[nfv, fvAtPeriod] = nfvFlow(2000, 0.05, 4, 0).

Fitting a better psychometric curve (Wichmann and Hill 2001)

Code
GitHub: Matlab code and Python package

Introduction
Psychometric curves quantify the relationship between stimulus parameters and an observer's subjective experience. For an introduction to the concept of psychometric curves in signal analysis theory and a function and instructions to fit a basic curve in Matlab, see Introduction to psychometric curves and fitting a simple psychometric function. This article will focus on fitting a more complex curve with additional parameters that allow for the subjects fallibility, which improves estimation of the key parameter, discrimination sensitivity.

The psychometric function used here is from the paper The psychometric function:I. Fitting, sampling, and goodness of fit, F. Wichmann and N. Hill 2001. This paper, and its companion paper, provide detailed theory and methods for fitting a psychometric function, assessing its goodness of fit, and finding confidence intervals for the ranges of the parameters estimated. This article will demonstrate how to fit Wichmann and Hill's psychometric function in Matlab using the same data from a fictional 2AFC experiment in this post. This will allow an interested reader to make a direct comparison between the Wichmann and Hill curve and the approach using Matlab's glmfit function with a built in binomial distribution and a logit link function.