From: <te...@us...> - 2009-03-24 17:40:51
|
Revision: 5636 http://octave.svn.sourceforge.net/octave/?rev=5636&view=rev Author: tealev Date: 2009-03-24 17:40:40 +0000 (Tue, 24 Mar 2009) Log Message: ----------- Confidence level of a sampled multinomial distribution. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX Added Paths: ----------- trunk/octave-forge/main/statistics/inst/cl_multinom.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2009-03-24 17:14:47 UTC (rev 5635) +++ trunk/octave-forge/main/statistics/INDEX 2009-03-24 17:40:40 UTC (rev 5636) @@ -25,6 +25,7 @@ jsucdf jsupdf random + cl_multinom Descriptive statistics nansum nanmax Added: trunk/octave-forge/main/statistics/inst/cl_multinom.m =================================================================== --- trunk/octave-forge/main/statistics/inst/cl_multinom.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/cl_multinom.m 2009-03-24 17:40:40 UTC (rev 5636) @@ -0,0 +1,124 @@ +# -*- texinfo -*- +# +# @deftypefn {Function File} {@var{CL}} cl_multinom( @var{x},@var{N},@var{b},@var{calculation_type} ) - Confidence level of multinomial portions +# Returns confidence level of multinomial parameters estimated @math{ p = x / sum(x) } with predefined confidence interval @var{b}. +# Finite population is also considered. +# +# This function calculates the level of confidence at which the samples represent the true distribution +# given that there is a predefined tolerance (confidence interval). +# This is the upside down case of the typical excercises at which we want to get the confidence interval +# given the confidence level (and the estimated parameters of the underlying distribution). +# But once we accept (lets say at elections) that we have a standard predefined +# maximal acceptable error rate (e.g. @var{b}=0.02 ) in the estimation and we just want to know that how sure we +# can be that the measured proportions are the same as in the +# entire population (ie. the expected value and mean of the samples are roghly the same) we need to use this function. +# +# @subheading Arguments +# @itemize @bullet +# @item @var{x} : int vector : sample frequencies bins +# @item @var{N} : int : Population size that was sampled by x. If N<sum(x), infinite number assumed +# @item @var{b} : real, vector : confidence interval +# if vector, it should be the size of x containing confence interval for each cells +# if scalar, each cell will have the same value of b unless it is zero or -1 +# if value is 0, b=.02 is assumed which is standard choice at elections +# otherwise it is calculated in a way that one sample in a cell alteration defines the confidence interval +# @item @var{calculation_type} : string : (Optional), described below +# "bromaghin" (default) - do not change it unless you have a good reason to do so +# "cochran" +# "agresti_cull" this is not exactly the solution at reference given below but an adjustment of the solutions above +# @end itemize +# +# @subheading Returns +# Confidence level. +# +# @subheading Example +# CL = cl_multinom( [27;43;19;11], 10000, 0.05 ) +# returns 0.69 confidence level. +# +# @subheading References +# +# "bromaghin" calculation type (default) is based on +# is based on the article +# Jeffrey F. Bromaghin, "Sample Size Determination for Interval Estimation of Multinomial Probabilities", The American Statistician vol 47, 1993, pp 203-206. +# +# "cochran" calculation type +# is based on article +# Robert T. Tortora, "A Note on Sample Size Estimation for Multinomial Populations", The American Statistician, , Vol 32. 1978, pp 100-102. +# +# "agresti_cull" calculation type +# is based on article in which Quesenberry Hurst and Goodman result is combined +# A. Agresti and B.A. Coull, "Approximate is better than \"exact\" for interval estimation of binomial portions", The American Statistician, Vol. 52, 1998, pp 119-126 +# +# @end deftypefn + + +# Copyright (C) 2009 Levente Torok / Tor...@gm... +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version.^ +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License" +# along with this program; If not, see <http://www.gnu.org/licenses/>. +# + +function CL = cl_multinom( x, N, b, calculation_type ) + + k = rows(x); + nn = sum(x); + p = x / nn; + + if (nargin < 3) + b = .05; + endif + if (isscalar( b )) + if (b==0) b=0.02; endif + b = ones( rows(x), 1 ) * b; + + if (b<0) b=1 ./ max( x, 1 ); endif + endif + bb = b .* b; + + if (N==nn) + CL = 1; + return; + endif + + if (N<nn) + fpc = 1; + else + fpc = (N-1) / (N-nn); # finite population correction tag + endif + + beta = p.*(1-p); + + if ( nargin < 4 ) + calculation_type = "bromaghin"; + endif + + switch calculation_type + case {"cochran"} + t = sqrt( fpc * nn * bb ./ beta ) + alpha = ( 1 - normcdf( t )) * 2 + + case {"bromaghin"} + t = sqrt( fpc * (nn * 2 * bb )./ ( beta - 2 * bb + sqrt( beta .* beta - bb .* ( 4*beta - 1 ))) ); + alpha = ( 1 - normcdf( t )) * 2; + + case {"agresti_cull"} + ts = fpc * nn * bb ./ beta ; + if ( k<=2 ) + alpha = 1 - chi2cdf( ts, k-1 ); % adjusted Wilson interval + else + alpha = 1 - chi2cdf( ts/k, 1 ); % Goodman interval with Bonferroni argument + endif + endswitch + + CL = 1 - max( alpha ); + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2009-09-30 11:16:59
|
Revision: 6283 http://octave.svn.sourceforge.net/octave/?rev=6283&view=rev Author: hauberg Date: 2009-09-30 11:16:48 +0000 (Wed, 30 Sep 2009) Log Message: ----------- Support Von Mises distributions Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX Added Paths: ----------- trunk/octave-forge/main/statistics/inst/vmpdf.m trunk/octave-forge/main/statistics/inst/vmrnd.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2009-09-30 08:18:02 UTC (rev 6282) +++ trunk/octave-forge/main/statistics/INDEX 2009-09-30 11:16:48 UTC (rev 6283) @@ -12,6 +12,7 @@ hygestat lognstat mvnpdf mvnrnd mvncdf + vmpdf vmrnd mvtcdf normstat nbinstat Added: trunk/octave-forge/main/statistics/inst/vmpdf.m =================================================================== --- trunk/octave-forge/main/statistics/inst/vmpdf.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/vmpdf.m 2009-09-30 11:16:48 UTC (rev 6283) @@ -0,0 +1,46 @@ +## Copyright (C) 2009 Soren Hauberg +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} @var{theta} = vmpdf (@var{x}, @var{mu}, @var{k}) +## Evaluates the Von Mises probability density function. +## +## The Von Mises distribution has probability density function +## @example +## f (@var{x}) = exp (@var{k} * cos (@var{x} - @var{mu})) / @var{Z} , +## @end example +## where @var{Z} is a normalisation constant. By default, @var{mu} is 0 and +## @var{k} is 1. +## @seealso{vmrnd} +## @end deftypefn + +function p = vmpdf (x, mu = 0, k = 1) + ## Check input + if (!isreal (x)) + error ("vmpdf: first input must be real"); + endif + + if (!isreal (mu)) + error ("vmpdf: second input must be a scalar"); + endif + + if (!isreal (k) || k <= 0) + error ("vmpdf: third input must be a real positive scalar"); + endif + + ## Evaluate PDF + Z = 2 * pi * besseli (0, k); + p = exp (k * cos (x-mu)) / Z; +endfunction Added: trunk/octave-forge/main/statistics/inst/vmrnd.m =================================================================== --- trunk/octave-forge/main/statistics/inst/vmrnd.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/vmrnd.m 2009-09-30 11:16:48 UTC (rev 6283) @@ -0,0 +1,76 @@ +## Copyright (C) 2009 Soren Hauberg +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} @var{theta} = vmrnd (@var{mu}, @var{k}) +## @deftypefnx{Function File} @var{theta} = vmrnd (@var{mu}, @var{k}, @var{sz}) +## Draw random angles from a Von Mises distribution with mean @var{mu} and +## concentration @var{k}. +## +## The Von Mises distribution has probability density function +## @example +## f (@var{x}) = exp (@var{k} * cos (@var{x} - @var{mu})) / @var{Z} , +## @end example +## where @var{Z} is a normalisation constant. +## +## The output, @var{theta}, is a matrix of size @var{sz} containing random angles +## drawn from the given Von Mises distribution. By default, @var{mu} is 0 +## and @var{k} is 1. +## @seealso{vmpdf} +## @end deftypefn + +function theta = vmrnd (mu = 0, k = 1, sz = 1) + ## Check input + if (!isreal (mu)) + error ("vmrnd: first input must be a scalar"); + endif + + if (!isreal (k) || k <= 0) + error ("vmrnd: second input must be a real positive scalar"); + endif + + if (isscalar (sz)) + sz = [sz, sz]; + elseif (!isvector (sz)) + error ("vmrnd: third input must be a scalar or a vector"); + endif + + ## Simulate! + if (k < 1e-6) + ## k is small: sample uniformly on circle + theta = 2 * pi * rand (sz) - pi; + + else + a = 1 + sqrt (1 + 4 * k.^2); + b = (a - sqrt (2 * a)) / (2 * k); + r = (1 + b^2) / (2 * b); + + N = prod (sz); + notdone = true (N, 1); + while (any (notdone)) + u (:, notdone) = rand (3, N); + + z (notdone) = cos (pi * u (1, notdone)); + f (notdone) = (1 + r * z (notdone)) ./ (r + z (notdone)); + c (notdone) = k * (r - f (notdone)); + + notdone = (u (2, :) >= c .* (2 - c)) & (log (c) - log (u (2, :)) + 1 - c < 0); + N = sum (notdone); + endwhile + + theta = mu + sign (u (3, :) - 0.5) .* acos (f); + theta = reshape (theta, sz); + endif +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2010-11-26 22:07:59
|
Revision: 7956 http://octave.svn.sourceforge.net/octave/?rev=7956&view=rev Author: hauberg Date: 2010-11-26 22:07:52 +0000 (Fri, 26 Nov 2010) Log Message: ----------- New function: combnk Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX Added Paths: ----------- trunk/octave-forge/main/statistics/inst/combnk.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2010-11-26 16:56:40 UTC (rev 7955) +++ trunk/octave-forge/main/statistics/INDEX 2010-11-26 22:07:52 UTC (rev 7956) @@ -42,6 +42,7 @@ trimmean zscore tabulate + combnk Experimental design fullfact ff2n Linear regression Added: trunk/octave-forge/main/statistics/inst/combnk.m =================================================================== --- trunk/octave-forge/main/statistics/inst/combnk.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/combnk.m 2010-11-26 22:07:52 UTC (rev 7956) @@ -0,0 +1,82 @@ +## Copyright (C) 2010 Soren Hauberg +## +## This program is free software: you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{c} =} combnk (@var{data}, @var{k}) +## Return all combinations of @var{k} elements in @var{data}. +## @end deftypefn + +function retval = combnk (data, k) + ## Check input + if (! isvector (data)) + error ("combnk: first input argument must be a vector"); + endif + + if (!isreal (k) || k != round (k) || k < 0) + error ("combnk: second input argument must be a non-negative integer"); + endif + + ## Simple checks + n = numel (data); + if (k == 0 || k > n) + retval = resize (data, 0, k); + elseif (k == n) + retval = data (:).'; + else + retval = __combnk__ (data, k); + endif + + ## For some odd reason Matlab seems to treat strings differently compared to other data-types... + if (ischar (data)) + retval = flipud (retval); + endif +endfunction + +function retval = __combnk__ (data, k) + ## Recursion stopping criteria + if (k == 1) + retval = data (:); + else + ## Process data + n = numel (data); + retval = []; + for j = 1:n + C = __combnk__ (data ((j+1):end), k-1); + C = cat (2, repmat (data (j), rows (C), 1), C); + retval = [retval; C]; + endfor + endif +endfunction + +%!demo +%! c = combnk (1:5, 2); +%! disp ("All pairs of integers between 1 and 5:"); +%! disp (c); + +%!test +%! c = combnk (1:3, 2); +%! assert (c, [1, 2; 1, 3; 2, 3]); + +%!test +%! c = combnk (1:3, 6); +%! assert (isempty (c)); + +%!test +%! c = combnk ({1, 2, 3}, 2); +%! assert (c, {1, 2; 1, 3; 2, 3}); + +%!test +%! c = combnk ("hello", 2); +%! assert (c, ["lo"; "lo"; "ll"; "eo"; "el"; "el"; "ho"; "hl"; "hl"; "he"]); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2011-11-14 14:47:10
|
Revision: 9085 http://octave.svn.sourceforge.net/octave/?rev=9085&view=rev Author: asnelt Date: 2011-11-14 14:47:03 +0000 (Mon, 14 Nov 2011) Log Message: ----------- Initial commit. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX Added Paths: ----------- trunk/octave-forge/main/statistics/inst/monotone_smooth.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2011-11-13 21:14:56 UTC (rev 9084) +++ trunk/octave-forge/main/statistics/INDEX 2011-11-14 14:47:03 UTC (rev 9085) @@ -45,8 +45,9 @@ combnk Experimental design fullfact ff2n -Linear regression +Regression anovan + monotone_smooth princomp regress Plots Added: trunk/octave-forge/main/statistics/inst/monotone_smooth.m =================================================================== --- trunk/octave-forge/main/statistics/inst/monotone_smooth.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/monotone_smooth.m 2011-11-14 14:47:03 UTC (rev 9085) @@ -0,0 +1,147 @@ +## Copyright (C) 2011 Nir Krakauer +## +## This program is free software: you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{yy} =} monotone_smooth (@var{x}, @var{y}, @var{h}) +## Produce a smooth monotone increasing approximation to a sampled functional dependence y(x) using a kernel method. +## (An Epanechnikov smoothing kernel is applied to y(x); this is integrated to yield the monotone increasing form. See Reference 1 for details.) +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{x} is a vector of values of the independent variable. +## +## @item +## @var{y} is a vector of values of the dependent variable, of the same size as @var{x}. For best performance, it is recommended that the @var{y} already be fairly smooth, e.g. by applying a kernel smoothing to the original values if they are noisy. +## +## @item +## @var{h} is the kernel bandwidth to use. If @var{h} is not given, a "reasonable" value is computed. +## +## @end itemize +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{yy} is the vector of smooth monotone increasing function values at @var{x}. +## +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## x = 0:0.1:10; +## y = (x .^ 2) + 3 * randn(size(x)); %typically non-monotonic from the added noise +## ys = ([y(1) y(1:(end-1))] + y + [y(2:end) y(end)])/3; %crudely smoothed via moving average, but still typically non-monotonic +## yy = monotone_smooth(x, ys); %yy is monotone increasing in x +## plot(x, y, '+', x, ys, x, yy) +## @end group +## +## @end example +## +## @subheading References +## +## @enumerate +## @item +## Holger Dette, Natalie Neumeyer and Kay F. Pilz (2006), A simple nonparametric estimator of a strictly monotone regression function, @cite{Bernoulli}, 12:469-490 +## @item +## Regine Scheder (2007), R Package 'monoProc', Version 1.0-6, @url{http://cran.r-project.org/web/packages/monoProc/monoProc.pdf} (The implementation here is based on the monoProc function mono.1d) +## @end enumerate +## @end deftypefn + +## Author: Nir Krakauer <nkr...@cc...> +## Description: Nonparametric monotone increasing regression + +function yy = monotone_smooth(x, y, h) + +if nargin < 2 || ~isnumeric(x) || ~isnumeric(y) + print_usage (); +endif + +n = numel(x); + +%reshape x to be a vector +x = x(:); + +%set filter bandwidth at a reasonable default value, if not specified +if nargin < 3 || ~isscalar(h) + s = std(x); + h = s / (n^0.2); +end + +x_min = min(x); +x_max = max(x); + +y_min = min(y); +y_max = max(y); + +%transform range of x to [0, 1] +xl = (x - x_min) / (x_max - x_min); + +yy = ones(size(y)); + +%Epanechnikov smoothing kernel (with finite support) +%K_epanech_kernel = @(z) (3/4) * ((1 - z).^2) .* (abs(z) < 1); + + +K_epanech_int = @(z) mean(((abs(z) < 1)/2) - (3/4) * (z .* (abs(z) < 1) - (1/3) * (z.^3) .* (abs(z) < 1)) + (z < -1)); + +%integral of kernels up to t +monotone_inverse = @(t) K_epanech_int((y - t) / h); + + + +%find the value of the monotone smooth function at each point in x +niter_max = 150; %maximum number of iterations for estimating each value (should not be reached in most cases) +for l = 1:n + + tmax = y_max; + tmin = y_min; + wmin = monotone_inverse(tmin); + wmax = monotone_inverse(tmax); + if (wmax == wmin) + yy(l) = tmin; + else + wt = xl(l); + iter_max_reached = 1; + for i = 1:niter_max + wt_scaled = (wt - wmin) / (wmax - wmin); + tn = tmin + wt_scaled * (tmax - tmin) ; + wn = monotone_inverse(tn); + wn_scaled = (wn - wmin) / (wmax - wmin); + %if (abs(wt-wn) < 1E-4) || (tn < (y_min-0.1)) || (tn > (y_max+0.1)) %%criterion for break in the R code -- replaced by the following line to hopefully be less dependent on the scale of y + if (abs(wt_scaled-wn_scaled) < 1E-4) || (wt_scaled < -0.1) || (wt_scaled > 1.1) + iter_max_reached = 0; + break + end + if wn > wt + tmax = tn; + wmax = wn; + else + tmin = tn; + wmin = wn; + end + end + if iter_max_reached + warning(['at x = ' num2str(x(l)) ', maximum number of iterations reached without convergence; approximation may not be optimal']) + end + yy(l) = tmin + (wt - wmin) * (tmax - tmin) / (wmax - wmin); + end +end + + + This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2011-12-16 20:48:23
|
Revision: 9416 http://octave.svn.sourceforge.net/octave/?rev=9416&view=rev Author: asnelt Date: 2011-12-16 20:48:17 +0000 (Fri, 16 Dec 2011) Log Message: ----------- Applied bugfix by Kyle Winfree to repanova. Modified Paths: -------------- trunk/octave-forge/main/statistics/NEWS trunk/octave-forge/main/statistics/inst/repanova.m Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2011-12-16 20:33:08 UTC (rev 9415) +++ trunk/octave-forge/main/statistics/NEWS 2011-12-16 20:48:17 UTC (rev 9416) @@ -5,7 +5,10 @@ monotone_smooth + ** Bug fixes on the function: + repanova + Summary of important user-visible changes for statistics 1.1.0: ------------------------------------------------------------------- Modified: trunk/octave-forge/main/statistics/inst/repanova.m =================================================================== --- trunk/octave-forge/main/statistics/inst/repanova.m 2011-12-16 20:33:08 UTC (rev 9415) +++ trunk/octave-forge/main/statistics/inst/repanova.m 2011-12-16 20:48:17 UTC (rev 9416) @@ -33,60 +33,68 @@ function [p, table, st] = repanova(varargin) - if (nargin == 2) - X = varargin{1}; - condition = varargin{2}; - option = 'cell'; - elseif (nargin == 3) - X = varargin{1}; - condition = varargin{2}; - option = varargin{3}; - else - print_usage; - endif +switch nargin + case 0 + error('Too few inputs.'); + case 1 + X = varargin{1}; + for c = 1:size(X, 2) + condition{c} = ['time', num2str(c)]; + end + option = 'cell'; + case 2 + X = varargin{1}; + condition = varargin{2}; + option = 'cell'; + case 3 + X = varargin{1}; + condition = varargin{2}; + option = varargin{3}; + otherwise + error('Too many inputs.'); +end + % Find the means of the subjects and measures, ignoring any NaNs + u_subjects = nanmean(X,2); + u_measures = nanmean(X,1); + u_grand = nansum(nansum(X)) / (size(X,1) * size(X,2)); + % Differences between rows will be reflected in SS subjects, differences + % between columns will be reflected in SS_within subjects. + N = size(X,1); % number of subjects + J = size(X,2); % number of samples per subject + SS_measures = N * nansum((u_measures - u_grand).^2); + SS_subjects = J * nansum((u_subjects - u_grand).^2); + SS_total = nansum(nansum((X - u_grand).^2)); + SS_error = SS_total - SS_measures - SS_subjects; + df_measures = J - 1; + df_subjects = N - 1; + df_grand = (N*J) - 1; + df_error = df_grand - df_measures - df_subjects; + MS_measures = SS_measures / df_measures; + MS_subjects = SS_subjects / df_subjects; + MS_error = SS_error / df_error; % variation expected as a result of sampling error alone + F = MS_measures / MS_error; + p = 1 - fcdf(F, df_measures, df_error); % Probability of F given equal means. - % Find the means of the subjects and measures, ignoring any NaNs - u_subjects = nanmean(X,2); - u_measures = nanmean(X,1); - u_grand = nansum(nansum(X)) / (size(X,1) * size(X,2)); - % Differences between rows will be reflected in SS subjects, differences - % between columns will be reflected in SS_within subjects. - N = size(X,1); % number of subjects - J = size(X,2); % number of samples per subject - SS_measures = N * nansum((u_measures - u_grand).^2); - SS_subjects = J * nansum((u_subjects - u_grand).^2); - SS_total = nansum(nansum((X - u_grand).^2)); - SS_error = SS_total - SS_measures - SS_subjects; - df_measures = J - 1; - df_subjects = N - 1; - df_grand = (N*J) - 1; - df_error = df_grand - df_measures - df_subjects; - MS_measures = SS_measures / df_measures; - MS_subjects = SS_subjects / df_subjects; - MS_error = SS_error / df_error; % variation expected as a result of sampling error alone - F = MS_measures / MS_error; - p = 1 - fcdf(F,df_measures,df_subjects); % Probability of F given equal means. + if strcmp(option, 'string') + table = [sprintf('\nSource\tSS\tdf\tMS\tF\tProb > F'), ... + sprintf('\nSubject\t%g\t%i\t%g', SS_subjects, df_subjects, MS_subjects), ... + sprintf('\nMeasure\t%g\t%i\t%g\t%g\t%g', SS_measures, df_measures, MS_measures, F, p), ... + sprintf('\nError\t%g\t%i\t%g', SS_error, df_error, MS_error), ... + sprintf('\n')]; + else + table = {'Source', 'Partial SS', 'df', 'MS', 'F', 'Prob > F'; ... + 'Subject', SS_subjects, df_subjects, MS_subjects, '', ''; ... + 'Measure', SS_measures, df_measures, MS_measures, F, p}; + end - if strcmp(option, 'string') - table = [sprintf('\nSource\tSS\tdf\tMS\tF\tProb > F'), ... - sprintf('\nSubject\t%g\t%i\t%g', SS_subjects, df_subjects, MS_subjects), ... - sprintf('\nMeasure\t%g\t%i\t%g\t%g\t%g', SS_measures, df_measures, MS_measures, F, p), ... - sprintf('\nError\t%g\t%i\t%g', SS_error, df_error, MS_error), ... - sprintf('\n')]; - else - table = {'Source', 'Partial SS', 'df', 'MS', 'F', 'Prob > F'; ... - 'Subject', SS_subjects, df_subjects, MS_subjects, '', ''; ... - 'Measure', SS_measures, df_measures, MS_measures, F, p}; - end + st.gnames = condition'; % this is the same struct format used in anova1 + st.n = repmat(N, 1, J); + st.source = 'anova1'; % it cannot be assumed that 'repanova' is a supported source for multcompare + st.means = u_measures; + st.df = df_error; + st.s = sqrt(MS_error); +end - st.gnames = condition'; % this is the same struct format used in anova1 - st.n = repmat(N, 1, J); - st.source = 'anova1'; % it cannot be assumed that 'repanova' is a supported source for multcompare - st.means = u_measures; - st.df = df_error; - st.s = sqrt(MS_error); -endfunction - % This function was created with guidance from the following websites: % http://courses.washington.edu/stat217/rmANOVA.html % http://grants.hhp.coe.uh.edu/doconnor/PEP6305/Topic%20010%20Repeated%20Measures.htm This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <car...@us...> - 2012-04-17 14:16:29
|
Revision: 10266 http://octave.svn.sourceforge.net/octave/?rev=10266&view=rev Author: carandraug Date: 2012-04-17 14:16:23 +0000 (Tue, 17 Apr 2012) Log Message: ----------- statistics: removed dependency from miscellaneous package (only dependended because of cell2csv and csv2cell now in io) Modified Paths: -------------- trunk/octave-forge/main/statistics/DESCRIPTION trunk/octave-forge/main/statistics/NEWS Modified: trunk/octave-forge/main/statistics/DESCRIPTION =================================================================== --- trunk/octave-forge/main/statistics/DESCRIPTION 2012-04-17 14:08:58 UTC (rev 10265) +++ trunk/octave-forge/main/statistics/DESCRIPTION 2012-04-17 14:16:23 UTC (rev 10266) @@ -6,7 +6,7 @@ Title: Statistics Description: Additional statistics functions for Octave. Categories: Statistics -Depends: octave (>= 2.9.7), miscellaneous, io (>= 1.0.18) +Depends: octave (>= 2.9.7), io (>= 1.0.18) Autoload: yes License: GPLv3+, public domain Url: http://octave.sf.net Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-04-17 14:08:58 UTC (rev 10265) +++ trunk/octave-forge/main/statistics/NEWS 2012-04-17 14:16:23 UTC (rev 10266) @@ -2,9 +2,8 @@ ------------------------------------------------------------------- ** The package is now dependent on the io package (version 1.0.18 or - later) since some of the functions that were previously part of - the miscellaneous package (already a dependency of statistics), - have been moved to io. + later) since the functions that it depended of from miscellaneous + package have been moved to io. Summary of important user-visible changes for statistics 1.1.1: ------------------------------------------------------------------- This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2012-05-09 19:54:11
|
Revision: 10391 http://octave.svn.sourceforge.net/octave/?rev=10391&view=rev Author: asnelt Date: 2012-05-09 19:54:04 +0000 (Wed, 09 May 2012) Log Message: ----------- Initial commit. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX Added Paths: ----------- trunk/octave-forge/main/statistics/inst/copularnd.m trunk/octave-forge/main/statistics/inst/mvtrnd.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2012-05-09 19:12:01 UTC (rev 10390) +++ trunk/octave-forge/main/statistics/INDEX 2012-05-09 19:54:04 UTC (rev 10391) @@ -5,7 +5,7 @@ binostat chi2stat cl_multinom - copulacdf copulapdf + copulacdf copulapdf copularnd expstat fstat gamlike @@ -17,7 +17,7 @@ lognstat mvnpdf mvnrnd mvncdf mnpdf mnrnd - mvtcdf + mvtcdf mvtrnd nbinstat normalise_distribution normstat Added: trunk/octave-forge/main/statistics/inst/copularnd.m =================================================================== --- trunk/octave-forge/main/statistics/inst/copularnd.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/copularnd.m 2012-05-09 19:54:04 UTC (rev 10391) @@ -0,0 +1,281 @@ +## Copyright (C) 2012 Arno Onken +## +## This program is free software: you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{x} =} copularnd (@var{family}, @var{theta}, @var{n}) +## @deftypefnx {Function File} {} copularnd (@var{family}, @var{theta}, @var{n}, @var{d}) +## @deftypefnx {Function File} {} copularnd ('t', @var{theta}, @var{nu}, @var{n}) +## Generate random samples from a copula family. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{family} is the copula family name. Currently, @var{family} can be +## @code{'Gaussian'} for the Gaussian family, @code{'t'} for the Student's t +## family, or @code{'Clayton'} for the Clayton family. +## +## @item +## @var{theta} is the parameter of the copula. For the Gaussian and Student's t +## copula, @var{theta} must be a correlation matrix. For bivariate copulas +## @var{theta} can also be a correlation coefficient. For the Clayton family, +## @var{theta} must be a vector with the same number of elements as samples to +## be generated or be scalar. +## +## @item +## @var{nu} is the degrees of freedom for the Student's t family. @var{nu} must +## be a vector with the same number of elements as samples to be generated or +## be scalar. +## +## @item +## @var{n} is the number of rows of the matrix to be generated. @var{n} must be +## a non-negative integer and corresponds to the number of samples to be +## generated. +## +## @item +## @var{d} is the number of columns of the matrix to be generated. @var{d} must +## be a positive integer and corresponds to the dimension of the copula. +## @end itemize +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{x} is a matrix of random samples from the copula with @var{n} samples +## of distribution dimension @var{d}. +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## theta = 0.5; +## x = copularnd ("Gaussian", theta); +## @end group +## +## @group +## theta = 0.5; +## nu = 2; +## x = copularnd ("t", theta, nu); +## @end group +## +## @group +## theta = 0.5; +## n = 2; +## x = copularnd ("Clayton", theta, n); +## @end group +## @end example +## +## @subheading References +## +## @enumerate +## @item +## Roger B. Nelsen. @cite{An Introduction to Copulas}. Springer, New York, +## second edition, 2006. +## @end enumerate +## @end deftypefn + +## Author: Arno Onken <as...@as...> +## Description: Random samples from a copula family + +function x = copularnd (family, theta, nu, n) + + # Check arguments + if (nargin < 2) + print_usage (); + endif + + if (! ischar (family)) + error ("copularnd: family must be one of 'Gaussian', 't', and 'Clayton'"); + endif + + lower_family = lower (family); + + # Check family and copula parameters + switch (lower_family) + + case {"gaussian"} + # Gaussian family + if (isscalar (theta)) + # Expand a scalar to a correlation matrix + theta = [1, theta; theta, 1]; + endif + if (! ismatrix (theta) || any (diag (theta) != 1) || any (any (theta != theta')) || min (eig (theta)) <= 0) + error ("copularnd: theta must be a correlation matrix"); + endif + if (nargin > 3) + d = n; + if (! isscalar (d) || d != size (theta, 1)) + error ("copularnd: d must correspond to dimension of theta"); + endif + else + d = size (theta, 1); + endif + if (nargin < 3) + n = 1; + else + n = nu; + if (! isscalar (n) || (n < 0) || round (n) != n) + error ("copularnd: n must be a non-negative integer"); + endif + endif + + case {"t"} + # Student's t family + if (nargin < 3) + print_usage (); + endif + if (isscalar (theta)) + # Expand a scalar to a correlation matrix + theta = [1, theta; theta, 1]; + endif + if (! ismatrix (theta) || any (diag (theta) != 1) || any (any (theta != theta')) || min (eig (theta)) <= 0) + error ("copularnd: theta must be a correlation matrix"); + endif + if (! isscalar (nu) && (! isvector (nu) || length (nu) != n)) + error ("copularnd: nu must be a vector with the same number of rows as x or be scalar"); + endif + nu = nu(:); + if (nargin < 4) + n = 1; + else + if (! isscalar (n) || (n < 0) || round (n) != n) + error ("copularnd: n must be a non-negative integer"); + endif + endif + + case {"clayton"} + # Archimedian one parameter family + if (nargin < 4) + # Default is bivariate + d = 2; + else + d = n; + if (! isscalar (d) || (d < 2) || round (d) != d) + error ("copularnd: d must be an integer greater than 1"); + endif + endif + if (nargin < 3) + # Default is one sample + n = 1; + else + n = nu; + if (! isscalar (n) || (n < 0) || round (n) != n) + error ("copularnd: n must be a non-negative integer"); + endif + endif + if (! isvector (theta) || (! isscalar (theta) && size (theta, 1) != n)) + error ("copularnd: theta must be a column vector with the number of rows equal to n or be scalar"); + endif + if (n > 1 && isscalar (theta)) + theta = repmat (theta, n, 1); + endif + + otherwise + error ("copularnd: unknown copula family '%s'", family); + + endswitch + + if (n == 0) + # Input is empty + x = zeros (0, d); + else + + # Draw random samples according to family + switch (lower_family) + + case {"gaussian"} + # The Gaussian family + x = normcdf (mvnrnd (zeros (1, d), theta, n), 0, 1); + # No parameter bounds check + k = []; + + case {"t"} + # The Student's t family + x = tcdf (mvtrnd (theta, nu, n), nu); + # No parameter bounds check + k = []; + + case {"clayton"} + # The Clayton family + u = rand (n, d); + if (d == 2) + x = zeros (n, 2); + # Conditional distribution method for the bivariate case which also + # works for theta < 0 + x(:, 1) = u(:, 1); + x(:, 2) = (1 + u(:, 1) .^ (-theta) .* (u(:, 2) .^ (-theta ./ (1 + theta)) - 1)) .^ (-1 ./ theta); + else + # Apply the algorithm by Marshall and Olkin: + # Frailty distribution for Clayton copula is gamma + y = randg (1 ./ theta, n, 1); + x = (1 - log (u) ./ repmat (y, 1, d)) .^ (-1 ./ repmat (theta, 1, d)); + endif + k = find (theta == 0); + if (any (k)) + # Produkt copula at columns k + x(k, :) = u(k, :); + endif + # Continue argument check + if (d == 2) + k = find (! (theta >= -1) | ! (theta < inf)); + else + k = find (! (theta >= 0) | ! (theta < inf)); + endif + + endswitch + + # Out of bounds parameters + if (any (k)) + x(k, :) = NaN; + endif + + endif + +endfunction + +%!test +%! theta = 0.5; +%! x = copularnd ("Gaussian", theta); +%! assert (size (x), [1, 2]); +%! assert (all ((x >= 0) & (x <= 1))); + +%!test +%! theta = 0.5; +%! nu = 2; +%! x = copularnd ("t", theta, nu); +%! assert (size (x), [1, 2]); +%! assert (all ((x >= 0) & (x <= 1))); + +%!test +%! theta = 0.5; +%! x = copularnd ("Clayton", theta); +%! assert (size (x), [1, 2]); +%! assert (all ((x >= 0) & (x <= 1))); + +%!test +%! theta = 0.5; +%! n = 2; +%! x = copularnd ("Clayton", theta, n); +%! assert (size (x), [n, 2]); +%! assert (all ((x >= 0) & (x <= 1))); + +%!test +%! theta = [1; 2]; +%! n = 2; +%! d = 3; +%! x = copularnd ("Clayton", theta, n, d); +%! assert (size (x), [n, d]); +%! assert (all ((x >= 0) & (x <= 1))); Added: trunk/octave-forge/main/statistics/inst/mvtrnd.m =================================================================== --- trunk/octave-forge/main/statistics/inst/mvtrnd.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/mvtrnd.m 2012-05-09 19:54:04 UTC (rev 10391) @@ -0,0 +1,137 @@ +## Copyright (C) 2012 Arno Onken <as...@as...> +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{x} =} mvtrnd (@var{sigma}, @var{nu}) +## @deftypefnx {Function File} {@var{x} =} mvtrnd (@var{sigma}, @var{nu}, @var{n}) +## Generate random samples from the multivariate t-distribution. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{sigma} is the matrix of correlation coefficients. If there are any +## non-unit diagonal elements then @var{sigma} will be normalized. +## +## @item +## @var{nu} is the degrees of freedom for the multivariate t-distribution. +## @var{nu} must be a vector with the same number of elements as samples to be +## generated or be scalar. +## +## @item +## @var{n} is the number of rows of the matrix to be generated. @var{n} must be +## a non-negative integer and corresponds to the number of samples to be +## generated. +## @end itemize +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{x} is a matrix of random samples from the multivariate t-distribution +## with @var{n} row samples. +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## sigma = [1, 0.5; 0.5, 1]; +## nu = 3; +## n = 10; +## x = mvtrnd (sigma, nu, n); +## @end group +## +## @group +## sigma = [1, 0.5; 0.5, 1]; +## nu = [2; 3]; +## n = 2; +## x = mvtrnd (sigma, nu, 2); +## @end group +## @end example +## +## @subheading References +## +## @enumerate +## @item +## Wendy L. Martinez and Angel R. Martinez. @cite{Computational Statistics +## Handbook with MATLAB}. Appendix E, pages 547-557, Chapman & Hall/CRC, 2001. +## +## @item +## Samuel Kotz and Saralees Nadarajah. @cite{Multivariate t Distributions and +## Their Applications}. Cambridge University Press, Cambridge, 2004. +## @end enumerate +## @end deftypefn + +## Author: Arno Onken <as...@as...> +## Description: Random samples from the multivariate t-distribution + +function x = mvtrnd (sigma, nu, n) + + # Check arguments + if (nargin < 2) + print_usage (); + endif + + if (! ismatrix (sigma) || any (any (sigma != sigma')) || min (eig (sigma)) <= 0) + error ("mvtrnd: sigma must be a positive definite matrix"); + endif + + if (!isvector (nu) || any (nu <= 0)) + error ("mvtrnd: nu must be a positive scalar or vector"); + endif + nu = nu(:); + + if (nargin > 2) + if (! isscalar (n) || n < 0 | round (n) != n) + error ("mvtrnd: n must be a non-negative integer") + endif + if (isscalar (nu)) + nu = nu * ones (n, 1); + else + if (length (nu) != n) + error ("mvtrnd: n must match the length of nu") + endif + endif + else + n = length (nu); + endif + + # Normalize sigma + if (any (diag (sigma) != 1)) + sigma = sigma ./ sqrt (diag (sigma) * diag (sigma)'); + endif + + # Dimension + d = size (sigma, 1); + # Draw samples + y = mvnrnd (zeros (1, d), sigma, n); + u = repmat (chi2rnd (nu), 1, d); + x = y .* sqrt (repmat (nu, 1, d) ./ u); +endfunction + +%!test +%! sigma = [1, 0.5; 0.5, 1]; +%! nu = 3; +%! n = 10; +%! x = mvtrnd (sigma, nu, n); +%! assert (size (x), [10, 2]); + +%!test +%! sigma = [1, 0.5; 0.5, 1]; +%! nu = [2; 3]; +%! n = 2; +%! x = mvtrnd (sigma, nu, 2); +%! assert (size (x), [2, 2]); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <jpi...@us...> - 2012-09-18 21:35:48
|
Revision: 11049 http://octave.svn.sourceforge.net/octave/?rev=11049&view=rev Author: jpicarbajal Date: 2012-09-18 21:35:42 +0000 (Tue, 18 Sep 2012) Log Message: ----------- statistics:adding dendogram function. Working but not finshed. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX trunk/octave-forge/main/statistics/NEWS trunk/octave-forge/main/statistics/inst/linkage.m Added Paths: ----------- trunk/octave-forge/main/statistics/inst/dendogram.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2012-09-18 18:43:56 UTC (rev 11048) +++ trunk/octave-forge/main/statistics/INDEX 2012-09-18 21:35:42 UTC (rev 11049) @@ -57,6 +57,7 @@ normplot histfit repanova + dendogram Models hmmestimate hmmgenerate hmmviterbi Hypothesis testing Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-09-18 18:43:56 UTC (rev 11048) +++ trunk/octave-forge/main/statistics/NEWS 2012-09-18 21:35:42 UTC (rev 11049) @@ -3,7 +3,7 @@ ** The following functions are new: - regress_gp + regress_gp dendogram ** The interface of the following functions has been modified: Added: trunk/octave-forge/main/statistics/inst/dendogram.m =================================================================== --- trunk/octave-forge/main/statistics/inst/dendogram.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/dendogram.m 2012-09-18 21:35:42 UTC (rev 11049) @@ -0,0 +1,114 @@ +%% Copyright (c) 2012 Juan Pablo Carbajal <car...@if...> +%% +%% This program is free software: you can redistribute it and/or modify +%% it under the terms of the GNU General Public License as published by +%% the Free Software Foundation, either version 3 of the License, or +%% any later version. +%% +%% This program is distributed in the hope that it will be useful, +%% but WITHOUT ANY WARRANTY; without even the implied warranty of +%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +%% GNU General Public License for more details. +%% +%% You should have received a copy of the GNU General Public License +%% along with this program. If not, see <http://www.gnu.org/licenses/>. + +%% -*- texinfo -*- +%% @deftypefn {Function File} {@var{h} = } dendogram (@var{tree}) +%% Plots a dendogram using the output of function @command{linkage}. +%% +%% TODO: This function needs your help to be finished! Leafs must be ordered to +%% prevent corssing of edges. +%% @seealso{linkage} +%% @end deftypefn + +## TODO +# Order leafs nodes to avoid edge crossing. + +function h = dendogram (tree) + + [m d] = size (tree); + if d != 3 + error ("Input data must be a tree as returned by function linkage.") + end + n = m + 1; + + nc = max(tree(:,1:2)(:)); + + % Vector with the horizontal and vertical position of each cluster + p = zeros (nc,2); + +##### This is an ugly hack and is a partial solution. If you know how to oder +##### the leafs do it. Please do not copy from dendogram.m of Mathworks! + + %% Order the leafs + % find which clusters contain th eleafs + tf = ismember (tree(:,1:2),1:n); + [idx_i,idx_j] = find (tf); + ind0 = find(tf(:)); + + % Assign a position between 1:n to the leafs according to clustering + % I am trying to avoid crossings. Better way? sure!. + i = 1; + j = 0; + while j < n + + % If the current leaf hasn't a position assigned, assign one. + if p(tree(idx_i(i),idx_j(i)),1) == 0 + j+=1; + p(tree(idx_i(i),idx_j(i)),1) = j; + end + + % Check if the current cluster has another leaf. If it has and the leaf + % doesn't have a position, assign one. + next = tree(idx_i(i),setdiff(1:2,idx_j(i))); + if next <= n && p(next,1) == 0 + j+=1; + p(next,1) = j; + end + + % Find the id of the current cluster and check if any other cluster + % that will be merged with this one has a leaf. + c_id = idx_i(i) + n; + [idx,~] = find (tree(:,1:2)==c_id); + idx2 = find (tree(idx,1:2)<=n); + if ~isempty (idx2) && p(tree(idx,idx2),1) == 0 + j+=1; + p(tree(idx,idx2),1) = j; + end + i+=1; + end + +##### End of the hack + + % Compute the horizontal position, begin-end + % and vertical position of all clusters. + for i = 1:m + p(n+i,1) = mean (p(tree(i,1:2),1)); + p(n+i,2) = tree(i,3); + x(i,1:2) = p(tree(i,1:2),1); + end + + + + clf + % plot horizontal lines + tmp = line (x', tree(:,[3 3])'); + + % plot vertical lines + for i=1:nc + [ind,~] = find (tree(:,1:2)==i); + tmp = line (p([i; i],1),[p(i,2); tree(ind,3)]); + end + + xticks = 1:n; + xl_txt = arrayfun (@num2str, tree(:,1:2)(ind0),"uniformoutput",false); + set (gca,"xticklabel",xl_txt,"xtick",xticks); + axis ([0.5 n+0.5 0 max(tree(:,3))+0.1*min(tree(:,3))]); + +endfunction + +%!demo +%! y = [4 5; 2 6; 3 7; 8 9; 1 10]; +%! y(:,3) = 1:5; +%! dendogram(y); Modified: trunk/octave-forge/main/statistics/inst/linkage.m =================================================================== --- trunk/octave-forge/main/statistics/inst/linkage.m 2012-09-18 18:43:56 UTC (rev 11048) +++ trunk/octave-forge/main/statistics/inst/linkage.m 2012-09-18 21:35:42 UTC (rev 11049) @@ -34,8 +34,8 @@ ## cluster and numbering those from 1 to n. Then it merges two ## clusters, chosen according to @var{method}, to create a new cluster ## numbered n+1, and so on until all observations are grouped into -## a single cluster numbered 2*n-1. Row m of the -## (m-1)x3 output matrix relates to cluster n+m: the first +## a single cluster numbered 2(n-1). Row k of the +## (m-1)x3 output matrix relates to cluster n+k: the first ## two columns are the numbers of the two component clusters and column ## 3 contains their distance. ## This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2011-12-22 13:22:15
|
Revision: 9456 http://octave.svn.sourceforge.net/octave/?rev=9456&view=rev Author: asnelt Date: 2011-12-22 13:22:04 +0000 (Thu, 22 Dec 2011) Log Message: ----------- Initial commit. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX trunk/octave-forge/main/statistics/NEWS Added Paths: ----------- trunk/octave-forge/main/statistics/inst/kmeans.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2011-12-21 21:33:00 UTC (rev 9455) +++ trunk/octave-forge/main/statistics/INDEX 2011-12-22 13:22:04 UTC (rev 9456) @@ -62,8 +62,9 @@ Fitting gamfit Clustering + kmeans + linkage pdist - linkage squareform Reading and Writing caseread Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2011-12-21 21:33:00 UTC (rev 9455) +++ trunk/octave-forge/main/statistics/NEWS 2011-12-22 13:22:04 UTC (rev 9456) @@ -1,9 +1,9 @@ Summary of important user-visible changes for statistics 1XXXXXXXX: ------------------------------------------------------------------- - ** The following functions are new in 1.1.0: + ** The following functions are new in 1.1.X: - monotone_smooth + monotone_smooth kmeans ** Bug fixes on the function: Added: trunk/octave-forge/main/statistics/inst/kmeans.m =================================================================== --- trunk/octave-forge/main/statistics/inst/kmeans.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/kmeans.m 2011-12-22 13:22:04 UTC (rev 9456) @@ -0,0 +1,84 @@ +## Copyright (C) 2011 Soren Hauberg +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, write to the Free Software +## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + +function [classes, centers] = kmeans (data, k, varargin) + ## Input checking + if (!ismatrix (data) || !isreal (data)) + error ("kmeans: first input argument must be a DxN real data matrix"); + endif + if (!isscalar (k)) + error ("kmeans: second input argument must be a scalar"); + endif + + [N, D] = size (data); + + ## (so far) Harcoded options + maxiter = Inf; + start = "sample"; + + ## Find initial clusters + switch (lower (start)) + case "sample" + idx = randperm (N) (1:k); + centers = data (idx, :); + otherwise + error ("kmeans: unsupported initial clustering parameter"); + endswitch + + ## Run the algorithm + D = zeros (N, k); + iterations = 0; + prevcenters = centers; + while (true) + ## Compute distances + for i = 1:k + D (:, i) = sum (( data - repmat (centers (i, :), N, 1)).^2, 2); + endfor + + ## Classify + [tmp, classes] = min (D, [], 2); + + ## Recompute centers + for i = 1:k + centers (i, :) = mean (data (classes == i, :)); + endfor + + ## Check for convergence + iterations++; + if (all (centers (:) == prevcenters (:)) || iterations >= maxiter) + break; + endif + prevcenters = centers; + endwhile +endfunction + +%!demo +%! ## Generate a two-cluster problem +%! C1 = randn (100, 2) + 1; +%! C2 = randn (100, 2) - 1; +%! data = [C1; C2]; +%! +%! ## Perform clustering +%! [idx, centers] = kmeans (data, 2); +%! +%! ## Plot the result +%! figure +%! plot (data (idx==1, 1), data (idx==1, 2), 'ro'); +%! hold on +%! plot (data (idx==2, 1), data (idx==2, 2), 'bs'); +%! plot (centers (:, 1), centers (:, 2), 'kv', 'markersize', 10); +%! hold off + This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2011-12-23 20:49:10
|
Revision: 9460 http://octave.svn.sourceforge.net/octave/?rev=9460&view=rev Author: asnelt Date: 2011-12-23 20:49:04 +0000 (Fri, 23 Dec 2011) Log Message: ----------- Initial commit. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX trunk/octave-forge/main/statistics/NEWS Added Paths: ----------- trunk/octave-forge/main/statistics/inst/jackknife.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2011-12-23 18:29:45 UTC (rev 9459) +++ trunk/octave-forge/main/statistics/INDEX 2011-12-23 20:49:04 UTC (rev 9460) @@ -43,6 +43,7 @@ zscore tabulate combnk + jackknife Experimental design fullfact ff2n Regression Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2011-12-23 18:29:45 UTC (rev 9459) +++ trunk/octave-forge/main/statistics/NEWS 2011-12-23 20:49:04 UTC (rev 9460) @@ -3,7 +3,7 @@ ** The following functions are new in 1.1.X: - monotone_smooth kmeans + monotone_smooth kmeans jackknife ** Bug fixes on the function: Added: trunk/octave-forge/main/statistics/inst/jackknife.m =================================================================== --- trunk/octave-forge/main/statistics/inst/jackknife.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/jackknife.m 2011-12-23 20:49:04 UTC (rev 9460) @@ -0,0 +1,141 @@ +## Copyright (C) 2011 Alexander Klein +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} { [ @var{t}, @var{v} ] } = jackknife ( @var{E}, @var{x}, ... ) +## Compute jackknife estimates of a parameter taking one or more given samples as parameters, +## as well as its variance. +## In particular, @var{E} is the estimator to be jackknifed as a function name, handle, +## or inline function, and @var{x} is the sample for which the estimate is to be taken. +## @var{t} will be the estimate for the parameter, and @var{v} the estimate of its variance. +## +## Depending on the number of samples to be used, the estimator must have the appropriate form: +## If only one sample is used, then the estimator need not be concerned with cell arrays, +## for example jackknifing the standard deviation of a sample can be performed with +## @code{ [ @var{t}, @var{v} ] = jackknife(@@std, rand (100, 1))}. +## If, however, more than one sample is to be used, the samples must all be of equal size, +## and the estimator must address them as elements of a cell-array, +## in which they are aggregated in their order of appearance: +## @code{ [ @var{t}, @var{v} ] = jackknife(@@(x) std(x@{1@})/var(x@{2@}), rand (100, 1), randn (100, 1)}. +## +## If all goes well, a theoretical value @var{P} for the parameter is already known, +## and @var{n} is the sample size, then +## @code{ (@var{t}-@var{P})/sqrt(@var{v}) } should follow a t-distribution with @var{n}-1 degrees of freedom. +## +## Jackknifing is a well known method to reduce bias; further details can be found in: +## @itemize @bullet +## @item Rupert G. Miller: The jackknife-a review; Biometrika (1974) 61(1): 1-15; doi:10.1093/biomet/61.1.1 +## @item Rupert G. Miller: Jackknifing Variances; Ann. Math. Statist. Volume 39, Number 2 (1968), 567-582; doi:10.1214/aoms/1177698418 +## @item M. H. Quenouille: Notes on Bias in Estimation; Biometrika Vol. 43, No. 3/4 (Dec., 1956), pp. 353-360; doi:10.1093/biomet/43.3-4.353 +## @end itemize +## @end deftypefn + +## Author: Alexander Klein <ale...@ma...> +## Created: 2011-11-25 + +function [ theta_tilde, var_theta_tilde ] = jackknife ( anEstimator, varargin ) + + + ## Convert function name to handle if necessary, or throw + ## an error. + if ( !strcmp ( typeinfo ( anEstimator ), "function handle" ) ) + + if ( isascii ( anEstimator ) ) + + anEstimator = str2func ( anEstimator ); + + else + + error ( "Estimators must be passed as function names or handles!" ); + end + end + + + ## Simple jackknifing can be done with a single vector argument, and + ## first and foremost with a function that does not care about + ## cell-arrays. + if ( length ( varargin ) == 1 && isnumeric ( varargin { 1 } ) ) + + aSample = varargin { 1 }; + + g = length ( aSample ); + indices = 1 : g; + + theta_hat = anEstimator ( aSample ); + + theta_hat_minus_i = zeros ( 1, g ); + + for k = 1 : g + theta_hat_minus_i ( k ) = anEstimator ( aSample ( [ 1 : k - 1, k + 1 : g ] ) ); + end + + theta_tilde = g * theta_hat - ( g - 1 ) * mean ( theta_hat_minus_i ); + + var_theta_tilde = sumsq ( ( g * theta_hat - ( g - 1 ) * theta_hat_minus_i ) - theta_tilde ) / ( g * ( g - 1 ) ); + + ## More complicated input requires more work, however. + else + + g = cellfun ( @(x) length ( x ), varargin ); + + if ( any ( g - g ( 1 ) ) ) + + error ( "All passed data must be of equal length!" ); + end + + g = g ( 1 ); + + theta_hat = anEstimator ( varargin ); + + theta_hat_minus_i = zeros ( 1, g ); + + flags = 1 - eye ( 1, g ); + + for k = 1 : g + + theta_hat_minus_i ( k ) = anEstimator ( cellfun ( @(x) x( [ 1 : k - 1, k + 1 : g ] ), varargin, "UniformOutput", false ) ); + end + + theta_tilde = g * theta_hat - ( g - 1 ) * mean ( theta_hat_minus_i ); + + var_theta_tilde = sumsq ( ( g * theta_hat - ( g - 1 ) * theta_hat_minus_i ) - theta_tilde ) / ( g * ( g - 1 ) ); + + end +endfunction + + +%!test +%! ##Example from Quenouille, Table 1 +%! d=[0.18 4.00 1.04 0.85 2.14 1.01 3.01 2.33 1.57 2.19]; +%! t = jackknife ( @(x) 1/mean(x), d ); +%! assert ( t, 0.5240, 1e-5 ); + +%!demo +%! for k = 1:1000 +%! x=rand(10,1); +%! s(k)=std(x); +%! j(k)=jackknife(@std,x); +%! end +%! figure();hist([s',j'], 0:sqrt(1/12)/10:2*sqrt(1/12)) + +%!demo +%! for k = 1:1000 +%! x=randn(1,50); +%! y=rand(1,50); +%! [j(k),v(k)]=jackknife(@(x) std(x{1})/std(x{2}),y,x); +%! end +%! t=(j-sqrt(1/12))./sqrt(v); +%! figure();plot(sort(tcdf(t,49)),"-;Almost linear mapping indicates good fit with t-distribution.;") This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2012-01-19 18:13:22
|
Revision: 9545 http://octave.svn.sourceforge.net/octave/?rev=9545&view=rev Author: asnelt Date: 2012-01-19 18:13:16 +0000 (Thu, 19 Jan 2012) Log Message: ----------- Removed zscore function since it is now part of GNU octave core. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX trunk/octave-forge/main/statistics/NEWS Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2012-01-19 16:21:36 UTC (rev 9544) +++ trunk/octave-forge/main/statistics/INDEX 2012-01-19 18:13:16 UTC (rev 9545) @@ -40,7 +40,6 @@ harmmean mad trimmean - zscore tabulate combnk jackknife Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-01-19 16:21:36 UTC (rev 9544) +++ trunk/octave-forge/main/statistics/NEWS 2012-01-19 18:13:16 UTC (rev 9545) @@ -9,6 +9,11 @@ repanova + ** The following functions were removed since equivalents are now + part of GNU octave core: + + zscore + Summary of important user-visible changes for statistics 1.1.0: ------------------------------------------------------------------- This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <car...@us...> - 2012-03-14 22:29:07
|
Revision: 9887 http://octave.svn.sourceforge.net/octave/?rev=9887&view=rev Author: carandraug Date: 2012-03-14 22:29:00 +0000 (Wed, 14 Mar 2012) Log Message: ----------- normalise_distribution: workaround due to octave bug #34765 was made for releases 3.4.X only. However, regression in octave is not yet fixed as of 3.6.1. Expanding workaround for all releases after 3.4.0 until fixed in core. Modified Paths: -------------- trunk/octave-forge/main/statistics/NEWS trunk/octave-forge/main/statistics/inst/normalise_distribution.m Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-03-14 22:17:16 UTC (rev 9886) +++ trunk/octave-forge/main/statistics/NEWS 2012-03-14 22:29:00 UTC (rev 9887) @@ -7,6 +7,7 @@ ** Bug fixes on the function: + normalise_distribution repanova ** The following functions were removed since equivalents are now Modified: trunk/octave-forge/main/statistics/inst/normalise_distribution.m =================================================================== --- trunk/octave-forge/main/statistics/inst/normalise_distribution.m 2012-03-14 22:17:16 UTC (rev 9886) +++ trunk/octave-forge/main/statistics/inst/normalise_distribution.m 2012-03-14 22:29:00 UTC (rev 9887) @@ -1,23 +1,17 @@ -## Copyright (C) 2011 Alexander Klein -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. +## Copyright (C) 2011 Alexander Klein <ale...@ma...> ## -## Author: Alexander Klein <ale...@ma...> -## Created: 2011-09-13 +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## TODO: +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn{Function File} {@var{NORMALISED} =} normalise_distribution (@var{DATA}) @@ -198,7 +192,11 @@ ## differently depending on octave version. This applies the fix for all ## 3.4 releases but it probably only appeared on 3.4.3 (can someone check?) ## See https://savannah.gnu.org/bugs/index.php?34765 - if (compare_versions (OCTAVE_VERSION, "3.4", "<") || compare_versions (OCTAVE_VERSION, "3.6", ">")) + ## Turns out that the bug was not completely fixed the first time and is + ## still present in 3.6.1 + ## FIXME Once package dependency increases beyond an octave version that + ## has this fixed, remove this + if (compare_versions (OCTAVE_VERSION, "3.4", "<")) ## this is how it should work f_remap = @( k ) ( normal ( k ) ); normalised ( :, k ) = arrayfun ( f_remap, target_indices ); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <car...@us...> - 2012-03-14 22:44:12
|
Revision: 9888 http://octave.svn.sourceforge.net/octave/?rev=9888&view=rev Author: carandraug Date: 2012-03-14 22:44:06 +0000 (Wed, 14 Mar 2012) Log Message: ----------- combnk: fixed bug due to changes (probably bug fix) in octave-core. Handle cell input correctly Modified Paths: -------------- trunk/octave-forge/main/statistics/NEWS trunk/octave-forge/main/statistics/inst/combnk.m Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-03-14 22:29:00 UTC (rev 9887) +++ trunk/octave-forge/main/statistics/NEWS 2012-03-14 22:44:06 UTC (rev 9888) @@ -5,9 +5,9 @@ monotone_smooth kmeans jackknife - ** Bug fixes on the function: + ** Bug fixes on the functions: - normalise_distribution + normalise_distribution combnk repanova ** The following functions were removed since equivalents are now Modified: trunk/octave-forge/main/statistics/inst/combnk.m =================================================================== --- trunk/octave-forge/main/statistics/inst/combnk.m 2012-03-14 22:29:00 UTC (rev 9887) +++ trunk/octave-forge/main/statistics/inst/combnk.m 2012-03-14 22:44:06 UTC (rev 9888) @@ -1,17 +1,17 @@ -## Copyright (C) 2010 Soren Hauberg +## Copyright (C) 2010 Soren Hauberg <so...@ha...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {@var{c} =} combnk (@var{data}, @var{k}) @@ -51,7 +51,11 @@ else ## Process data n = numel (data); - retval = []; + if iscell (data) + retval = {}; + else + retval = []; + endif for j = 1:n C = __combnk__ (data ((j+1):end), k-1); C = cat (2, repmat (data (j), rows (C), 1), C); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2012-04-12 10:13:35
|
Revision: 10202 http://octave.svn.sourceforge.net/octave/?rev=10202&view=rev Author: asnelt Date: 2012-04-12 10:13:25 +0000 (Thu, 12 Apr 2012) Log Message: ----------- Updates DESCRIPTION and NEWS for new package release. Modified Paths: -------------- trunk/octave-forge/main/statistics/DESCRIPTION trunk/octave-forge/main/statistics/NEWS Modified: trunk/octave-forge/main/statistics/DESCRIPTION =================================================================== --- trunk/octave-forge/main/statistics/DESCRIPTION 2012-04-12 09:49:33 UTC (rev 10201) +++ trunk/octave-forge/main/statistics/DESCRIPTION 2012-04-12 10:13:25 UTC (rev 10202) @@ -1,6 +1,6 @@ Name: Statistics -Version: 1.1.0 -Date: 2011-11-10 +Version: 1.1.1 +Date: 2012-04-12 Author: Various Authors Maintainer: Arno Onken <as...@as...> Title: Statistics Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-04-12 09:49:33 UTC (rev 10201) +++ trunk/octave-forge/main/statistics/NEWS 2012-04-12 10:13:25 UTC (rev 10202) @@ -1,7 +1,7 @@ -Summary of important user-visible changes for statistics 1XXXXXXXX: +Summary of important user-visible changes for statistics 1.1.1: ------------------------------------------------------------------- - ** The following functions are new in 1.1.X: + ** The following functions are new in 1.1.1: monotone_smooth kmeans jackknife @@ -15,7 +15,7 @@ zscore - ** boxplot.m now retunrs a structure with handles to the plot elemenets. + ** boxplot.m now returns a structure with handles to the plot elemenets. Summary of important user-visible changes for statistics 1.1.0: ------------------------------------------------------------------- This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <car...@us...> - 2012-04-12 13:47:47
|
Revision: 10203 http://octave.svn.sourceforge.net/octave/?rev=10203&view=rev Author: carandraug Date: 2012-04-12 13:47:32 +0000 (Thu, 12 Apr 2012) Log Message: ----------- maint: adjust copyright headers on statistics package. Move them to top of file, update to GPLv3, add e-mail address, list them on DESCRIPTION file Modified Paths: -------------- trunk/octave-forge/main/statistics/COPYING trunk/octave-forge/main/statistics/DESCRIPTION trunk/octave-forge/main/statistics/inst/anderson_darling_cdf.m trunk/octave-forge/main/statistics/inst/anderson_darling_test.m trunk/octave-forge/main/statistics/inst/anovan.m trunk/octave-forge/main/statistics/inst/betastat.m trunk/octave-forge/main/statistics/inst/binostat.m trunk/octave-forge/main/statistics/inst/boxplot.m trunk/octave-forge/main/statistics/inst/caseread.m trunk/octave-forge/main/statistics/inst/casewrite.m trunk/octave-forge/main/statistics/inst/chi2stat.m trunk/octave-forge/main/statistics/inst/cl_multinom.m trunk/octave-forge/main/statistics/inst/copulacdf.m trunk/octave-forge/main/statistics/inst/copulapdf.m trunk/octave-forge/main/statistics/inst/expstat.m trunk/octave-forge/main/statistics/inst/ff2n.m trunk/octave-forge/main/statistics/inst/fstat.m trunk/octave-forge/main/statistics/inst/fullfact.m trunk/octave-forge/main/statistics/inst/gamfit.m trunk/octave-forge/main/statistics/inst/gamlike.m trunk/octave-forge/main/statistics/inst/gamstat.m trunk/octave-forge/main/statistics/inst/geomean.m trunk/octave-forge/main/statistics/inst/geostat.m trunk/octave-forge/main/statistics/inst/harmmean.m trunk/octave-forge/main/statistics/inst/histfit.m trunk/octave-forge/main/statistics/inst/hmmestimate.m trunk/octave-forge/main/statistics/inst/hmmgenerate.m trunk/octave-forge/main/statistics/inst/hmmviterbi.m trunk/octave-forge/main/statistics/inst/hygestat.m trunk/octave-forge/main/statistics/inst/jackknife.m trunk/octave-forge/main/statistics/inst/jsucdf.m trunk/octave-forge/main/statistics/inst/jsupdf.m trunk/octave-forge/main/statistics/inst/kmeans.m trunk/octave-forge/main/statistics/inst/linkage.m trunk/octave-forge/main/statistics/inst/lognstat.m trunk/octave-forge/main/statistics/inst/mad.m trunk/octave-forge/main/statistics/inst/monotone_smooth.m trunk/octave-forge/main/statistics/inst/mvncdf.m trunk/octave-forge/main/statistics/inst/mvnpdf.m trunk/octave-forge/main/statistics/inst/mvnrnd.m trunk/octave-forge/main/statistics/inst/mvtcdf.m trunk/octave-forge/main/statistics/inst/nanmax.m trunk/octave-forge/main/statistics/inst/nanmean.m trunk/octave-forge/main/statistics/inst/nanmedian.m trunk/octave-forge/main/statistics/inst/nanmin.m trunk/octave-forge/main/statistics/inst/nanstd.m trunk/octave-forge/main/statistics/inst/nansum.m trunk/octave-forge/main/statistics/inst/nanvar.m trunk/octave-forge/main/statistics/inst/nbinstat.m trunk/octave-forge/main/statistics/inst/normplot.m trunk/octave-forge/main/statistics/inst/normstat.m trunk/octave-forge/main/statistics/inst/pdist.m trunk/octave-forge/main/statistics/inst/poisstat.m trunk/octave-forge/main/statistics/inst/princomp.m trunk/octave-forge/main/statistics/inst/private/tbl_delim.m trunk/octave-forge/main/statistics/inst/random.m trunk/octave-forge/main/statistics/inst/raylcdf.m trunk/octave-forge/main/statistics/inst/raylinv.m trunk/octave-forge/main/statistics/inst/raylpdf.m trunk/octave-forge/main/statistics/inst/raylrnd.m trunk/octave-forge/main/statistics/inst/raylstat.m trunk/octave-forge/main/statistics/inst/regress.m trunk/octave-forge/main/statistics/inst/repanova.m trunk/octave-forge/main/statistics/inst/squareform.m trunk/octave-forge/main/statistics/inst/tabulate.m trunk/octave-forge/main/statistics/inst/tblread.m trunk/octave-forge/main/statistics/inst/tblwrite.m trunk/octave-forge/main/statistics/inst/trimmean.m trunk/octave-forge/main/statistics/inst/tstat.m trunk/octave-forge/main/statistics/inst/unidstat.m trunk/octave-forge/main/statistics/inst/unifstat.m trunk/octave-forge/main/statistics/inst/vmpdf.m trunk/octave-forge/main/statistics/inst/vmrnd.m trunk/octave-forge/main/statistics/inst/wblstat.m Modified: trunk/octave-forge/main/statistics/COPYING =================================================================== --- trunk/octave-forge/main/statistics/COPYING 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/COPYING 2012-04-12 13:47:32 UTC (rev 10203) @@ -1 +1 @@ -See the individual files for licenses. +See individual files for licenses Modified: trunk/octave-forge/main/statistics/DESCRIPTION =================================================================== --- trunk/octave-forge/main/statistics/DESCRIPTION 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/DESCRIPTION 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,12 +1,12 @@ Name: Statistics Version: 1.1.1 Date: 2012-04-12 -Author: Various Authors +Author: various authors Maintainer: Arno Onken <as...@as...> Title: Statistics Description: Additional statistics functions for Octave. Categories: Statistics Depends: octave (>= 2.9.7), miscellaneous Autoload: yes -License: See individual files +License: GPLv3+, public domain Url: http://octave.sf.net Modified: trunk/octave-forge/main/statistics/inst/anderson_darling_cdf.m =================================================================== --- trunk/octave-forge/main/statistics/inst/anderson_darling_cdf.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/anderson_darling_cdf.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,3 +1,6 @@ +## Author: Paul Kienzle <pki...@us...> +## This program is granted to the public domain. + ## -*- texinfo -*- ## @deftypefn {Function File} @var{p} = anderson_darling_cdf (@var{A}, @var{n}) ## @@ -73,9 +76,6 @@ ## @seealso{anderson_darling_test} ## @end deftypefn -## Author: Paul Kienzle -## This code is granted to the public domain. - function y = anderson_darling_cdf(z,n) y = ADinf(z); y += ADerrfix(y,n); Modified: trunk/octave-forge/main/statistics/inst/anderson_darling_test.m =================================================================== --- trunk/octave-forge/main/statistics/inst/anderson_darling_test.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/anderson_darling_test.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,3 +1,6 @@ +## Author: Paul Kienzle <pki...@us...> +## This program is granted to the public domain. + ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{q}, @var{Asq}, @var{info}] = } = @ ## anderson_darling_test (@var{x}, @var{distribution}) @@ -51,9 +54,6 @@ ## @seealso{anderson_darling_cdf} ## @end deftypefn - -## Author: Paul Kienzle -## This program is granted to the public domain. function [q,Asq,info] = anderson_darling_test(x,dist) if size(x,1) == 1, x=x(:); end Modified: trunk/octave-forge/main/statistics/inst/anovan.m =================================================================== --- trunk/octave-forge/main/statistics/inst/anovan.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/anovan.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,18 +1,17 @@ -## Copyright (C) 2003-2005 Andy Adler +## Copyright (C) 2003-2005 Andy Adler <ad...@nc...> ## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2, or (at your option) -## any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This software is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this software; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{pval}, @var{f}, @var{df_b}, @var{df_e}] =} anovan (@var{data}, @var{grps}) @@ -77,7 +76,6 @@ ## g3 = {'may'; 'may'; 'may'; 'may'; 'june'; 'june'; 'june'; 'june'}; ## anovan( y', [g1',g2',g3']) ## % Fails because we always do interactions - function [PVAL, FSTAT, DF_B, DFE] = anovan (data, grps, varargin) Modified: trunk/octave-forge/main/statistics/inst/betastat.m =================================================================== --- trunk/octave-forge/main/statistics/inst/betastat.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/betastat.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{m}, @var{v}] =} betastat (@var{a}, @var{b}) Modified: trunk/octave-forge/main/statistics/inst/binostat.m =================================================================== --- trunk/octave-forge/main/statistics/inst/binostat.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/binostat.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{m}, @var{v}] =} binostat (@var{n}, @var{p}) Modified: trunk/octave-forge/main/statistics/inst/boxplot.m =================================================================== --- trunk/octave-forge/main/statistics/inst/boxplot.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/boxplot.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,20 @@ -## Copyright (C) 2002 Alberto Terruzzi +## Copyright (C) 2002 Alberto Terruzzi <t-a...@li...> +## Copyright (C) 2006 Alberto Pose <ap...@al...> +## Copyright (C) 2011 Pascal Dupuis <Pas...@wo...> +## Copyright (C) 2012 Juan Pablo Carbajal <car...@if...> ## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program; If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {@var{s} =} boxplot (@var{data}, @var{notched}, @ Modified: trunk/octave-forge/main/statistics/inst/caseread.m =================================================================== --- trunk/octave-forge/main/statistics/inst/caseread.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/caseread.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2008 Bill Denney +## Copyright (C) 2008 Bill Denney <bi...@de...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {@var{names} =} caseread (@var{filename}) Modified: trunk/octave-forge/main/statistics/inst/casewrite.m =================================================================== --- trunk/octave-forge/main/statistics/inst/casewrite.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/casewrite.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2008 Bill Denney +## Copyright (C) 2008 Bill Denney <bi...@de...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {} casewrite (@var{strmat}, @var{filename}) Modified: trunk/octave-forge/main/statistics/inst/chi2stat.m =================================================================== --- trunk/octave-forge/main/statistics/inst/chi2stat.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/chi2stat.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{m}, @var{v}] =} chi2stat (@var{n}) Modified: trunk/octave-forge/main/statistics/inst/cl_multinom.m =================================================================== --- trunk/octave-forge/main/statistics/inst/cl_multinom.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/cl_multinom.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,70 +1,70 @@ -# -*- texinfo -*- -# -# @deftypefn {Function File} {@var{CL} =} cl_multinom (@var{x}, @var{N}, @var{b}, @var{calculation_type} ) - Confidence level of multinomial portions -# Returns confidence level of multinomial parameters estimated @math{ p = x / sum(x) } with predefined confidence interval @var{b}. -# Finite population is also considered. -# -# This function calculates the level of confidence at which the samples represent the true distribution -# given that there is a predefined tolerance (confidence interval). -# This is the upside down case of the typical excercises at which we want to get the confidence interval -# given the confidence level (and the estimated parameters of the underlying distribution). -# But once we accept (lets say at elections) that we have a standard predefined -# maximal acceptable error rate (e.g. @var{b}=0.02 ) in the estimation and we just want to know that how sure we -# can be that the measured proportions are the same as in the -# entire population (ie. the expected value and mean of the samples are roghly the same) we need to use this function. -# -# @subheading Arguments -# @itemize @bullet -# @item @var{x} : int vector : sample frequencies bins -# @item @var{N} : int : Population size that was sampled by x. If N<sum(x), infinite number assumed -# @item @var{b} : real, vector : confidence interval -# if vector, it should be the size of x containing confence interval for each cells -# if scalar, each cell will have the same value of b unless it is zero or -1 -# if value is 0, b=.02 is assumed which is standard choice at elections -# otherwise it is calculated in a way that one sample in a cell alteration defines the confidence interval -# @item @var{calculation_type} : string : (Optional), described below -# "bromaghin" (default) - do not change it unless you have a good reason to do so -# "cochran" -# "agresti_cull" this is not exactly the solution at reference given below but an adjustment of the solutions above -# @end itemize -# -# @subheading Returns -# Confidence level. -# -# @subheading Example -# CL = cl_multinom( [27;43;19;11], 10000, 0.05 ) -# returns 0.69 confidence level. -# -# @subheading References -# -# "bromaghin" calculation type (default) is based on -# is based on the article -# Jeffrey F. Bromaghin, "Sample Size Determination for Interval Estimation of Multinomial Probabilities", The American Statistician vol 47, 1993, pp 203-206. -# -# "cochran" calculation type -# is based on article -# Robert T. Tortora, "A Note on Sample Size Estimation for Multinomial Populations", The American Statistician, , Vol 32. 1978, pp 100-102. -# -# "agresti_cull" calculation type -# is based on article in which Quesenberry Hurst and Goodman result is combined -# A. Agresti and B.A. Coull, "Approximate is better than \"exact\" for interval estimation of binomial portions", The American Statistician, Vol. 52, 1998, pp 119-126 -# -# @end deftypefn +## Copyright (C) 2009 Levente Torok <Tor...@gm...> +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. -# Copyright (C) 2009 Levente Torok / Tor...@gm... -# This program is free software; you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation; either version 2 of the License, or -# (at your option) any later version.^ -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License" -# along with this program; If not, see <http://www.gnu.org/licenses/>. -# +## -*- texinfo -*- +## +## @deftypefn {Function File} {@var{CL} =} cl_multinom (@var{x}, @var{N}, @var{b}, @var{calculation_type} ) - Confidence level of multinomial portions +## Returns confidence level of multinomial parameters estimated @math{ p = x / sum(x) } with predefined confidence interval @var{b}. +## Finite population is also considered. +## +## This function calculates the level of confidence at which the samples represent the true distribution +## given that there is a predefined tolerance (confidence interval). +## This is the upside down case of the typical excercises at which we want to get the confidence interval +## given the confidence level (and the estimated parameters of the underlying distribution). +## But once we accept (lets say at elections) that we have a standard predefined +## maximal acceptable error rate (e.g. @var{b}=0.02 ) in the estimation and we just want to know that how sure we +## can be that the measured proportions are the same as in the +## entire population (ie. the expected value and mean of the samples are roghly the same) we need to use this function. +## +## @subheading Arguments +## @itemize @bullet +## @item @var{x} : int vector : sample frequencies bins +## @item @var{N} : int : Population size that was sampled by x. If N<sum(x), infinite number assumed +## @item @var{b} : real, vector : confidence interval +## if vector, it should be the size of x containing confence interval for each cells +## if scalar, each cell will have the same value of b unless it is zero or -1 +## if value is 0, b=.02 is assumed which is standard choice at elections +## otherwise it is calculated in a way that one sample in a cell alteration defines the confidence interval +## @item @var{calculation_type} : string : (Optional), described below +## "bromaghin" (default) - do not change it unless you have a good reason to do so +## "cochran" +## "agresti_cull" this is not exactly the solution at reference given below but an adjustment of the solutions above +## @end itemize +## +## @subheading Returns +## Confidence level. +## +## @subheading Example +## CL = cl_multinom( [27;43;19;11], 10000, 0.05 ) +## returns 0.69 confidence level. +## +## @subheading References +## +## "bromaghin" calculation type (default) is based on +## is based on the article +## Jeffrey F. Bromaghin, "Sample Size Determination for Interval Estimation of Multinomial Probabilities", The American Statistician vol 47, 1993, pp 203-206. +## +## "cochran" calculation type +## is based on article +## Robert T. Tortora, "A Note on Sample Size Estimation for Multinomial Populations", The American Statistician, , Vol 32. 1978, pp 100-102. +## +## "agresti_cull" calculation type +## is based on article in which Quesenberry Hurst and Goodman result is combined +## A. Agresti and B.A. Coull, "Approximate is better than \"exact\" for interval estimation of binomial portions", The American Statistician, Vol. 52, 1998, pp 119-126 +## +## @end deftypefn function CL = cl_multinom( x, N, b = .05, calculation_type = "bromaghin") Modified: trunk/octave-forge/main/statistics/inst/copulacdf.m =================================================================== --- trunk/octave-forge/main/statistics/inst/copulacdf.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/copulacdf.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2008 Arno Onken +## Copyright (C) 2008 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {@var{p} =} copulacdf (@var{family}, @var{x}, @var{theta}) Modified: trunk/octave-forge/main/statistics/inst/copulapdf.m =================================================================== --- trunk/octave-forge/main/statistics/inst/copulapdf.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/copulapdf.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2008 Arno Onken +## Copyright (C) 2008 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {@var{p} =} copulapdf (@var{family}, @var{x}, @var{theta}) Modified: trunk/octave-forge/main/statistics/inst/expstat.m =================================================================== --- trunk/octave-forge/main/statistics/inst/expstat.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/expstat.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{m}, @var{v}] =} expstat (@var{l}) Modified: trunk/octave-forge/main/statistics/inst/ff2n.m =================================================================== --- trunk/octave-forge/main/statistics/inst/ff2n.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/ff2n.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,3 +1,6 @@ +## Author: Paul Kienzle <pki...@us...> +## This program is granted to the public domain. + ## -*- texinfo -*- ## @deftypefn {Function File} ff2n (@var{n}) ## Full-factor design with n binary terms. @@ -5,9 +8,6 @@ ## @seealso {fullfact} ## @end deftypefn -## Author: Paul Kienzle <pki...@us...> -## This program is public domain - function A=ff2n(n) A = fullfact (2 * ones (1,n)) - 1; endfunction Modified: trunk/octave-forge/main/statistics/inst/fstat.m =================================================================== --- trunk/octave-forge/main/statistics/inst/fstat.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/fstat.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{mn}, @var{v}] =} fstat (@var{m}, @var{n}) Modified: trunk/octave-forge/main/statistics/inst/fullfact.m =================================================================== --- trunk/octave-forge/main/statistics/inst/fullfact.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/fullfact.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,3 +1,6 @@ +## Author: Paul Kienzle <pki...@us...> +## This program is granted to the public domain. + ## -*- texinfo -*- ## @deftypefn {Function File} fullfact (@var{N}) ## Full factorial design. @@ -10,9 +13,6 @@ ## ## @end deftypefn -## Author: Paul Kienzle <pki...@us...> -## This program is public domain. - function A = fullfact(n) if length(n) == 1 % combinatorial design with n either/or choices Modified: trunk/octave-forge/main/statistics/inst/gamfit.m =================================================================== --- trunk/octave-forge/main/statistics/inst/gamfit.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/gamfit.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,3 +1,6 @@ +## Author: Martijn van Oosterhout <kl...@sv...> +## This program is granted to the public domain. + ## -*- texinfo -*- ## @deftypefn {Function File} {} [A B] = gamfit (@var{R}) ## Finds the maximumlikelihood estimator for the Gamma distribution for R @@ -15,9 +18,6 @@ ## Matlab. To work for Matlab, value of b needs to be inverted in a few ## places (marked with **) -## Written by Martijn van Oosterhout <kl...@sv...> (Nov 2006) -## This code is public domain - function res = gamfit(R) if (nargin != 1) Modified: trunk/octave-forge/main/statistics/inst/gamlike.m =================================================================== --- trunk/octave-forge/main/statistics/inst/gamlike.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/gamlike.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,3 +1,6 @@ +## Author: Martijn van Oosterhout <kl...@sv...> +## This program is granted to the public domain. + ## -*- texinfo -*- ## @deftypefn {Function File} {@var{X} =} gamlike ([@var{A} @var{B}], @var{R}) ## Calculates the negative log-likelihood function for the Gamma @@ -5,9 +8,6 @@ ## @seealso{gampdf, gaminv, gamrnd, gamfit} ## @end deftypefn -## Written by Martijn van Oosterhout <kl...@sv...> (Nov 2006) -## This code is public domain - function res = gamlike(P,K) if (nargin != 2) Modified: trunk/octave-forge/main/statistics/inst/gamstat.m =================================================================== --- trunk/octave-forge/main/statistics/inst/gamstat.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/gamstat.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{m}, @var{v}] =} gamstat (@var{a}, @var{b}) Modified: trunk/octave-forge/main/statistics/inst/geomean.m =================================================================== --- trunk/octave-forge/main/statistics/inst/geomean.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/geomean.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2001 Paul Kienzle +## Copyright (C) 2001 Paul Kienzle <pki...@us...> ## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program; If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} geomean (@var{x}) Modified: trunk/octave-forge/main/statistics/inst/geostat.m =================================================================== --- trunk/octave-forge/main/statistics/inst/geostat.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/geostat.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{m}, @var{v}] =} geostat (@var{p}) Modified: trunk/octave-forge/main/statistics/inst/harmmean.m =================================================================== --- trunk/octave-forge/main/statistics/inst/harmmean.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/harmmean.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2001 Paul Kienzle +## Copyright (C) 2001 Paul Kienzle <pki...@us...> ## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program; If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} harmmean (@var{x}) Modified: trunk/octave-forge/main/statistics/inst/histfit.m =================================================================== --- trunk/octave-forge/main/statistics/inst/histfit.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/histfit.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2003 Alberto Terruzzi +## Copyright (C) 2003 Alberto Terruzzi <t-a...@li...> ## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program; If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} histfit (@var{data}, @var{nbins}) Modified: trunk/octave-forge/main/statistics/inst/hmmestimate.m =================================================================== --- trunk/octave-forge/main/statistics/inst/hmmestimate.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/hmmestimate.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{transprobest}, @var{outprobest}] =} hmmestimate (@var{sequence}, @var{states}) Modified: trunk/octave-forge/main/statistics/inst/hmmgenerate.m =================================================================== --- trunk/octave-forge/main/statistics/inst/hmmgenerate.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/hmmgenerate.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{sequence}, @var{states}] =} hmmgenerate (@var{len}, @var{transprob}, @var{outprob}) Modified: trunk/octave-forge/main/statistics/inst/hmmviterbi.m =================================================================== --- trunk/octave-forge/main/statistics/inst/hmmviterbi.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/hmmviterbi.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {@var{vpath} =} hmmviterbi (@var{sequence}, @var{transprob}, @var{outprob}) Modified: trunk/octave-forge/main/statistics/inst/hygestat.m =================================================================== --- trunk/octave-forge/main/statistics/inst/hygestat.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/hygestat.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,17 +1,17 @@ -## Copyright (C) 2006, 2007 Arno Onken +## Copyright (C) 2006, 2007 Arno Onken <as...@as...> ## -## This program is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. ## -## You should have received a copy of the GNU General Public License -## along with this program. If not, see <http://www.gnu.org/licenses/>. +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{mn}, @var{v}] =} hygestat (@var{t}, @var{m}, @var{n}) Modified: trunk/octave-forge/main/statistics/inst/jackknife.m =================================================================== --- trunk/octave-forge/main/statistics/inst/jackknife.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/jackknife.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,18 +1,17 @@ -## Copyright (C) 2011 Alexander Klein -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. +## Copyright (C) 2011 Alexander Klein <ale...@ma...> +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn{Function File} {@var{jackstat} =} jackknife (@var{E}, @var{x}, @dots{}) Modified: trunk/octave-forge/main/statistics/inst/jsucdf.m =================================================================== --- trunk/octave-forge/main/statistics/inst/jsucdf.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/jsucdf.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,11 +1,17 @@ -## Copyright (C) 2006 Frederick (Rick) A Niles +## Copyright (C) 2006 Frederick (Rick) A Niles <ni...@ri...> ## -## This file is intended to be used with this software. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2, or (at your option) -## any later version. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {} jsucdf (@var{x}, @var{alpha1}, @var{alpha2}) Modified: trunk/octave-forge/main/statistics/inst/jsupdf.m =================================================================== --- trunk/octave-forge/main/statistics/inst/jsupdf.m 2012-04-12 10:13:25 UTC (rev 10202) +++ trunk/octave-forge/main/statistics/inst/jsupdf.m 2012-04-12 13:47:32 UTC (rev 10203) @@ -1,11 +1,17 @@ -## Copyright (C) 2006 Frederick (Rick) A Niles +## Copyright (C) 2006 Frederick (Rick) A Niles <ni...@ri...> ## -## This file is intended to be used with this software. +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. ## -## This is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2, or (at your option) -## any later version. +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {} jsupdf (@var{x}, @var{alpha1}, @var{alpha2}) Modified: trunk/octave-forge/main/statistics/inst/kmeans.m ===============================... [truncated message content] |
From: <car...@us...> - 2012-04-17 14:09:10
|
Revision: 10265 http://octave.svn.sourceforge.net/octave/?rev=10265&view=rev Author: carandraug Date: 2012-04-17 14:08:58 +0000 (Tue, 17 Apr 2012) Log Message: ----------- statistics: now dependent on io package because csv2cell and cell2csv have been moved from miscellaneous to io package Modified Paths: -------------- trunk/octave-forge/main/statistics/DESCRIPTION trunk/octave-forge/main/statistics/NEWS Modified: trunk/octave-forge/main/statistics/DESCRIPTION =================================================================== --- trunk/octave-forge/main/statistics/DESCRIPTION 2012-04-17 14:04:28 UTC (rev 10264) +++ trunk/octave-forge/main/statistics/DESCRIPTION 2012-04-17 14:08:58 UTC (rev 10265) @@ -6,7 +6,7 @@ Title: Statistics Description: Additional statistics functions for Octave. Categories: Statistics -Depends: octave (>= 2.9.7), miscellaneous +Depends: octave (>= 2.9.7), miscellaneous, io (>= 1.0.18) Autoload: yes License: GPLv3+, public domain Url: http://octave.sf.net Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-04-17 14:04:28 UTC (rev 10264) +++ trunk/octave-forge/main/statistics/NEWS 2012-04-17 14:08:58 UTC (rev 10265) @@ -1,3 +1,11 @@ +Summary of important user-visible changes for statistics 1.1.2: +------------------------------------------------------------------- + + ** The package is now dependent on the io package (version 1.0.18 or + later) since some of the functions that were previously part of + the miscellaneous package (already a dependency of statistics), + have been moved to io. + Summary of important user-visible changes for statistics 1.1.1: ------------------------------------------------------------------- This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <car...@us...> - 2012-04-24 16:06:17
|
Revision: 10322 http://octave.svn.sourceforge.net/octave/?rev=10322&view=rev Author: carandraug Date: 2012-04-24 16:06:07 +0000 (Tue, 24 Apr 2012) Log Message: ----------- kmeans: patch by Daniel Ward <dw...@gm...> * expand the returned items to closer match the Matlab kmeans * 'emptyaction' property with the 'singleton' value * handle empty cluster better * throw an error if the user does not request and empty cluster handling, and there is an empty cluster Modified Paths: -------------- trunk/octave-forge/main/statistics/NEWS trunk/octave-forge/main/statistics/inst/kmeans.m Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-04-24 08:23:36 UTC (rev 10321) +++ trunk/octave-forge/main/statistics/NEWS 2012-04-24 16:06:07 UTC (rev 10322) @@ -5,6 +5,12 @@ later) since the functions that it depended of from miscellaneous package have been moved to io. + ** The function `kmeans' now accepts the 'emptyaction' property with + the 'singleton' value. This allows for the kmeans algorithm to handle + empty cluster better. It also throws an error if the user does not + request an empty cluster handling, and there is an empty cluster. + Plus, the returned items are now a closer match to Matlab. + Summary of important user-visible changes for statistics 1.1.1: ------------------------------------------------------------------- Modified: trunk/octave-forge/main/statistics/inst/kmeans.m =================================================================== --- trunk/octave-forge/main/statistics/inst/kmeans.m 2012-04-24 08:23:36 UTC (rev 10321) +++ trunk/octave-forge/main/statistics/inst/kmeans.m 2012-04-24 16:06:07 UTC (rev 10322) @@ -1,4 +1,5 @@ ## Copyright (C) 2011 Soren Hauberg <so...@ha...> +## Copyright (C) 2012 Daniel Ward <dw...@gm...> ## ## This program is free software; you can redistribute it and/or modify it under ## the terms of the GNU General Public License as published by the Free Software @@ -13,8 +14,26 @@ ## You should have received a copy of the GNU General Public License along with ## this program; if not, see <http://www.gnu.org/licenses/>. -function [classes, centers] = kmeans (data, k, varargin) - ## Input checking +function [classes, centers, sumd, D] = kmeans (data, k, varargin) + [reg, prop] = parseparams(varargin); + + #defaults for options + emptyaction = 'error'; + start = "sample"; + + #used for getting the number of samples + [N, D] = size (data); + + #used to hold the distances from each sample to each class + D = zeros (N, k); + + #used for convergence of the centroids + err = 1; + + #initial sum of distances + sumd = Inf; + + ## Input checking, validate the matrix and k if (!ismatrix (data) || !isreal (data)) error ("kmeans: first input argument must be a DxN real data matrix"); endif @@ -22,14 +41,20 @@ error ("kmeans: second input argument must be a scalar"); endif - [N, D] = size (data); + if(length(varargin) > 0) + ## check for the 'emptyaction' property + found = find(strcmpi(prop,'emptyaction') == 1); + switch(lower(prop{found+1})) + case 'singleton' + emptyaction = 'singleton'; + otherwise + error ("kmeans: unsupported empty cluster action parameter"); + endswitch + endif - ## (so far) Harcoded options - maxiter = Inf; - start = "sample"; - - ## Find initial clusters - switch (lower (start)) + ## check for the 'start' property + #found = find(strcmpi(prop,'start') == 1); + switch (lower(start)) case "sample" idx = randperm (N) (1:k); centers = data (idx, :); @@ -37,31 +62,70 @@ error ("kmeans: unsupported initial clustering parameter"); endswitch + ## Run the algorithm - D = zeros (N, k); - iterations = 0; - prevcenters = centers; - while (true) + while err > .001 ## Compute distances for i = 1:k - D (:, i) = sum (( data - repmat (centers (i, :), N, 1)).^2, 2); + D (:, i) = sumsq(data - repmat(centers(i, :), N, 1), 2); endfor ## Classify [tmp, classes] = min (D, [], 2); - ## Recompute centers + ## Calcualte new centroids for i = 1:k + + ## Check for empty clusters + if sum(classes==i)==0 || length(mean(data(classes == i, :))) == 0 + + switch(emptyaction) + ## if 'singleton', then find the point that is the + ## farthest and add it to the empty cluster + case 'singleton' + classes(maxCostSampleIndex(data,centers(i,:))) = i; + + ## if 'error' then throw the error + otherwise + error ("kmeans: empty cluster created"); + endswitch + + endif ## end check for empty clusters + + ## update the centroids centers (i, :) = mean (data (classes == i, :)); + endfor - ## Check for convergence - iterations++; - if (all (centers (:) == prevcenters (:)) || iterations >= maxiter) - break; + ## calculate the differnece in the sum of distances + err = sumd - objCost(data,classes,centers); + ## update the current sum of distances + sumd = objCost(data,classes,centers); + endwhile +endfunction + +## calculate the sum of distances +function obj = objCost(data,classes,centers) + [N D] = size(data); + sum = 0; + + for i=1:N + sum = sum + sumsq(data(i,:)-centers(classes(i),:)); + endfor + + obj = sum; +endfunction + +function index = maxCostSampleIndex(data,centers) + [n m] = size(data); + cost = 0; + index = 1; + for j = 1:n + if cost < sumsq(data(j,:) - centers) + cost = sumsq(data(j,:) - centers); + index = j; endif - prevcenters = centers; - endwhile + endfor endfunction %!demo @@ -80,4 +144,3 @@ %! plot (data (idx==2, 1), data (idx==2, 2), 'bs'); %! plot (centers (:, 1), centers (:, 2), 'kv', 'markersize', 10); %! hold off - This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2012-04-30 21:50:19
|
Revision: 10338 http://octave.svn.sourceforge.net/octave/?rev=10338&view=rev Author: asnelt Date: 2012-04-30 21:50:12 +0000 (Mon, 30 Apr 2012) Log Message: ----------- Initial commit. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX trunk/octave-forge/main/statistics/NEWS Added Paths: ----------- trunk/octave-forge/main/statistics/inst/mnrnd.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2012-04-30 17:13:19 UTC (rev 10337) +++ trunk/octave-forge/main/statistics/INDEX 2012-04-30 21:50:12 UTC (rev 10338) @@ -4,30 +4,31 @@ betastat binostat chi2stat + cl_multinom copulacdf copulapdf expstat fstat + gamlike gamstat geostat hygestat + jsucdf + jsupdf lognstat mvnpdf mvnrnd mvncdf + mnrnd + mvtcdf + nbinstat normalise_distribution - vmpdf vmrnd - mvtcdf normstat - nbinstat poisstat + random raylcdf raylinv raylpdf raylrnd raylstat tstat unidstat unifstat + vmpdf vmrnd wblstat - gamlike - jsucdf - jsupdf - random - cl_multinom Descriptive statistics nansum nanmax Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-04-30 17:13:19 UTC (rev 10337) +++ trunk/octave-forge/main/statistics/NEWS 2012-04-30 21:50:12 UTC (rev 10338) @@ -1,6 +1,10 @@ Summary of important user-visible changes for statistics 1.1.2: ------------------------------------------------------------------- + ** The following functions are new in 1.1.2: + + mnrnd + ** The package is now dependent on the io package (version 1.0.18 or later) since the functions that it depended of from miscellaneous package have been moved to io. Added: trunk/octave-forge/main/statistics/inst/mnrnd.m =================================================================== --- trunk/octave-forge/main/statistics/inst/mnrnd.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/mnrnd.m 2012-04-30 21:50:12 UTC (rev 10338) @@ -0,0 +1,182 @@ +## Copyright (C) 2012 Arno Onken +## +## This program is free software: you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{x} =} mnrnd (@var{n}, @var{p}) +## @deftypefnx {Function File} {@var{x} =} mnrnd (@var{n}, @var{p}, @var{s}) +## Generate random samples from the multinomial distribution. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{n} is the first parameter of the multinomial distribution. @var{n} can +## be scalar or a vector containing the number of trials of each multinomial +## sample. The elements of @var{n} must be non-negative integers. +## +## @item +## @var{n} is the second parameter of the multinomial distribution. @var{p} can +## be a vector with the probabilities of the categories or a matrix with each +## row containing the probabilities of a multinomial sample. If @var{p} has +## more than one row and @var{n} is non-scalar, then the number of rows of +## @var{p} must match the number of elements of @var{n}. +## +## @item +## @var{s} is the number of multinomial samples to be generated. @var{s} must +## be a non-negative integer. If @var{s} is specified, then @var{n} must be +## scalar and @var{p} must be a vector. +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{x} is a matrix of random samples from the multinomial distribution with +## corresponding parameters @var{n} and @var{p}. Each row corresponds to one +## multinomial sample. The number of columns, therefore, corresponds to the +## number of columns of @var{p}. If @var{s} is not specified, then the number +## of rows of @var{x} is the maximum of the number of elements of @var{n} and +## the number of rows of @var{p}. If a row of @var{p} does not sum to @code{1}, +## then the corresponding row of @var{x} will contain only @code{NaN} values. +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## n = 10; +## p = [0.2 0.5 0.3]; +## x = mnrnd (n, p); +## @end group +## +## @group +## n = 10 * ones (3, 1); +## p = [0.2 0.5 0.3]; +## x = mnrnd (n, p); +## @end group +## +## @group +## n = (1:2)'; +## p = [0.2 0.5 0.3; 0.1 0.1 0.8]; +## x = mnrnd (n, p); +## @end group +## @end example +## +## @subheading References +## +## @enumerate +## @item +## Wendy L. Martinez and Angel R. Martinez. @cite{Computational Statistics +## Handbook with MATLAB}. Appendix E, pages 547-557, Chapman & Hall/CRC, 2001. +## +## @item +## Athanasios Papoulis. @cite{Probability, Random Variables, and Stochastic +## Processes}. pages 104 and 148, McGraw-Hill, New York, second edition, 1984. +## @end enumerate +## @end deftypefn + +## Author: Arno Onken <as...@as...> +## Description: Random samples from the multinomial distribution + +function x = mnrnd (n, p, s) + + # Check arguments + if (nargin == 3) + if (! isscalar (n) || n < 0 || round (n) != n) + error ("mnrnd: n must be a non-negative integer"); + endif + if (! isvector (p) || any (p < 0 | p > 1)) + error ("mnrnd: p must be a vector of probabilities"); + endif + if (! isscalar (s) || s < 0 || round (s) != s) + error ("mnrnd: s must be a non-negative integer"); + endif + elseif (nargin == 2) + if (isvector (p) && size (p, 1) > 1) + p = p'; + endif + if (! isvector (n) || any (n < 0 | round (n) != n) || size (n, 2) > 1) + error ("mnrnd: n must be a non-negative integer column vector"); + endif + if (! ismatrix (p) || isempty (p) || any (p < 0 | p > 1)) + error ("mnrnd: p must be a non-empty matrix with rows of probabilities"); + endif + if (! isscalar (n) && size (p, 1) > 1 && length (n) != size (p, 1)) + error ("mnrnd: the length of n must match the number of rows of p"); + endif + else + print_usage (); + endif + + # Adjust input sizes + if (nargin == 3) + n = n * ones (s, 1); + p = repmat (p(:)', s, 1); + elseif (nargin == 2) + if (isscalar (n) && size (p, 1) > 1) + n = n * ones (size (p, 1), 1); + elseif (size (p, 1) == 1) + p = repmat (p, length (n), 1); + endif + endif + sz = size (p); + + # Upper bounds of categories + ub = cumsum (p, 2); + # Lower bounds of categories + lb = [zeros(sz(1), 1) ub(:, 1:(end-1))]; + + # Draw multinomial samples + x = zeros (sz); + for i = 1:sz(1) + # Check if the probabilities sum to 1 + if (ub(i, end) == 1) + # Draw uniform random numbers + r = repmat (rand (n(i), 1), 1, sz(2)); + # Compare the random numbers of r to the cumulated probabilities of p and + # count the number of samples for each category + x(i, :) = sum (r <= repmat (ub(i, :), n(i), 1) & r > repmat (lb(i, :), n(i), 1), 1); + else + x(i, :) = NaN; + endif + end + +endfunction + +%!test +%! n = 10; +%! p = [0.2 0.5 0.3]; +%! x = mnrnd (n, p); +%! assert (size (x), size (p)); +%! assert (all (x >= 0)); +%! assert (all (round (x) == x)); +%! assert (sum (x) == n); + +%!test +%! n = 10 * ones (3, 1); +%! p = [0.2 0.5 0.3]; +%! x = mnrnd (n, p); +%! assert (size (x), [length(n) length(p)]); +%! assert (all (x >= 0)); +%! assert (all (round (x) == x)); +%! assert (all (sum (x, 2) == n)); + +%!test +%! n = (1:2)'; +%! p = [0.2 0.5 0.3; 0.1 0.1 0.8]; +%! x = mnrnd (n, p); +%! assert (size (x), size (p)); +%! assert (all (x >= 0)); +%! assert (all (round (x) == x)); +%! assert (all (sum (x, 2) == n)); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2012-05-01 16:20:49
|
Revision: 10346 http://octave.svn.sourceforge.net/octave/?rev=10346&view=rev Author: asnelt Date: 2012-05-01 16:20:43 +0000 (Tue, 01 May 2012) Log Message: ----------- Initial commit. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX trunk/octave-forge/main/statistics/NEWS Added Paths: ----------- trunk/octave-forge/main/statistics/inst/mnpdf.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2012-05-01 16:18:26 UTC (rev 10345) +++ trunk/octave-forge/main/statistics/INDEX 2012-05-01 16:20:43 UTC (rev 10346) @@ -16,7 +16,7 @@ jsupdf lognstat mvnpdf mvnrnd mvncdf - mnrnd + mnpdf mnrnd mvtcdf nbinstat normalise_distribution Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-05-01 16:18:26 UTC (rev 10345) +++ trunk/octave-forge/main/statistics/NEWS 2012-05-01 16:20:43 UTC (rev 10346) @@ -3,7 +3,7 @@ ** The following functions are new in 1.1.2: - mnrnd + mnpdf mnrnd ** The package is now dependent on the io package (version 1.0.18 or later) since the functions that it depended of from miscellaneous Added: trunk/octave-forge/main/statistics/inst/mnpdf.m =================================================================== --- trunk/octave-forge/main/statistics/inst/mnpdf.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/mnpdf.m 2012-05-01 16:20:43 UTC (rev 10346) @@ -0,0 +1,132 @@ +## Copyright (C) 2012 Arno Onken +## +## This program is free software: you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{y} =} mnpdf (@var{x}, @var{p}) +## Compute the probability density function of the multinomial distribution. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{x} is vector with a single sample of a multinomial distribution with +## parameter @var{p} or a matrix of random samples from multinomial +## distributions. In the latter case, each row of @var{x} is a sample from a +## multinomial distribution with the corresponding row of @var{p} being its +## parameter. +## +## @item +## @var{p} is a vector with the probabilities of the categories or a matrix +## with each row containing the probabilities of a multinomial sample. +## @end itemize +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{y} is a vector of probabilites of the random samples @var{x} from the +## multinomial distribution with corresponding parameter @var{p}. The parameter +## @var{n} of the multinomial distribution is the sum of the elements of each +## row of @var{x}. The length of @var{y} is the number of columns of @var{x}. +## If a row of @var{p} does not sum to @code{1}, then the corresponding element +## of @var{y} will be @code{NaN}. +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## x = [1, 4, 2]; +## p = [0.2, 0.5, 0.3]; +## y = mnpdf (x, p); +## @end group +## +## @group +## x = [1, 4, 2; 1, 0, 9]; +## p = [0.2, 0.5, 0.3; 0.1, 0.1, 0.8]; +## y = mnpdf (x, p); +## @end group +## @end example +## +## @subheading References +## +## @enumerate +## @item +## Wendy L. Martinez and Angel R. Martinez. @cite{Computational Statistics +## Handbook with MATLAB}. Appendix E, pages 547-557, Chapman & Hall/CRC, 2001. +## +## @item +## Merran Evans, Nicholas Hastings and Brian Peacock. @cite{Statistical +## Distributions}. pages 134-136, Wiley, New York, third edition, 2000. +## @end enumerate +## @end deftypefn + +## Author: Arno Onken <as...@as...> +## Description: PDF of the multinomial distribution + +function y = mnpdf (x, p) + + # Check arguments + if (nargin != 2) + print_usage (); + endif + + if (! ismatrix (x) || any (x(:) < 0 | round (x(:) != x(:)))) + error ("mnpdf: x must be a matrix of non-negative integer values"); + endif + if (! ismatrix (p) || any (p(:) < 0)) + error ("mnpdf: p must be a non-empty matrix with rows of probabilities"); + endif + + # Adjust input sizes + if (! isvector (x) || ! isvector (p)) + if (isvector (x)) + x = x(:)'; + endif + if (isvector (p)) + p = p(:)'; + endif + if (size (x, 1) == 1 && size (p, 1) > 1) + x = repmat (x, size (p, 1), 1); + elseif (size (x, 1) > 1 && size (p, 1) == 1) + p = repmat (p, size (x, 1), 1); + endif + endif + # Continue argument check + if (any (size (x) != size (p))) + error ("mnpdf: x and p must have compatible sizes"); + endif + + # Count total number of elements of each multinomial sample + n = sum (x, 2); + # Compute probability density function of the multinomial distribution + t = x .* log (p); + t(x == 0) = 0; + y = exp (gammaln (n+1) - sum (gammaln (x+1), 2) + sum (t, 2)); + y(sum (p, 2) != 1) = NaN; + +endfunction + +%!test +%! x = [1, 4, 2]; +%! p = [0.2, 0.5, 0.3]; +%! y = mnpdf (x, p); +%! assert (y, 0.11812, 0.001); + +%!test +%! x = [1, 4, 2; 1, 0, 9]; +%! p = [0.2, 0.5, 0.3; 0.1, 0.1, 0.8]; +%! y = mnpdf (x, p); +%! assert (y, [0.11812; 0.13422], 0.001); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2012-05-06 14:33:01
|
Revision: 10366 http://octave.svn.sourceforge.net/octave/?rev=10366&view=rev Author: asnelt Date: 2012-05-06 14:32:55 +0000 (Sun, 06 May 2012) Log Message: ----------- Relaxed condition for invalid rows in mnpdf and mnrnd so that larger numbers of categories are supported. Modified Paths: -------------- trunk/octave-forge/main/statistics/NEWS trunk/octave-forge/main/statistics/inst/mnpdf.m trunk/octave-forge/main/statistics/inst/mnrnd.m Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-05-05 13:03:58 UTC (rev 10365) +++ trunk/octave-forge/main/statistics/NEWS 2012-05-06 14:32:55 UTC (rev 10366) @@ -1,3 +1,9 @@ +Summary of important user-visible changes for statistics 1.1.3: +------------------------------------------------------------------- + + ** The functions mnpdf and mnrnd are now also usable for greater numbers + of categories for which the rows do not exactly sum to 1. + Summary of important user-visible changes for statistics 1.1.2: ------------------------------------------------------------------- Modified: trunk/octave-forge/main/statistics/inst/mnpdf.m =================================================================== --- trunk/octave-forge/main/statistics/inst/mnpdf.m 2012-05-05 13:03:58 UTC (rev 10365) +++ trunk/octave-forge/main/statistics/inst/mnpdf.m 2012-05-06 14:32:55 UTC (rev 10366) @@ -115,7 +115,9 @@ t = x .* log (p); t(x == 0) = 0; y = exp (gammaln (n+1) - sum (gammaln (x+1), 2) + sum (t, 2)); - y(sum (p, 2) != 1) = NaN; + # Set invalid rows to NaN + k = (abs (sum (p, 2) - 1) > 1e-6); + y(k) = NaN; endfunction Modified: trunk/octave-forge/main/statistics/inst/mnrnd.m =================================================================== --- trunk/octave-forge/main/statistics/inst/mnrnd.m 2012-05-05 13:03:58 UTC (rev 10365) +++ trunk/octave-forge/main/statistics/inst/mnrnd.m 2012-05-06 14:32:55 UTC (rev 10366) @@ -135,23 +135,24 @@ # Upper bounds of categories ub = cumsum (p, 2); + # Make sure that the greatest upper bound is 1 + gub = ub(:, end); + ub(:, end) = 1; # Lower bounds of categories lb = [zeros(sz(1), 1) ub(:, 1:(end-1))]; # Draw multinomial samples x = zeros (sz); for i = 1:sz(1) - # Check if the probabilities sum to 1 - if (ub(i, end) == 1) - # Draw uniform random numbers - r = repmat (rand (n(i), 1), 1, sz(2)); - # Compare the random numbers of r to the cumulated probabilities of p and - # count the number of samples for each category - x(i, :) = sum (r <= repmat (ub(i, :), n(i), 1) & r > repmat (lb(i, :), n(i), 1), 1); - else - x(i, :) = NaN; - endif + # Draw uniform random numbers + r = repmat (rand (n(i), 1), 1, sz(2)); + # Compare the random numbers of r to the cumulated probabilities of p and + # count the number of samples for each category + x(i, :) = sum (r <= repmat (ub(i, :), n(i), 1) & r > repmat (lb(i, :), n(i), 1), 1); endfor + # Set invalid rows to NaN + k = (abs (gub - 1) > 1e-6); + x(k, :) = NaN; endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2012-05-09 19:56:26
|
Revision: 10392 http://octave.svn.sourceforge.net/octave/?rev=10392&view=rev Author: asnelt Date: 2012-05-09 19:56:20 +0000 (Wed, 09 May 2012) Log Message: ----------- Update for release 1.1.3. Modified Paths: -------------- trunk/octave-forge/main/statistics/DESCRIPTION trunk/octave-forge/main/statistics/NEWS Modified: trunk/octave-forge/main/statistics/DESCRIPTION =================================================================== --- trunk/octave-forge/main/statistics/DESCRIPTION 2012-05-09 19:54:04 UTC (rev 10391) +++ trunk/octave-forge/main/statistics/DESCRIPTION 2012-05-09 19:56:20 UTC (rev 10392) @@ -1,6 +1,6 @@ Name: Statistics -Version: 1.1.2 -Date: 2012-05-01 +Version: 1.1.3 +Date: 2012-05-09 Author: various authors Maintainer: Arno Onken <as...@as...> Title: Statistics Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-05-09 19:54:04 UTC (rev 10391) +++ trunk/octave-forge/main/statistics/NEWS 2012-05-09 19:56:20 UTC (rev 10392) @@ -1,6 +1,10 @@ Summary of important user-visible changes for statistics 1.1.3: ------------------------------------------------------------------- + ** The following functions are new in 1.1.3: + + copularnd mvtrnd + ** The functions mnpdf and mnrnd are now also usable for greater numbers of categories for which the rows do not exactly sum to 1. This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <jpi...@us...> - 2012-06-06 17:16:02
|
Revision: 10571 http://octave.svn.sourceforge.net/octave/?rev=10571&view=rev Author: jpicarbajal Date: 2012-06-06 17:15:51 +0000 (Wed, 06 Jun 2012) Log Message: ----------- statistics: adding regression based on gaussian processes. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX trunk/octave-forge/main/statistics/NEWS Added Paths: ----------- trunk/octave-forge/main/statistics/inst/regress_gp.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2012-06-05 14:20:07 UTC (rev 10570) +++ trunk/octave-forge/main/statistics/INDEX 2012-06-06 17:15:51 UTC (rev 10571) @@ -51,6 +51,7 @@ monotone_smooth princomp regress + regress_gp Plots boxplot normplot Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-06-05 14:20:07 UTC (rev 10570) +++ trunk/octave-forge/main/statistics/NEWS 2012-06-06 17:15:51 UTC (rev 10571) @@ -1,3 +1,10 @@ +Summary of important user-visible changes for statistics 1.2.0: +------------------------------------------------------------------- + + ** The following functions are new in 1.2.0: + + regress_gp + Summary of important user-visible changes for statistics 1.1.3: ------------------------------------------------------------------- Added: trunk/octave-forge/main/statistics/inst/regress_gp.m =================================================================== --- trunk/octave-forge/main/statistics/inst/regress_gp.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/regress_gp.m 2012-06-06 17:15:51 UTC (rev 10571) @@ -0,0 +1,123 @@ +## Copyright (c) 2012 Juan Pablo Carbajal <car...@if...> +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{m}, @var{K}] =} regress_gp (@var{x}, @var{y}, @var{Sp}) +## @deftypefnx {Function File} {[@dots{} @var{yi} @var{dy}] =} sqp (@dots{}, @var{xi}) +## Linear scalar regression using gaussian processes. +## +## It estimates the model @var{y} = @var{x}'*m for @var{x} R^D and @var{y} in R. +## The information about errors of the predictions (interpolation/extrapolation) is given +## by the covarianve matrix @var{K}. If D==1 the inputs must be column vectors, +## if D>1 then @var{x} is n-by-D, with n the number of data points. @var{Sp} defines +## the prior covariance of @var{m}, it should be a (D+1)-by-(D+1) positive definite matrix, +## if it is empty, the default is @code{Sp = 100*eye(size(x,2)+1)}. +## +## If @var{xi} inputs are provided, the model is evaluated and returned in @var{yi}. +## The estimation of the variation of @var{yi} are given in @var{dy}. +## +## Run @code{demo regress_gp} to see an examples. +## +## The function is a direc implementation of the formulae in pages 11-12 of +## Gaussian Processes for Machine Learning. Carl Edward Rasmussen and @ +## Christopher K. I. Williams. The MIT Press, 2006. ISBN 0-262-18253-X. +## available online at @url{http://gaussianprocess.org/gpml/}. +## +## @seealso{regress} +## @end deftypefn + +function [wm K yi dy] = regress_gp (x,y,Sp=[],xi=[]) + + if isempty(Sp) + Sp = 100*eye(size(x,2)+1); + end + + x = [ones(1,size(x,1)); x']; + + A = x*x' + inv (Sp); + K = inv (A); + wm = K*x*y; + + yi =[]; + dy =[]; + if !isempty (xi); + xi = [ones(size(xi,1),1) xi]; + yi = xi*wm; + dy = diag (xi*K*xi'); + end + +endfunction + +%!demo +%! % 1D Data +%! x = 2*rand (5,1)-1; +%! y = 2*x -1 + 0.3*randn (5,1); +%! +%! % Points for interpolation/extrapolation +%! xi = linspace (-2,2,10)'; +%! +%! [m K yi dy] = regress_gp (x,y,[],xi); +%! +%! plot (x,y,'xk',xi,yi,'r-',xi,yi+[-dy +dy],'b-'); + +%!demo +%! % 2D Data +%! x = 2*rand (4,2)-1; +%! y = 2*x(:,1)-3*x(:,2) -1 + 1*randn (4,1); +%! +%! % Mesh for interpolation/extrapolation +%! [xi yi] = meshgrid (linspace (-1,1,10)); +%! +%! [m K zi dz] = regress_gp (x,y,[],[xi(:) yi(:)]); +%! zi = reshape (zi, 10,10); +%! dz = reshape (dz,10,10); +%! +%! plot3 (x(:,1),x(:,2),y,'.g','markersize',8); +%! hold on; +%! h = mesh (xi,yi,zi,zeros(10,10)); +%! set(h,'facecolor','none'); +%! h = mesh (xi,yi,zi+dz,ones(10,10)); +%! set(h,'facecolor','none'); +%! h = mesh (xi,yi,zi-dz,ones(10,10)); +%! set(h,'facecolor','none'); +%! hold off +%! axis tight +%! view(80,25) + +%!demo +%! % Projection over basis function +%! pp = [2 2 0.3 1]; +%! n = 10; +%! x = 2*rand (n,1)-1; +%! y = polyval(pp,x) + 0.3*randn (n,1); +%! +%! % Powers +%! px = [sqrt(abs(x)) x x.^2 x.^3]; +%! +%! % Points for interpolation/extrapolation +%! xi = linspace (-1,1,100)'; +%! pxi = [sqrt(abs(xi)) xi xi.^2 xi.^3]; +%! +%! Sp = 100*eye(size(px,2)+1); +%! Sp(2,2) = 1; # We don't believe the sqrt is present +%! [m K yi dy] = regress_gp (px,y,Sp,pxi); +%! disp(m) +%! +%! plot (x,y,'xk;Data;',xi,yi,'r-;Estimation;',xi,polyval(pp,xi),'g-;True;'); +%! axis tight +%! axis manual +%! hold on +%! plot (xi,yi+[-dy +dy],'b-'); +%! hold off This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <jpi...@us...> - 2012-06-07 08:26:52
|
Revision: 10575 http://octave.svn.sourceforge.net/octave/?rev=10575&view=rev Author: jpicarbajal Date: 2012-06-07 08:26:41 +0000 (Thu, 07 Jun 2012) Log Message: ----------- statistics: mvnrnd: improving doc. Adding tolerance in check for positive definite. Using broadcasting. Using chol 'upper'. Perturnbing matrix to estabilize chol. Argument tol has a defualt value. Modified Paths: -------------- trunk/octave-forge/main/statistics/NEWS trunk/octave-forge/main/statistics/inst/mvnrnd.m Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-06-06 21:22:02 UTC (rev 10574) +++ trunk/octave-forge/main/statistics/NEWS 2012-06-07 08:26:41 UTC (rev 10575) @@ -1,10 +1,14 @@ Summary of important user-visible changes for statistics 1.2.0: ------------------------------------------------------------------- - ** The following functions are new in 1.2.0: + ** The following functions are new: regress_gp + + ** The interface of the following functions has been modified: + mvnrnd + Summary of important user-visible changes for statistics 1.1.3: ------------------------------------------------------------------- Modified: trunk/octave-forge/main/statistics/inst/mvnrnd.m =================================================================== --- trunk/octave-forge/main/statistics/inst/mvnrnd.m 2012-06-06 21:22:02 UTC (rev 10574) +++ trunk/octave-forge/main/statistics/inst/mvnrnd.m 2012-06-07 08:26:41 UTC (rev 10575) @@ -16,12 +16,20 @@ ## -*- texinfo -*- ## @deftypefn {Function File} @var{s} = mvnrnd (@var{mu}, @var{Sigma}) ## @deftypefnx{Function File} @var{s} = mvnrnd (@var{mu}, @var{Sigma}, @var{n}) +## @deftypefnx{Function File} @var{s} = mvnrnd (@dots{}, @var{tol}) ## Draw @var{n} random @var{d}-dimensional vectors from a multivariate Gaussian ## distribution with mean @var{mu}(@var{n}x@var{d}) and covariance matrix ## @var{Sigma}(@var{d}x@var{d}). +## +## If @var{n} is given then @var{mu} can be a 1-by-@var{d} vector. If it is not +## given @var{mu} must be @var{n}-by-@var{d}. +## +## If the argument @var{tol} is given the eigenvalues of @var{Sigma} are checked +## for positivity against -100*tol. the default value of tol is @code{eps*norm (Sigma, "fro")}. +## ## @end deftypefn -function s = mvnrnd(mu,Sigma,K) +function s = mvnrnd_jpi (mu, Sigma, K, tol=eps*norm (Sigma, "fro")) % Iain Murray 2003 -- I got sick of this simple thing not being in Octave and % locking up a stats-toolbox license in Matlab for no good @@ -32,36 +40,56 @@ % * Add GPL notice. % * Add docs for argument K + % 2012 Juan Pablo Carbajal <car...@if...> + % * Uses Octave 3.6.2 boradcast. + % * Stabilizes chol by perturbing Sigma with a epsilon multiple of the identity. + % The effect on the generated samples is to add additional independent noise + % of variance epsilon. Ref: GPML Rasmussen & Williams. 2006. pp 200-201 + % * Improved doc. + % * Added tolerance to the positive definite check + % * Used chol with option 'upper'. + + warning ("off", "Octave:broadcast"); + % If mu is column vector and Sigma not a scalar then assume user didn't read % help but let them off and flip mu. Don't be more liberal than this or it will % encourage errors (eg what should you do if mu is square?). - if ((size(mu,2)==1)&&(size(Sigma)~=[1,1])) - mu=mu'; + if (size (mu, 2) == 1) && (size (Sigma) != [1,1]) + mu = mu'; end - if nargin==3 - mu=repmat(mu,K,1); + [n,d] = size (mu); + + if nargin >= 3 + n = K; end - [n,d]=size(mu); - - if (size(Sigma)~=[d,d]) - error('Sigma must have dimensions dxd where mu is nxd.'); + if size (Sigma) != [d,d] + error ('Sigma must have dimensions dxd where mu is nxd.'); end try - U=chol(Sigma); + U = chol (Sigma + tol*eye (d),"upper"); catch - [E,Lambda]=eig(Sigma); - if (min(diag(Lambda))<0),error('Sigma must be positive semi-definite.'),end - U = sqrt(Lambda)*E'; + [E , Lambda] = eig (Sigma); + + if min (diag (Lambda)) < -100*tol + error('Sigma must be positive semi-definite. Lowest eigenvaluen %g', ... + min (diag (Lambda))); + else + Lambda(Lambda<0) = 0; + end + warning ("mvnrnd:InvalidInput","Cholesky factorization failed. Using diagonalized matrix.") + U = sqrt (Lambda) * E'; end s = randn(n,d)*U + mu; + + warning ("on", "Octave:broadcast"); endfunction % {{{ END OF CODE --- Guess I should provide an explanation: -% +% % We can draw from axis aligned unit Gaussians with randn(d) % x ~ A*exp(-0.5*x'*x) % We can then rotate this distribution using @@ -76,7 +104,7 @@ % Sigma = U'*U % For a given Sigma we can use the chol function to find the corresponding U, % draw x and find y. We can adjust for a non-zero mean by just adding it on. -% +% % But the Cholsky decomposition function doesn't always work... % Consider Sigma=[1 1;1 1]. Now inv(Sigma) doesn't actually exist, but Matlab's % mvnrnd provides samples with this covariance st x(1)~N(0,1) x(2)=x(1). The @@ -94,7 +122,7 @@ % so we can give up. % % Paul Kienzle adds: -% Where it exists, chol(Sigma) is numerically well behaved. chol(hilb(12)) +% Where it exists, chol(Sigma) is numerically well behaved. chol(hilb(12)) % for doubles and for 100 digit floating point differ in the last digit. % Where chol(Sigma) doesn't exist, X*sqrt(Lambda)*E' will be somewhat % accurate. For example, the elements of sqrt(Lambda)*E' for hilb(12), @@ -102,4 +130,3 @@ % tested using the TNT+JAMA for eig and chol templates, and qlib for % 100 digit precision. % }}} - This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <car...@us...> - 2012-07-19 03:22:06
|
Revision: 10754 http://octave.svn.sourceforge.net/octave/?rev=10754&view=rev Author: carandraug Date: 2012-07-19 03:22:00 +0000 (Thu, 19 Jul 2012) Log Message: ----------- kmeans: bugfix to deal with clusters that have only one element (patch by Christoph Kling <ck...@un...>) Modified Paths: -------------- trunk/octave-forge/main/statistics/NEWS trunk/octave-forge/main/statistics/inst/kmeans.m Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-07-18 16:16:30 UTC (rev 10753) +++ trunk/octave-forge/main/statistics/NEWS 2012-07-19 03:22:00 UTC (rev 10754) @@ -9,6 +9,9 @@ mvnrnd + ** `kmeans' has been fixed to deal with clusters that contain only + one element. + Summary of important user-visible changes for statistics 1.1.3: ------------------------------------------------------------------- Modified: trunk/octave-forge/main/statistics/inst/kmeans.m =================================================================== --- trunk/octave-forge/main/statistics/inst/kmeans.m 2012-07-18 16:16:30 UTC (rev 10753) +++ trunk/octave-forge/main/statistics/inst/kmeans.m 2012-07-19 03:22:00 UTC (rev 10754) @@ -33,20 +33,20 @@ ## used to hold the distances from each sample to each class D = zeros (nRows, k); - + #used for convergence of the centroids err = 1; - + #initial sum of distances sumd = Inf; - + ## Input checking, validate the matrix and k if (!isnumeric (data) || !ismatrix (data) || !isreal (data)) error ("kmeans: first input argument must be a DxN real data matrix"); elseif (!isscalar (k)) error ("kmeans: second input argument must be a scalar"); endif - + if (length (varargin) > 0) ## check for the 'emptyaction' property found = find (strcmpi (prop, "emptyaction") == 1); @@ -57,7 +57,7 @@ error ("kmeans: unsupported empty cluster action parameter"); endswitch endif - + ## check for the 'start' property switch (lower (start)) case "sample" @@ -66,27 +66,30 @@ otherwise error ("kmeans: unsupported initial clustering parameter"); endswitch - + ## Run the algorithm while err > .001 ## Compute distances for i = 1:k D (:, i) = sumsq (data - repmat (centers(i, :), nRows, 1), 2); endfor - + ## Classify [tmp, classes] = min (D, [], 2); - - ## Calcualte new centroids + + ## Calculate new centroids for i = 1:k + ## Get binary vector indicating membership in cluster i + membership = (classes == i); ## Check for empty clusters - if (sum (classes == i) ==0 || length (mean (data(classes == i, :))) == 0) - + if (sum (membership) == 0) switch emptyaction ## if 'singleton', then find the point that is the ## farthest and add it to the empty cluster case 'singleton' - classes(maxCostSampleIndex (data, centers(i,:))) = i; + idx=maxCostSampleIndex (data, centers(i,:)); + classes(idx) = i; + membership(idx)=1; ## if 'error' then throw the error otherwise error ("kmeans: empty cluster created"); @@ -94,10 +97,11 @@ endif ## end check for empty clusters ## update the centroids - centers(i, :) = mean (data(classes == i, :)); + members = data(membership, :); + centers(i, :) = sum(members,1)/size(members,1); endfor - ## calculate the differnece in the sum of distances + ## calculate the difference in the sum of distances err = sumd - objCost (data, classes, centers); ## update the current sum of distances sumd = objCost (data, classes, centers); @@ -112,11 +116,11 @@ endfor endfunction -function index = maxCostSampleIndex (data, centers) +function idx = maxCostSampleIndex (data, centers) cost = 0; - for index = 1:rows (data) - if cost < sumsq (data(index,:) - centers) - cost = sumsq (data(index,:) - centers); + for idx = 1:rows (data) + if cost < sumsq (data(idx,:) - centers) + cost = sumsq (data(idx,:) - centers); endif endfor endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |