http://www.agnld.uni-potsdam.de/~marwan/matlab-tutorials/html/statistics.html#37

Statistics with Matlab

Contents

Matlab provides commands and operations for the computation of basic statistics.

randn('seed', 1);
x = randn(1000,1);
y = randn(1000,1);

The mean of a series can be computed by

mean(x)

ans =

    0.0020

If there are some NaNs in the vector x, this command will also be NaN. For this case we can use the command nanmean.

Further commands are

median

median(x)

ans =

    0.0197

standard deviation

std(x)

ans =

    1.0096

variance

var(x)

ans =

    1.0193

largest value

max(x)

ans =

    3.0206

smalles value

min(x)

ans =

   -3.4333

value range

range(x)

ans =

    6.4539

interquartile range (robust estimate for the spread of the data)

iqr(x)

ans =

    1.3076

covariance matrix

cov(x,y)

ans =

    1.0193   -0.0011
   -0.0011    0.9783

correlation matrix

corrcoef(x,y)

ans =

    1.0000   -0.0011
   -0.0011    1.0000

Sometimes it is desired to know the position of the kargest or smalles value in the vector. The commands max and min can also provide this information:

[a i] = max(x)

a =

    3.0206


i =

   317

Now the position of the largest value a in the vector x is in variable i.

To sort the values in the vector x we use the sort command:

xs = sort(x);

or with the index vector of the sorted values

[xs i] = sort(x);

A histogram of the data can be computed with

hist(x)

The default settings for a histogram are ten intervalls. To change this to another number of intervalls, e. g. 50, we call

hist(x, 50)

The number of counts h per intervall i can be accessed by

[h i] = hist(x,50);

Now we would like to know the 5% and 95% percentiles of the distribution of the series x.

prctile(x,[5 95])

ans =

   -1.6467    1.6592

The minimum, maximum, lower and upper quartile and median values can be represented by a box and whisker plot

boxplot([x y])

Where mean and standard deviation are the first two moments of the distribution, skewness and kurtosis are related to the third and fourth moments of the distributions. The skewness reflects the shape or asymmetry of the distribution: if it is negative, the data spread out more to the left of the mean, if it is positive, the data spread out more to the right.

skewness(x)

ans =

   -0.0786

The kurtosis represents how outlier-prone the distribution is. The kurtosis of the normal distribution is 3. Distributions that are more outlier-prone than the normal distribution have kurtosis greater than 3; distributions that are less outlier-prone have kurtosis less than 3.

kurtosis(x)

ans =

    3.0258

All moments of a distribution can be derived with the command moment.

To assess statistical significance, the bootstrap statistics is sometimes helpful. With bootstrapping we can get new realisations of the original data series by a random resampling. The Matlab command bootstrp can be applied for a bootstrap statistic. To compute a bootstrap statistic of the mean of our vector x by using 500 new realisations, we call

m = bootstrp(500, 'mean', x);

To show the variation of the mean across all bootstrap realisations, we plot the histogram

hist(m)

As expected, the mean spreads around the origin.

Linear regression models can be useful for the study of relations between two data series. Matlab provides different commands to estimate linear regression coefficients and corresponding statistics.

Least squares fit can be performed by the command regress. Assume a linear system

randn('seed',1);
x = [1:10:500]';
y = 2.5 * x + 80 + 50 * randn(50,1);

plot(x,y,'.')

and estimate the linear fit and the confidence intervals within 95% with

[b, bint] = regress(y, [ones(length(x),1), x], 0.05)

b =

   78.1034
    2.4832


bint =

   50.3461  105.8608
    2.3859    2.5805

In order to get not only an estimate for the regression coefficient, but also the offset, we added a vector consisting of ones in the second argument of regress. The result b yiels the offset 85 and the regression coefficient 2.5, the corresponding confidence intervals are stored in bint.

hold on
plot(1:500, b(2) * [1:500] + b(1))

Residuals can be automatically computed with

[b, bint, residuals] = regress(y, [ones(length(x),1), x], 0.05);
clf, plot(residuals), ylabel('Residuals')

For nonlinear systems, the command polyfit is more appropriate. We consider the following nonlinear system

randn('seed', 1);
y = sin(x/50) ./ x + .002 * randn(50,1);
clf
plot(x,y,'.')

and fit a polynom of order 5

p = polyfit(x, y, 5);

hold on
plot(x,polyval(p,x))

Warning: Polynomial is badly conditioned. Add points with distinct X
         values, reduce the degree of the polynomial, or try centering
         and scaling as described in HELP POLYFIT.

Error bounds can be derived by

[p s] = polyfit(x, y, 5);
[y_fit delta] = polyval(p, x, s);

Warning: Polynomial is badly conditioned. Add points with distinct X
         values, reduce the degree of the polynomial, or try centering
         and scaling as described in HELP POLYFIT.

The variable delta contains the estimate of the standard deviation of the error in predicting a future observation at x by P(x). As error bounds we use 2 delta corresponding roughly to the 95% confidence interval.

plot(x, y_fit + 2 * delta, ':', x, y_fit - 2 * delta, ':')

Exercise

1. Consider the file ecg.dat and compute and interprete the mean, standard deviation, skewness and kurtosis. Plot and interprete the histogram (chose an appropriate number of intervals).

2. The file eeg.dat contains EEG measurements on 59 channels (arranged as columns). Compute a matrix of correlation coefficients between all channels.

3. Consider the first 120 data points in the file ecg.dat. Fit polynoms of different orders to this sequence of data. Compute bootstrap estimates of the last coefficient of the polynom and present the result in a histogram.

Published with MATLAB® 7.3

This toolbox is a very simple polynomial manipulation package.
I originally wrote it some years ago as a means of learning how to work with objects in Matlab. When I decided to publish this package, I completely rewrote it from scratch just this week, as the old one had some bugs in it.

Why publish this at all? Because I think writing such a tool is a marvelous way to learn about OOP and the use of classes in Matlab. A second reason is that this toolbox is actually of use to me on occasion, whenever I need to manipulate simple polynomials in one or several variables. If, as is the case with me and you don't own the symbolic toolbox, you may find it of interest too. I have included the functions orthpoly and gaussquadrule, as neat applications of the sympoly tools.

I've also put in a helper document that discusses things I felt important in writing a toolbox like this.

If you wish to use these tools, they are quite easy to use. A few quick examples:

% This creates 3 sympoly objects in your
% workspace: x, y, z.
sympoly x y z

% Add 1 to x, put the result in p.
p = x+1

% arbitrary expressions
q = (x-1)^3 + x*y*z - (x+1)*(z-1)

There are many other examples in the ReadMe file, some involving arrays of sympolys, as well as many more.

Whatever use you do find for this toolbox, have fun with it. I did. I even get some use from it occasionally.

If by some amazing chance you do find any bugs, please e-mail me.

code.zip

Polyfitn is an extension of polyfit, allowing the user to create models with more than one independent variable. It also allows the user to specify a general model, for example, a quadratic model, with constant and quadratic terms, but no linear term.
For example, to fit a polynomial model to points selected from a cosine curve, we will only need the even ordered terms.
x = -2:.1:2;
y = cos(x);
p = polyfitn(x,y,'constant x^2 x^4 x^6');

p.Coefficients
ans =
    [0.99996 -0.49968 0.041242 -0.0012079]
The coefficients won't be exact of course, as I used only a finite number of terms for what is essentially a truncated Taylor series, and I had only a finite amount of points to build the model from. The exact first 4 coefficients for a cosine series would have been:
>> [1 -1/2 1/24 -1/720]
ans =
            1 -0.5 0.041667 -0.0013889

so we got the expected result.

Of course, polyfitn works in higher dimensions, as it was this problem it was really designed to solve.

x = rand(100,1);
y = rand(100,1);
z = exp(x+y) + randn(100,1)/100;
p = polyfitn([x,y],z,3);

The result can be converted into a symbolic form to view the model more simply. Here I'll use my sympoly toolbox, but there is also a polyn2sym function provided.

polyn2sympoly(p)
ans =
    0.43896*X1^3 + 1.4919*X1^2*X2 + 0.041084*X1^2 + 1.4615*X1*X2^2 - 0.095977*X1*X2 + 1.2799*X1 + 0.56912*X2^3 - 0.15306*X2^2 + 1.361*X2 + 0.94819

And of course, parameter error estimates are generated for those who wish to determine the significance of the terms generated.

I've also provided tools for evaluating these models, and to differentiate a model.

A caveat - beware the use of high order polynomials to fit your data. Just because a low order model works, a high order model is not necessarily better. High order polynomials often suffer from severe ringing between the data points. Always plot your data. Think about the model you will be building. Then plot the resulting model. Use your eyes to validate the result, more so than simply looking at an r-squared coefficient (although I do return that parameter too.)

If you do find that a high order polynomial mode is necessary because your curve is simply too complicated, consider using a regression or smoothing spline model instead.

Hello,

I have data for two independent variables and one dependent variable (obtained from experiment). I need to fit this data using linear regression to a 10 coefficient third degree polynomial equation - for the engineers among you, this is the standard equation for specifying refrigeration compressor performance. This can be done quite easily in EES by entering data into a table and selecting the option "Linear Regression" from the Tables menu. How do I achieve the same in Matlab? My equation is of the form X = C1 + C2.(S) + C3.(D) + C4·(S2) + C5 · (S·D) + C6 · (D2 ) + C7 · (S3 ) + C8 · (D·S2 ) +C9 · (S·D2 ) + C10 · (D3 ) Here, X is the dependent variable and S and D are the independent variables. The numbers next to S and D indicate the power to which they are raised. C1 to C10 are the coefficients that need to be calculated.

Using the 'regress' function gives me a constant of 0, and warns me that my design matrix is rank deficient. Using the curve fitting toolbox (cftool - polynomial option) gives me ridiculous values for the coefficients (p00 = -6.436e15). When I try to input a custom equation in the cftool, it is switching to non-linear regression and asks me to input guess values for the coefficients, which I don't want to do. What other functions are available that I might use to perform this regression and how do I implement them. Any suggestions/ help/ recommendations would be greatly appreciated.

Thanks very much.

Gautam.

---------------------------------------------------------------------------------------------------

Please post your regression equation and 15-10 rows of your X, S and D data.

It’s difficult to determine what the problem may be without that information.

---------------------------------------------------------------------------------------------------

>> X = [-10 101; -10 114; -10 129; -10 144; -5 96; -5 106; -5 119; -5 133; 5 92; 5 97; 5 108; 5 122];

>> Y = [72.03; 67.90; 62.77; 57.55; 88.53; 84.70; 79.91; 73.73; 123.43; 120.91; 115.44; 107.83];

>> x1 = X(:,1);

>> x2 = X(:,2);

>> cftool

I used a custom equation, Y = f(x1, x2) and the custom equation is:

a + b*x1 + c*x2 + d*(x1^2) + e*(x1*x2) + f*(x2^2) + g*(x1^3) + h*(x1^2*x2) + i*(x1*x2^2) + j*(x2^3)

What turned up in the Results box was this:

Complex value computed by model function, fitting cannot continue. Try using or tightening upper and lower bounds on coefficients.

Any help?

---------------------------------------------------------------------------------------------------

You have only three distinct values in the first column of your X matrix, so you won't be able to estimate constant, linear, quadratic, and cubic effects for that column. Here's how you can use all terms in the cubic model except x1^3:

fitlm(X,Y,[0 0;1 0;0 1;1 1;2 0;0 2;0 3;1 2;2 1])

Alternatively, if you know that a third order term is appropriate, you will need to collect some data at another value of x1.

---------------------------------------------------------------------------------------------------

Thank you very much Mr. Lane, for the comment. I certainly need to have the x1^3 term - I am not allowed to modify the equation - it is the AHRI 540-2004 Standard. And, the data is indistinct because it is a measurement of evaporator temperature (at the suction of a compressor) and that is how compressor data is taken and performance evaluated. Either the evaporator or the condenser temperature will have to be held constant, while the other is measured. Is there any way I can get around that?

---------------------------------------------------------------------------------------------------

To Gautam: Tom has told you the essential difficulty with your problem as you have presented it. The difficulty is not with your ten-coefficient third degree polynomial. It is with the shortage of data to be fitted to it. The quantity you called 'x1' occurs in only three distinct values: -10, -5, and +5. That is just not enough, given the complexity of your polynomial.

Look at it this way. Suppose your problem involves only a single variable and is approximated with a third degree polynomial with four adjustable coefficients. It is obvious that three data values of that single variable would not be enough to uniquely determine the four coefficients. You would have only three equations and four unknowns, and there are infinitely many coefficient combinations that would give you an exact fit to that limited data.

In other words, Gautam, your data is too easy to fit as far as the x1 variable is concerned. You need much more variation in x1 to produce a meaningful set of ten coefficients. Actually you are also not providing enough variation in the other variable, x2. You have only four values of x2 for each value of x1 and that is insufficient to give enough information to determine the coefficients with sufficient reliability. It is too easy to fit. You really need much greater variation in both variables than you have provided. That is what is behind the messages you have received about "rank deficiency". Ideally the messages ought to have complained with the declaration: "Data! Data! Give me more data! I am starving from too little data!"

---------------------------------------------------------------------------------------------------

Thank you very much Mr. Stafford for clearing that up for me. I'm a newcomer to Matlab and could not grasp completely the meaning behind the warnings. I tried out fitlm with a much bigger data set and it did give me the exact same values EES did, which has around 2% variation with the published coefficients. I cannot thank you enough.

---------------------------------------------------------------------------------------------------

There are numerous ways to do this. The most easiest of them all is to use the Curve Fitting Tool with the correct options. Try this link for help:

http://www.mathworks.com/help/curvefit/polynomial.html

Using the Statistics Toolbox you can do that by specifying the polyijf modelspec:

http://www.mathworks.com/help/stats/fitlm.html#inputarg_modelspec

Another way is to solve a linear system in MATLAB. Assuming X, S and D are vectors in your MATLAB workspace. Create random data and compute the coefficients in C:

S = randn(100,1);
D = randn(100,1),
X = randn(100,1);

M = [ones(length(S),1), S, D, S.^2, S.*D, D.^2, S.^3, D.*S.^2, S.*D.^2, D.^3]
C = M\X

---------------------------------------------------------------------------------------------------

Shashank, thanks for the prompt response. In addition to the cftool I posted above, I also tried the below:

Using the same data (X and Y), I used the fitlm function with 'poly33' as the argument:

mdl = fitlm(X,Y,'poly33')

It did work, but the coefficients it gave me are not the coefficients I should be getting. The actual values of the coefficients are:

beta = [119.61 4.37 -0.39 -0.04 0.09 -0.01 0.00 0.02 0.00 0.00];

However, when I used fitlm, the value of the constant (119.61 above) is obtained as 0, which should not be the case.

As for the last method you mentioned (using the backslash operator), it gave me a warning that the rank of my matrix was deficient and still gave me wrong values for the values of the coefficients (and different from those obtained using the other functions).

Some more suggestions please?

---------------------------------------------------------------------------------------------------

I'm posting some of the attempts I made below:

X = [-10 101; -10 114; -10 129; -10 144; -5 96; -5 106; -5 119; -5 133; 5 92; 5 97; 5 108; 5 122];

>> Y = [72.03; 67.90; 62.77; 57.55; 88.53; 84.70; 79.91; 73.73; 123.43; 120.91; 115.44; 107.83];

>> x1 = X(:,1);

>> x2 = X(:,2);

>> mdl = fitlm(X,Y,'poly33')

Warning: Regression design matrix is rank deficient to within machine precision. > In TermsRegression>TermsRegression.checkDesignRank at 98 In LinearModel.LinearModel>LinearModel.fit at 868 In fitlm at 117

mdl =

Linear regression model: y ~ 1 + x1^2 + x1*x2 + x2^2 + x1^3 + (x1^2):x2 + x1:(x2^2) + x2^3

Estimated Coefficients: Estimate SE tStat pValue _________ ________ _______ ______

    (Intercept)              0             0          NaN         NaN
    x1                 -10.143        5.3752       -1.887     0.19979
    x2                 -0.3125        1.1432     -0.27335     0.81022
    x1^2                5.5413        1.6852       3.2882    0.081361
    x1:x2            0.0043522      0.026071      0.16694     0.88277
    x2^2            2.6105e-05      0.010348    0.0025228     0.99822
    x1^3               0.55158       0.16763       3.2905     0.08126
    x1^2:x2        -5.4944e-05    0.00021022     -0.26136     0.81827
    x1:x2^2        -8.3639e-05    0.00012017     -0.69599     0.55844
    x2^3             -4.12e-06     3.109e-05     -0.13252      0.9067

Number of observations: 12, Error degrees of freedom: 3 Root Mean Squared Error: 0.165 R-squared: 1, Adjusted R-Squared 1 F-statistic vs. constant model: 2.76e+04, p-value = 3.3e-07

Attempt 2:

 A = [1	-10	101	100	-1010	10201	-1000	10100	-102010	1030301
1	-5	96	25	-480	9216	-125	2400	-46080	884736
1	5	92	25	460	8464	125	2300	42320	778688
1	-10	114	100	-1140	12996	-1000	11400	-129960	1481544
1	-5	106	25	-530	11236	-125	2650	-56180	1191016
1	5	97	25	485	9409	125	2425	47045	912673
1	-10	129	100	-1290	16641	-1000	12900	-166410	2146689
1	-5	119	25	-595	14161	-125	2975	-70805	1685159
1	5	108	25	540	11664	125	2700	58320	1259712
1	-10	144	100	-1440	20736	-1000	14400	-207360	2985984
1	-5	133	25	-665	17689	-125	3325	-88445	2352637
1	5	122	25	610	14884	125	3050	74420	1815848
];

b = A\Y Warning: Rank deficient, rank = 9, tol = 4.463946e-09.

b =

         0
  -62.6154
   18.9421
  -35.6208
    2.9627
   -0.0958
   -3.8184
   -0.0236
   -0.0136
    0.0001

Attempt 3:

 sf = fit([x1 x2],Y, 'poly33', 'Robust', 'Bi')

Warning: Equation is badly conditioned. Remove repeated data points or try centering and scaling. > In Warning>Warning.throw at 30 In fit>iLinearFit at 690 In fit>iFit at 396 In fit at 108 Warning: Iteration limit reached for robust fitting. > In Warning>Warning.throw at 30 In fit>iRobustLinearFit at 796 In fit>iFit at 404 In fit at 108

 Linear model Poly33:
     sf(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2 + p30*x^3 + 
                    p21*x^2*y + p12*x*y^2 + p03*y^3
     Coefficients (with 95% confidence bounds):
       p00 =  -2.736e+14  (-6.205e+17, 6.199e+17)
       p10 =  -2.736e+13  (-6.205e+16, 6.199e+16)
       p01 =      0.3112  (-346.1, 346.7)
       p20 =   1.094e+13  (-2.48e+16, 2.482e+16)
       p11 =   -0.009459  (-5.569, 5.551)
       p02 =   -0.005727  (-3.308, 3.296)
       p30 =   1.094e+12  (-2.48e+15, 2.482e+15)
       p21 =  -0.0001193  (-0.4058, 0.4055)
       p12 =  -1.846e-05  (-0.02545, 0.02542)
       p03 =   1.358e-05  (-0.01007, 0.01009)

Attempt 4:

 sf = fit([x1 x2],Y, 'poly33', 'Robust', 'Bi','Normalize','on')
Warning: Equation is badly conditioned. Remove repeated data points or try centering and scaling. 
> In Warning>Warning.throw at 30
  In fit>iLinearFit at 690
  In fit>iFit at 396
  In fit at 108 
Warning: Iteration limit reached for robust fitting. 
> In Warning>Warning.throw at 30
  In fit>iRobustLinearFit at 796
  In fit>iFit at 404
  In fit at 108 

 Linear model Poly33:
     sf(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2 + p30*x^3 + 
                    p21*x^2*y + p12*x*y^2 + p03*y^3
       where x is normalized by mean -3.333 and std 6.513
       and where y is normalized by mean 113.4 and std 16.34
     Coefficients (with 95% confidence bounds):
       p00 =   2.802e+14  (-4.173e+15, 4.734e+15)
       p10 =    1.15e+15  (-1.712e+16, 1.942e+16)
       p01 =      -6.475  (-15.32, 2.371)
       p20 =       1.005  (-3.659, 5.67)
       p11 =      -1.507  (-6.289, 3.275)
       p02 =     -0.3495  (-5.084, 4.385)
       p30 =  -8.362e+14  (-1.413e+16, 1.245e+16)
       p21 =     -0.2157  (-8.185, 7.754)
       p12 =     -0.6978  (-12.13, 10.73)
       p03 =     -0.3148  (-7.039, 6.409)

Is there anything I haven't yet tried out or any mistakes I made? Anything at all would be helpful. Thank you.

---------------------------------------------------------------------------------------------------