From 11dfcf441e03e762efb448ee9d613a04005d7246 Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Sun, 30 Jun 2024 22:26:39 +0300 Subject: [PATCH 01/21] Add Octave to languages.yml --- lib/linguist/languages.yml | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/lib/linguist/languages.yml b/lib/linguist/languages.yml index aae547bd83..962d75754c 100644 --- a/lib/linguist/languages.yml +++ b/lib/linguist/languages.yml @@ -4020,15 +4020,11 @@ M4Sugar: MATLAB: type: programming color: "#e16737" - aliases: - - octave extensions: - ".matlab" - ".m" tm_scope: source.matlab ace_mode: matlab - codemirror_mode: octave - codemirror_mime_type: text/x-octave language_id: 225 MAXScript: type: programming @@ -4884,6 +4880,14 @@ Objective-J: tm_scope: source.js.objj ace_mode: text language_id: 259 +Octave: + type: programming + color: "#2299c4" + extensions: + - ".m" + tm_scope: source.octave + ace_mode: octave + language_id: 312 Odin: type: programming color: "#60AFFE" From d1be7b0babd07215182cf0159c84ff0c8137551b Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Sun, 30 Jun 2024 22:30:21 +0300 Subject: [PATCH 02/21] Update heuristics.yml to distinguish Octave from MATLAB --- lib/linguist/heuristics.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/linguist/heuristics.yml b/lib/linguist/heuristics.yml index 8a6fcea078..5cfa8197eb 100644 --- a/lib/linguist/heuristics.yml +++ b/lib/linguist/heuristics.yml @@ -417,6 +417,8 @@ disambiguations: rules: - language: Objective-C named_pattern: objectivec + - language: Octave + pattern: '^\s*##' - language: Mercury pattern: ':- module' - language: MUF From e64f599d4ddc85d990a3d91b4aae81a9686a9c8d Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Tue, 9 Jul 2024 16:21:00 +0300 Subject: [PATCH 03/21] Add sample file for Octave language --- samples/Octave/accumarray.m | 466 ++++++++++++++++++++++++++++++++++++ 1 file changed, 466 insertions(+) create mode 100644 samples/Octave/accumarray.m diff --git a/samples/Octave/accumarray.m b/samples/Octave/accumarray.m new file mode 100644 index 0000000000..4f850c12ea --- /dev/null +++ b/samples/Octave/accumarray.m @@ -0,0 +1,466 @@ +######################################################################## +## +## Copyright (C) 2007-2024 The Octave Project Developers +## +## See the file COPYRIGHT.md in the top-level directory of this +## distribution or . +## +## This file is part of Octave. +## +## 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 3 of the License, or +## (at your option) any later version. +## +## Octave 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 +## . +## +######################################################################## + +## -*- texinfo -*- +## @deftypefn {} {@var{A} =} accumarray (@var{subs}, @var{vals}) +## @deftypefnx {} {@var{A} =} accumarray (@var{subs}, @var{vals}, @var{sz}) +## @deftypefnx {} {@var{A} =} accumarray (@var{subs}, @var{vals}, @var{sz}, @var{fcn}) +## @deftypefnx {} {@var{A} =} accumarray (@var{subs}, @var{vals}, @var{sz}, @var{fcn}, @var{fillval}) +## @deftypefnx {} {@var{A} =} accumarray (@var{subs}, @var{vals}, @var{sz}, @var{fcn}, @var{fillval}, @var{issparse}) +## +## Create an array by accumulating the elements of a vector into the +## positions defined by their subscripts. +## +## The subscripts are defined by the rows of the matrix @var{subs} and the +## values by @var{vals}. Each row of @var{subs} corresponds to one of the +## values in @var{vals}. If @var{vals} is a scalar, it will be used for each +## of the row of @var{subs}. If @var{subs} is a cell array of vectors, all +## vectors must be of the same length, and the subscripts in the @var{k}th +## vector must correspond to the @var{k}th dimension of the result. +## +## The size of the matrix will be determined by the subscripts +## themselves. However, if @var{sz} is defined it determines the matrix +## size. The length of @var{sz} must correspond to the number of columns +## in @var{subs}. An exception is if @var{subs} has only one column, in +## which case @var{sz} may be the dimensions of a vector and the +## subscripts of @var{subs} are taken as the indices into it. +## +## The default action of @code{accumarray} is to sum the elements with +## the same subscripts. This behavior can be modified by defining the +## @var{fcn} function. This should be a function or function handle +## that accepts a column vector and returns a scalar. The result of the +## function should not depend on the order of the subscripts. +## +## The elements of the returned array that have no subscripts associated +## with them are set to zero. Defining @var{fillval} to some other value +## allows these values to be defined. This behavior changes, however, +## for certain values of @var{fcn}. If @var{fcn} is @code{@@min} +## (respectively, @code{@@max}) then the result will be filled with the +## minimum (respectively, maximum) integer if @var{vals} is of integral +## type, logical false (respectively, logical true) if @var{vals} is of +## logical type, zero if @var{fillval} is zero and all values are +## non-positive (respectively, non-negative), and NaN otherwise. +## +## By default @code{accumarray} returns a full matrix. If +## @var{issparse} is logically true, then a sparse matrix is returned +## instead. +## +## The following @code{accumarray} example constructs a frequency table +## that in the first column counts how many occurrences each number in +## the second column has, taken from the vector @var{x}. Note the usage +## of @code{unique} for assigning to all repeated elements of @var{x} +## the same index (@pxref{XREFunique,,@code{unique}}). +## +## @example +## @group +## @var{x} = [91, 92, 90, 92, 90, 89, 91, 89, 90, 100, 100, 100]; +## [@var{u}, ~, @var{j}] = unique (@var{x}); +## [accumarray(@var{j}', 1), @var{u}'] +## @result{} 2 89 +## 3 90 +## 2 91 +## 2 92 +## 3 100 +## @end group +## @end example +## +## Another example, where the result is a multi-dimensional 3-D array and +## the default value (zero) appears in the output: +## +## @example +## @group +## accumarray ([1, 1, 1; +## 2, 1, 2; +## 2, 3, 2; +## 2, 1, 2; +## 2, 3, 2], 101:105) +## @result{} ans(:,:,1) = [101, 0, 0; 0, 0, 0] +## @result{} ans(:,:,2) = [0, 0, 0; 206, 0, 208] +## @end group +## @end example +## +## The sparse option can be used as an alternative to the @code{sparse} +## constructor (@pxref{XREFsparse,,@code{sparse}}). Thus +## +## @example +## sparse (@var{i}, @var{j}, @var{sv}) +## @end example +## +## @noindent +## can be written with @code{accumarray} as +## +## @example +## accumarray ([@var{i}, @var{j}], @var{sv}', [], [], 0, true) +## @end example +## +## @noindent +## For repeated indices, @code{sparse} adds the corresponding value. To +## take the minimum instead, use @code{min} as an accumulator function: +## +## @example +## accumarray ([@var{i}, @var{j}], @var{sv}', [], @@min, 0, true) +## @end example +## +## The complexity of accumarray in general for the non-sparse case is +## generally O(M+N), where N is the number of subscripts and M is the +## maximum subscript (linearized in multi-dimensional case). If +## @var{fcn} is one of @code{@@sum} (default), @code{@@max}, +## @code{@@min} or @code{@@(x) @{x@}}, an optimized code path is used. +## Note that for general reduction function the interpreter overhead can +## play a major part and it may be more efficient to do multiple +## accumarray calls and compute the results in a vectorized manner. +## +## @seealso{accumdim, unique, sparse} +## @end deftypefn + +function A = accumarray (subs, vals, sz = [], fcn = [], fillval = [], issparse = []) + + if (nargin < 2) + print_usage (); + endif + + lenvals = length (vals); + + if (iscell (subs)) + subs = cellfun (@vec, subs, "uniformoutput", false); + ndims = numel (subs); + if (ndims == 1) + subs = subs{1}; + endif + + lensubs = cellfun (@length, subs); + + if (any (lensubs != lensubs(1)) || (lenvals > 1 && lenvals != lensubs(1))) + error ("accumarray: dimension mismatch"); + endif + + else + ndims = columns (subs); + if (lenvals > 1 && lenvals != rows (subs)) + error ("accumarray: dimension mismatch"); + endif + endif + + if (isempty (fcn)) + fcn = @sum; + elseif (! is_function_handle (fcn)) + error ("accumarray: FCN must be a function handle"); + endif + + if (isempty (fillval)) + fillval = 0; + endif + + if (isempty (issparse)) + issparse = false; + endif + + if (issparse) + + ## Sparse case. + ## Avoid linearizing the subscripts, because it could overflow. + + if (fillval != 0) + error ("accumarray: FILLVAL must be zero in the sparse case"); + endif + + ## Ensure subscripts are a two-column matrix. + if (iscell (subs)) + subs = [subs{:}]; + endif + + ## Validate dimensions. + if (ndims == 1) + subs(:,2) = 1; + elseif (ndims != 2) + error ("accumarray: in the sparse case, needs 1 or 2 subscripts"); + endif + + if (isnumeric (vals) || islogical (vals)) + vals = double (vals); + else + error ("accumarray: in the sparse case, values must be numeric or logical"); + endif + + if (fcn != @sum) + + ## Reduce values. This is not needed if we're about to sum them, + ## because "sparse" can do that. + + ## Sort indices. + [subs, idx] = sortrows (subs); + n = rows (subs); + ## Identify runs. + jdx = find (any (diff (subs, 1, 1), 2)); + jdx = [jdx; n]; + + vals = cellfun (fcn, mat2cell (vals(:)(idx), diff ([0; jdx]))); + subs = subs(jdx, :); + mode = "unique"; + else + mode = "sum"; + endif + + ## Form the sparse matrix. + if (isempty (sz)) + A = sparse (subs(:,1), subs(:,2), vals, mode); + elseif (length (sz) == 2) + + ## Row vector case + if (sz(1) == 1) + [i, j] = deal (subs(:,2), subs(:,1)); + else + [i, j] = deal (subs(:,1), subs(:,2)); + endif + A = sparse (i, j, vals, sz(1), sz(2), mode); + else + error ("accumarray: dimensions mismatch"); + endif + + else + + ## Linearize subscripts. + if (ndims > 1) + if (isempty (sz)) + if (iscell (subs)) + sz = cellfun ("max", subs); + else + sz = max (subs, [], 1); + endif + elseif (ndims != length (sz)) + error ("accumarray: dimensions mismatch"); + endif + + ## Convert multidimensional subscripts. + if (isnumeric (subs)) + subs = num2cell (subs, 1); + endif + subs = sub2ind (sz, subs{:}); # creates index cache + elseif (! isempty (sz) && length (sz) < 2) + error ("accumarray: needs at least 2 dimensions"); + elseif (! isindex (subs)) # creates index cache + error ("accumarray: indices must be positive integers"); + endif + + + ## Some built-in reductions handled efficiently. + + if (fcn == @sum) + ## Fast summation. + if (isempty (sz)) + A = __accumarray_sum__ (subs, vals); + else + A = __accumarray_sum__ (subs, vals, prod (sz)); + ## set proper shape. + A = reshape (A, sz); + endif + + ## we fill in nonzero fill value. + if (fillval != 0) + mask = true (size (A)); + mask(subs) = false; + A(mask) = fillval; + endif + elseif (fcn == @max) + ## Fast maximization. + + if (isinteger (vals)) + zero = intmin (vals); + elseif (islogical (vals)) + zero = false; + elseif (fillval == 0 && all (vals(:) >= 0)) + ## This is a common case - fillval is zero, all numbers + ## nonegative. + zero = 0; + else + zero = NaN; # Neutral value. + endif + + if (isempty (sz)) + A = __accumarray_max__ (subs, vals, zero); + else + A = __accumarray_max__ (subs, vals, zero, prod (sz)); + A = reshape (A, sz); + endif + + if (fillval != zero && ! (isnan (fillval) || isnan (zero))) + mask = true (size (A)); + mask(subs) = false; + A(mask) = fillval; + endif + elseif (fcn == @min) + ## Fast minimization. + + if (isinteger (vals)) + zero = intmax (vals); + elseif (islogical (vals)) + zero = true; + elseif (fillval == 0 && all (vals(:) <= 0)) + ## This is a common case - fillval is zero, all numbers + ## non-positive. + zero = 0; + else + zero = NaN; # Neutral value. + endif + + if (isempty (sz)) + A = __accumarray_min__ (subs, vals, zero); + else + A = __accumarray_min__ (subs, vals, zero, prod (sz)); + A = reshape (A, sz); + endif + + if (fillval != zero && ! (isnan (fillval) || isnan (zero))) + mask = true (size (A)); + mask(subs) = false; + A(mask) = fillval; + endif + else + + ## The general case. Reduce values. + n = rows (subs); + if (numel (vals) == 1) + vals = vals(ones (1, n), 1); + else + vals = vals(:); + endif + + ## Sort indices. + [subs, idx] = sort (subs); + ## Identify runs. + jdx = find (subs(1:n-1) != subs(2:n)); + if (n != 0) # bug #47287 + jdx = [jdx; n]; + endif + vals = mat2cell (vals(idx), diff ([0; jdx])); + ## Optimize the case when function is @(x) {x}, i.e., we just want + ## to collect the values to cells. + persistent simple_cell_str = func2str (@(x) {x}); + if (! strcmp (func2str (fcn), simple_cell_str)) + vals = cellfun (fcn, vals); + endif + + subs = subs(jdx); + + if (isempty (sz)) + sz = max (subs); + ## If subs is empty, sz will be too, and length will be 0, hence "<= 1" + if (length (sz) <= 1) + sz(2) = 1; + endif + endif + + ## Construct matrix of fillvals. + if (iscell (vals)) + A = cell (sz); + elseif (fillval == 0) + A = zeros (sz, class (vals)); + else + A = repmat (fillval, sz); + endif + + ## Set the reduced values. + A(subs) = vals; + endif + endif + +endfunction + + +%!assert (accumarray ([1; 2; 4; 2; 4], 101:105), [101; 206; 0; 208]) +%!assert (accumarray ([1 1 1; 2 1 2; 2 3 2; 2 1 2; 2 3 2], 101:105), +%! cat (3, [101 0 0; 0 0 0], [0 0 0; 206 0 208])) + +%!assert (accumarray ([1 1 1; 2 1 2; 2 3 2; 2 1 2; 2 3 2], 101:105, [], @(x) sin (sum (x))), +%! sin (cat (3, [101,0,0;0,0,0],[0,0,0;206,0,208]))) + +%!assert (accumarray ({[1 3 3 2 3 1 2 2 3 3 1 2], [3 4 2 1 4 3 4 2 2 4 3 4], [1 1 2 2 1 1 2 1 1 1 2 2]}, 101:112), +%! cat (3, [0 0 207 0; 0 108 0 0; 0 109 0 317], [0 0 111 0; 104 0 0 219; 0 103 0 0])) + +%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 101:105, [2 4], @max, NaN), +%! [101 NaN NaN NaN; 104 NaN 105 NaN]) + +%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 101:105, [], @prod), +%! [101 0 0; 10608 0 10815]) +%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 101:105, [2 4], @prod, 0, true), +%! sparse ([1 2 2], [1 1 3], [101 10608 10815], 2, 4)) +%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 1, [2 4]), [1 0 0 0; 2 0 2 0]) +%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 101:105, [2 4], @(x) length (x) > 1), +%! [false false false false; true false true false]) + +%!assert (accumarray ([1; 2], [3; 4], [2, 1], @min, [], 0), [3; 4]) +%!assert (accumarray ([1; 2], [3; 4], [2, 1], @min, [], 1), sparse ([3; 4])) +%!assert (accumarray ([1; 2], [3; 4], [1, 2], @min, [], 0), [3, 4]) +%!assert (accumarray ([1; 2], [3; 4], [1, 2], @min, [], 1), sparse ([3, 4])) + +%!test +%! A = accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 101:105, [2,4], @(x) {x}); +%! assert (A{2},[102; 104]); + +%!test +%! subs = ceil (rand (2000, 3)*10); +%! vals = rand (2000, 1); +%! assert (accumarray (subs, vals, [], @max), +%! accumarray (subs, vals, [], @(x) max (x))); + +%!test +%! subs = ceil (rand (2000, 1)*100); +%! vals = rand (2000, 1); +%! assert (accumarray (subs, vals, [100, 1], @min, NaN), +%! accumarray (subs, vals, [100, 1], @(x) min (x), NaN)); + +%!test +%! subs = ceil (rand (2000, 2)*30); +%! subsc = num2cell (subs, 1); +%! vals = rand (2000, 1); +%! assert (accumarray (subsc, vals, [], [], 0, true), +%! accumarray (subs, vals, [], [], 0, true)); + +%!test +%! subs = ceil (rand (2000, 3)*10); +%! subsc = num2cell (subs, 1); +%! vals = rand (2000, 1); +%! assert (accumarray (subsc, vals, [], @max), +%! accumarray (subs, vals, [], @max)); + +%!error accumarray (1:5) +%!error accumarray ([1,2,3],1:2) + +## Handle empty arrays +%!test <*47287> +%! ## min, max, and sum are special cases within accumarray so test them. +%! funcs = {@(x) length (x) > 1, @min, @max, @sum}; +%! for idx = 1:numel (funcs) +%! assert (accumarray (zeros (0, 1), [], [0 1] , funcs{idx}), zeros (0, 1)); +%! assert (accumarray (zeros (0, 1), [], [1 0] , funcs{idx}), zeros (1, 0)); +%! assert (accumarray (zeros (0, 1), [], [] , funcs{idx}), zeros (0, 1)); +%! endfor + +## Matlab returns an array of doubles even though FCN returns cells. In +## Octave, we do not have that bug, at least for this case. +%!assert (accumarray (zeros (0, 1), [], [0 1] , @(x) {x}), cell (0, 1)) + +%!error +%! accumarray ([1; 2; 3], [1; 2; 3], [3 1], '@(x) {x}') From 0e7bf18d77c4e5ac17b32f6e837fda9366f0ed6f Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Tue, 9 Jul 2024 14:25:25 +0100 Subject: [PATCH 04/21] Add sample files for Octave language --- samples/Octave/BetaDistribution.m | 824 +++++++++++++++++++++++++ samples/Octave/ClusterCriterion.m | 245 ++++++++ samples/Octave/anova1.m | 351 +++++++++++ samples/Octave/silhouette.m | 206 +++++++ samples/Octave/var.m | 956 ++++++++++++++++++++++++++++++ 5 files changed, 2582 insertions(+) create mode 100644 samples/Octave/BetaDistribution.m create mode 100644 samples/Octave/ClusterCriterion.m create mode 100644 samples/Octave/anova1.m create mode 100644 samples/Octave/silhouette.m create mode 100644 samples/Octave/var.m diff --git a/samples/Octave/BetaDistribution.m b/samples/Octave/BetaDistribution.m new file mode 100644 index 0000000000..c71509099d --- /dev/null +++ b/samples/Octave/BetaDistribution.m @@ -0,0 +1,824 @@ +## Copyright (C) 2024 Andreas Bertsatos +## +## This file is part of the statistics package for GNU Octave. +## +## 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 . + +classdef BetaDistribution + ## -*- texinfo -*- + ## @deftypefn {statistics} BetaDistribution + ## + ## Beta probability distribution object. + ## + ## A @code{BetaDistribution} object consists of parameters, a model + ## description, and sample data for a beta probability distribution. + ## + ## The beta distribution uses the following parameters. + ## + ## @multitable @columnfractions 0.25 0.48 0.27 + ## @headitem @var{Parameter} @tab @var{Description} @tab @var{Support} + ## + ## @item @qcode{a} @tab 1st Shape parameter @tab @math{α > 0} + ## @item @qcode{b} @tab 2nd Shape parameter @tab @math{β > 0} + ## @end multitable + ## + ## There are several ways to create a @code{BetaDistribution} object. + ## + ## @itemize + ## @item Fit a distribution to data using the @code{fitdist} function. + ## @item Create a distribution with specified parameter values using the + ## @code{makedist} function. + ## @item Use the constructor @qcode{BetaDistribution (@var{a}, @var{b})} + ## to create a beta distribution with specified parameter values. + ## @item Use the static method @qcode{BetaDistribution.fit (@var{x}, + ## @var{censor}, @var{freq}, @var{options})} to a distribution to data @var{x}. + ## @end itemize + ## + ## It is highly recommended to use @code{fitdist} and @code{makedist} + ## functions to create probability distribution objects, instead of the + ## constructor and the aforementioned static method. + ## + ## A @code{BetaDistribution} object contains the following properties, + ## which can be accessed using dot notation. + ## + ## @multitable @columnfractions 0.25 0.25 0.25 0.25 + ## @item @qcode{DistributionName} @tab @qcode{DistributionCode} @tab + ## @qcode{NumParameters} @tab @qcode{ParameterNames} + ## @item @qcode{ParameterDescription} @tab @qcode{ParameterValues} @tab + ## @qcode{ParameterValues} @tab @qcode{ParameterCI} + ## @item @qcode{ParameterIsFixed} @tab @qcode{Truncation} @tab + ## @qcode{IsTruncated} @tab @qcode{InputData} + ## @end multitable + ## + ## A @code{BetaDistribution} object contains the following methods: + ## @code{cdf}, @code{icdf}, @code{iqr}, @code{mean}, @code{median}, + ## @code{negloglik}, @code{paramci}, @code{pdf}, @code{plot}, @code{proflik}, + ## @code{random}, @code{std}, @code{truncate}, @code{var}. + ## + ## Further information about the Beta distribution can be found at + ## @url{https://en.wikipedia.org/wiki/Beta_distribution} + ## + ## @seealso{fitdist, makedist, betacdf, betainv, betapdf, betarnd, lognfit, + ## betalike, betastat} + ## @end deftypefn + + properties (Dependent = true) + a + b + endproperties + + properties (GetAccess = public, Constant = true) + CensoringAllowed = false; + DistributionName = "BetaDistribution"; + DistributionCode = "beta"; + NumParameters = 2; + ParameterNames = {"a", "b"}; + ParameterDescription = {"1st Shape", "2nd Shape"}; + endproperties + + properties (GetAccess = public, Constant = true) + ParameterRange = [realmin, realmin; Inf, Inf]; + ParameterLogCI = [true, true]; + endproperties + + properties (GetAccess = public , SetAccess = protected) + ParameterValues + ParameterCI + ParameterCovariance + ParameterIsFixed + Truncation + IsTruncated + InputData + endproperties + + methods (Hidden) + + function this = BetaDistribution (a, b) + if (nargin == 0) + a = 1; + b = 1; + endif + checkparams (a, b); + this.InputData = []; + this.IsTruncated = false; + this.ParameterValues = [a, b]; + this.ParameterIsFixed = [true, true]; + this.ParameterCovariance = zeros (this.NumParameters); + endfunction + + function display (this) + fprintf ("%s =\n", inputname(1)); + __disp__ (this, "beta distribution"); + endfunction + + function disp (this) + __disp__ (this, "beta distribution"); + endfunction + + function this = set.a (this, a) + checkparams (a, this.b); + this.InputData = []; + this.ParameterValues(1) = a; + this.ParameterIsFixed = [true, true]; + this.ParameterCovariance = zeros (this.NumParameters); + endfunction + + function a = get.a (this) + a = this.ParameterValues(1); + endfunction + + function this = set.b (this, b) + checkparams (this.a, b); + this.InputData = []; + this.ParameterValues(2) = b; + this.ParameterIsFixed = [true, true]; + this.ParameterCovariance = zeros (this.NumParameters); + endfunction + + function b = get.b (this) + b = this.ParameterValues(2); + endfunction + + endmethods + + methods (Access = public) + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{p} =} cdf (@var{pd}, @var{x}) + ## @deftypefnx {BetaDistribution} {@var{p} =} cdf (@var{pd}, @var{x}, @qcode{"upper"}) + ## + ## Compute the cumulative distribution function (CDF). + ## + ## @code{@var{p} = cdf (@var{pd}, @var{x})} computes the CDF of the + ## probability distribution object, @var{pd}, evaluated at the values in + ## @var{x}. + ## + ## @code{@var{p} = cdf (@dots{}, @qcode{"upper"})} returns the complement of + ## the CDF of the probability distribution object, @var{pd}, evaluated at + ## the values in @var{x}. + ## + ## @end deftypefn + function p = cdf (this, x, uflag) + if (! isscalar (this)) + error ("cdf: requires a scalar probability distribution."); + endif + ## Check for "upper" flag + if (nargin > 2 && strcmpi (uflag, "upper")) + utail = true; + elseif (nargin > 2 && ! strcmpi (uflag, "upper")) + error ("cdf: invalid argument for upper tail."); + else + utail = false; + endif + ## Do the computations + p = betacdf (x, this.a, this.b); + if (this.IsTruncated) + lx = this.Truncation(1); + lb = x < lx; + ux = this.Truncation(2); + ub = x > ux; + p(lb) = 0; + p(ub) = 1; + p(! (lb | ub)) -= betacdf (lx, this.a, this.b); + p(! (lb | ub)) /= diff (betacdf ([lx, ux], this.a, this.b)); + endif + ## Apply uflag + if (utail) + p = 1 - p; + endif + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{p} =} icdf (@var{pd}, @var{p}) + ## + ## Compute the cumulative distribution function (CDF). + ## + ## @code{@var{p} = icdf (@var{pd}, @var{x})} computes the quantile (the + ## inverse of the CDF) of the probability distribution object, @var{pd}, + ## evaluated at the values in @var{x}. + ## + ## @end deftypefn + function x = icdf (this, p) + if (! isscalar (this)) + error ("icdf: requires a scalar probability distribution."); + endif + if (this.IsTruncated) + lp = betacdf (this.Truncation(1), this.a, this.b); + up = betacdf (this.Truncation(2), this.a, this.b); + ## Adjust p values within range of p @ lower limit and p @ upper limit + is_nan = p < 0 | p > 1; + p(is_nan) = NaN; + np = lp + (up - lp) .* p; + x = betainv (np, this.a, this.b); + x(x < this.Truncation(1)) = this.Truncation(1); + x(x > this.Truncation(2)) = this.Truncation(2); + else + x = betainv (p, this.a, this.b); + endif + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{r} =} iqr (@var{pd}) + ## + ## Compute the interquartile range of a probability distribution. + ## + ## @code{@var{r} = iqr (@var{pd})} computes the interquartile range of the + ## probability distribution object, @var{pd}. + ## + ## @end deftypefn + function r = iqr (this) + if (! isscalar (this)) + error ("iqr: requires a scalar probability distribution."); + endif + r = diff (icdf (this, [0.25, 0.75])); + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{m} =} mean (@var{pd}) + ## + ## Compute the mean of a probability distribution. + ## + ## @code{@var{m} = mean (@var{pd})} computes the mean of the probability + ## distribution object, @var{pd}. + ## + ## @end deftypefn + function m = mean (this) + if (! isscalar (this)) + error ("mean: requires a scalar probability distribution."); + endif + if (this.IsTruncated) + fm = @(x) x .* pdf (this, x); + m = integral (fm, this.Truncation(1), this.Truncation(2)); + else + m = betastat (this.a, this.b); + endif + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{m} =} median (@var{pd}) + ## + ## Compute the median of a probability distribution. + ## + ## @code{@var{m} = median (@var{pd})} computes the median of the probability + ## distribution object, @var{pd}. + ## + ## @end deftypefn + function m = median (this) + if (! isscalar (this)) + error ("median: requires a scalar probability distribution."); + endif + if (this.IsTruncated) + lx = this.Truncation(1); + ux = this.Truncation(2); + Fa_b = betacdf ([lx, ux], this.a, this.b); + m = betainv (sum (Fa_b) / 2, this.a, this.b); + else + m = betainv (0.5, this.a, this.b); + endif + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{nlogL} =} negloglik (@var{pd}) + ## + ## Compute the negative loglikelihood of a probability distribution. + ## + ## @code{@var{m} = negloglik (@var{pd})} computes the negative loglikelihood + ## of the probability distribution object, @var{pd}. + ## + ## @end deftypefn + function nlogL = negloglik (this) + if (! isscalar (this)) + error ("negloglik: requires a scalar probability distribution."); + endif + if (isempty (this.InputData)) + nlogL = []; + return + endif + nlogL = - betalike ([this.a, this.b], this.InputData.data, ... + this.InputData.freq); + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{ci} =} paramci (@var{pd}) + ## @deftypefnx {BetaDistribution} {@var{ci} =} paramci (@var{pd}, @var{Name}, @var{Value}) + ## + ## Compute the confidence intervals for probability distribution parameters. + ## + ## @code{@var{ci} = paramci (@var{pd})} computes the lower and upper + ## boundaries of the 95% confidence interval for each parameter of the + ## probability distribution object, @var{pd}. + ## + ## @code{@var{ci} = paramci (@var{pd}, @var{Name}, @var{Value})} computes the + ## confidence intervals with additional options specified specified by + ## @qcode{Name-Value} pair arguments listed below. + ## + ## @multitable @columnfractions 0.18 0.02 0.8 + ## @headitem @var{Name} @tab @tab @var{Value} + ## + ## @item @qcode{"Alpha"} @tab @tab A scalar value in the range @math{(0,1)} + ## specifying the significance level for the confidence interval. The + ## default value 0.05 corresponds to a 95% confidence interval. + ## + ## @item @qcode{"Parameter"} @tab @tab A character vector or a cell array of + ## character vectors specifying the parameter names for which to compute + ## confidence intervals. By default, @code{paramci} computes confidence + ## intervals for all distribution parameters. + ## @end multitable + ## + ## @code{paramci} is meaningful only when @var{pd} is fitted to data, + ## otherwise an empty array, @qcode{[]}, is returned. + ## + ## @end deftypefn + function ci = paramci (this, varargin) + if (! isscalar (this)) + error ("paramci: requires a scalar probability distribution."); + endif + if (isempty (this.InputData)) + ci = [this.ParameterValues; this.ParameterValues]; + else + ci = __paramci__ (this, varargin{:}); + endif + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{y} =} pdf (@var{pd}, @var{x}) + ## + ## Compute the probability distribution function (PDF). + ## + ## @code{@var{y} = pdf (@var{pd}, @var{x})} computes the PDF of the + ## probability distribution object, @var{pd}, evaluated at the values in + ## @var{x}. + ## + ## @end deftypefn + function y = pdf (this, x) + if (! isscalar (this)) + error ("pdf: requires a scalar probability distribution."); + endif + y = betapdf (x, this.a, this.b); + if (this.IsTruncated) + lx = this.Truncation(1); + lb = x < lx; + ux = this.Truncation(2); + ub = x > ux; + y(lb | ub) = 0; + y(! (lb | ub)) /= diff (betacdf ([lx, ux], this.a, this.b)); + endif + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {} plot (@var{pd}) + ## @deftypefnx {BetaDistribution} {} plot (@var{pd}, @var{Name}, @var{Value}) + ## @deftypefnx {BetaDistribution} {@var{h} =} plot (@dots{}) + ## + ## Plot a probability distribution object. + ## + ## @code{plot (@var{pd}} plots a probability density function (PDF) of the + ## probability distribution object @var{pd}. If @var{pd} contains data, + ## which have been fitted by @code{fitdist}, the PDF is superimposed over a + ## histogram of the data. + ## + ## @code{plot (@var{pd}, @var{Name}, @var{Value})} specifies additional + ## options with the @qcode{Name-Value} pair arguments listed below. + ## + ## @multitable @columnfractions 0.18 0.02 0.8 + ## @headitem @tab @var{Name} @tab @var{Value} + ## + ## @item @qcode{"PlotType"} @tab @tab A character vector specifying the plot + ## type. @qcode{"pdf"} plots the probability density function (PDF). When + ## @var{pd} is fit to data, the PDF is superimposed on a histogram of the + ## data. @qcode{"cdf"} plots the cumulative density function (CDF). When + ## @var{pd} is fit to data, the CDF is superimposed over an empirical CDF. + ## @qcode{"probability"} plots a probability plot using a CDF of the data + ## and a CDF of the fitted probability distribution. This option is + ## available only when @var{pd} is fitted to data. + ## + ## @item @qcode{"Discrete"} @tab @tab A logical scalar to specify whether to + ## plot the PDF or CDF of a discrete distribution object as a line plot or a + ## stem plot, by specifying @qcode{false} or @qcode{true}, respectively. By + ## default, it is @qcode{true} for discrete distributions and @qcode{false} + ## for continuous distributions. When @var{pd} is a continuous distribution + ## object, option is ignored. + ## + ## @item @qcode{"Parent"} @tab @tab An axes graphics object for plot. If + ## not specified, the @code{plot} function plots into the current axes or + ## creates a new axes object if one does not exist. + ## @end multitable + ## + ## @code{@var{h} = plot (@dots{})} returns a graphics handle to the plotted + ## objects. + ## + ## @end deftypefn + function [varargout] = plot (this, varargin) + if (! isscalar (this)) + error ("plot: requires a scalar probability distribution."); + endif + h = __plot__ (this, false, varargin{:}); + if (nargout > 0) + varargout{1} = h; + endif + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {[@var{nlogL}, @var{param}] =} proflik (@var{pd}, @var{pnum}) + ## @deftypefnx {BetaDistribution} {[@var{nlogL}, @var{param}] =} proflik (@var{pd}, @var{pnum}, @qcode{"Display"}, @var{display}) + ## @deftypefnx {BetaDistribution} {[@var{nlogL}, @var{param}] =} proflik (@var{pd}, @var{pnum}, @var{setparam}) + ## @deftypefnx {BetaDistribution} {[@var{nlogL}, @var{param}] =} proflik (@var{pd}, @var{pnum}, @var{setparam}, @qcode{"Display"}, @var{display}) + ## + ## Profile likelihood function for a probability distribution object. + ## + ## @code{[@var{nlogL}, @var{param}] = proflik (@var{pd}, @var{pnum})} + ## returns a vector @var{nlogL} of negative loglikelihood values and a + ## vector @var{param} of corresponding parameter values for the parameter in + ## the position indicated by @var{pnum}. By default, @code{proflik} uses + ## the lower and upper bounds of the 95% confidence interval and computes + ## 100 equispaced values for the selected parameter. @var{pd} must be + ## fitted to data. + ## + ## @code{[@var{nlogL}, @var{param}] = proflik (@var{pd}, @var{pnum}, + ## @qcode{"Display"}, @qcode{"on"})} also plots the profile likelihood + ## against the default range of the selected parameter. + ## + ## @code{[@var{nlogL}, @var{param}] = proflik (@var{pd}, @var{pnum}, + ## @var{setparam})} defines a user-defined range of the selected parameter. + ## + ## @code{[@var{nlogL}, @var{param}] = proflik (@var{pd}, @var{pnum}, + ## @var{setparam}, @qcode{"Display"}, @qcode{"on"})} also plots the profile + ## likelihood against the user-defined range of the selected parameter. + ## + ## For the beta distribution, @qcode{@var{pnum} = 1} selects the parameter + ## @qcode{a} and @qcode{@var{pnum} = 2} selects the parameter @var{b}. + ## + ## When opted to display the profile likelihood plot, @code{proflik} also + ## plots the baseline loglikelihood computed at the lower bound of the 95% + ## confidence interval and estimated maximum likelihood. The latter might + ## not be observable if it is outside of the used-defined range of parameter + ## values. + ## + ## @end deftypefn + function [varargout] = proflik (this, pnum, varargin) + if (! isscalar (this)) + error ("proflik: requires a scalar probability distribution."); + endif + if (isempty (this.InputData)) + error ("proflik: no fitted data available."); + endif + [varargout{1:nargout}] = __proflik__ (this, pnum, varargin{:}); + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{y} =} random (@var{pd}) + ## @deftypefnx {BetaDistribution} {@var{y} =} random (@var{pd}, @var{rows}) + ## @deftypefnx {BetaDistribution} {@var{y} =} random (@var{pd}, @var{rows}, @var{cols}, @dots{}) + ## @deftypefnx {BetaDistribution} {@var{y} =} random (@var{pd}, [@var{sz}]) + ## + ## Generate random arrays from the probability distribution object. + ## + ## @code{@var{r} = random (@var{pd})} returns a random number from the + ## distribution object @var{pd}. + ## + ## When called with a single size argument, @code{betarnd} returns a square + ## matrix with the dimension specified. When called with more than one + ## scalar argument, the first two arguments are taken as the number of rows + ## and columns and any further arguments specify additional matrix + ## dimensions. The size may also be specified with a row vector of + ## dimensions, @var{sz}. + ## + ## @end deftypefn + function r = random (this, varargin) + if (! isscalar (this)) + error ("random: requires a scalar probability distribution."); + endif + if (this.IsTruncated) + sz = [varargin{:}]; + ps = prod (sz); + ## Get an estimate of how many more random numbers we need to randomly + ## pick the appropriate size from + lx = this.Truncation(1); + ux = this.Truncation(2); + ratio = 1 / diff (betacdf ([lx, ux], this.a, this.b)); + nsize = fix (2 * ratio * ps); # times 2 to be on the safe side + ## Generate the numbers and remove out-of-bound random samples + r = betarnd (this.a, this.b, nsize, 1); + r(r < lx | r > ux) = []; + ## Randomly select the required size and reshape to requested dimensions + idx = randperm (numel (r), ps); + r = reshape (r(idx), sz); + else + r = betarnd (this.a, this.b, varargin{:}); + endif + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{s} =} std (@var{pd}) + ## + ## Compute the standard deviation of a probability distribution. + ## + ## @code{@var{s} = std (@var{pd})} computes the standard deviation of the + ## probability distribution object, @var{pd}. + ## + ## @end deftypefn + function s = std (this) + if (! isscalar (this)) + error ("std: requires a scalar probability distribution."); + endif + v = var (this); + s = sqrt (v); + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{t} =} truncate (@var{pd}, @var{lower}, @var{upper}) + ## + ## Truncate a probability distribution. + ## + ## @code{@var{t} = truncate (@var{pd})} returns a probability distribution + ## @var{t}, which is the probability distribution @var{pd} truncated to the + ## specified interval with lower limit, @var{lower}, and upper limit, + ## @var{upper}. If @var{pd} is fitted to data with @code{fitdist}, the + ## returned probability distribution @var{t} is not fitted, does not contain + ## any data or estimated values, and it is as it has been created with the + ## @var{makedist} function, but it includes the truncation interval. + ## + ## @end deftypefn + function this = truncate (this, lower, upper) + if (! isscalar (this)) + error ("truncate: requires a scalar probability distribution."); + endif + if (nargin < 3) + error ("truncate: missing input argument."); + elseif (lower >= upper) + error ("truncate: invalid lower upper limits."); + endif + this.Truncation = [lower, upper]; + this.IsTruncated = true; + this.InputData = []; + this.ParameterIsFixed = [true, true]; + this.ParameterCovariance = zeros (this.NumParameters); + endfunction + + ## -*- texinfo -*- + ## @deftypefn {BetaDistribution} {@var{v} =} var (@var{pd}) + ## + ## Compute the variance of a probability distribution. + ## + ## @code{@var{v} = var (@var{pd})} computes the standard deviation of the + ## probability distribution object, @var{pd}. + ## + ## @end deftypefn + function v = var (this) + if (! isscalar (this)) + error ("var: requires a scalar probability distribution."); + endif + if (this.IsTruncated) + fm = @(x) x .* pdf (this, x); + m = integral (fm, this.Truncation(1), this.Truncation(2)); + fv = @(x) ((x - m) .^ 2) .* pdf (this, x); + v = integral (fv, this.Truncation(1), this.Truncation(2)); + else + [~, v] = betastat (this.a, this.b); + endif + endfunction + + endmethods + + methods (Static, Hidden) + + function pd = fit (x, varargin) + ## Check input arguments + if (nargin < 2) + alpha = 0.05; + else + alpha = varargin{1}; + endif + if (nargin < 3) + freq = []; + else + freq = varargin{2}; + endif + if (nargin < 4) + options.Display = "off"; + options.MaxFunEvals = 400; + options.MaxIter = 200; + options.TolX = 1e-6; + else + options = varargin{3}; + endif + ## Fit data + [phat, pci] = betafit (x, alpha, freq, options); + [~, acov] = betalike (phat, x, freq); + ## Create fitted distribution object + pd = BetaDistribution.makeFitted (phat, pci, acov, x, freq); + endfunction + + function pd = makeFitted (phat, pci, acov, x, freq) + a = phat(1); + b = phat(2); + pd = BetaDistribution (a, b); + pd.ParameterCI = pci; + pd.ParameterIsFixed = [false, false]; + pd.ParameterCovariance = acov; + pd.InputData = struct ("data", x, "cens", [], "freq", freq); + endfunction + + endmethods + +endclassdef + +function checkparams (a, b) + if (! (isscalar (a) && isnumeric (a) && isreal (a) && isfinite (a) + && a > 0)) + error ("BetaDistribution: A must be a positive real scalar.") + endif + if (! (isscalar (b) && isnumeric (b) && isreal (b) + && isfinite (b) && b > 0)) + error ("BetaDistribution: B must be a positive real scalar.") + endif +endfunction + +## Test output +%!shared pd, t +%! pd = BetaDistribution; +%! t = truncate (pd, 0.2, 0.8); +%!assert (cdf (pd, [0:0.2:1]), [0, 0.2, 0.4, 0.6, 0.8, 1], 1e-4); +%!assert (cdf (t, [0:0.2:1]), [0, 0, 0.3333, 0.6667, 1, 1], 1e-4); +%!assert (cdf (pd, [-1, 1, NaN]), [0, 1, NaN], 1e-4); +%!assert (cdf (t, [-1, 1, NaN]), [0, 1, NaN], 1e-4); +%!assert (icdf (pd, [0:0.2:1]), [0, 0.2, 0.4, 0.6, 0.8, 1], 1e-4); +%!assert (icdf (t, [0:0.2:1]), [0.2, 0.32, 0.44, 0.56, 0.68, 0.8], 1e-4); +%!assert (icdf (pd, [-1, 0.4:0.2:1, NaN]), [NaN, 0.4, 0.6, 0.8, 1, NaN], 1e-4); +%!assert (icdf (t, [-1, 0.4:0.2:1, NaN]), [NaN, 0.44, 0.56, 0.68, 0.8, NaN], 1e-4); +%!assert (iqr (pd), 0.5, 1e-4); +%!assert (iqr (t), 0.3, 1e-4); +%!assert (mean (pd), 0.5); +%!assert (mean (t), 0.5, 1e-6); +%!assert (median (pd), 0.5); +%!assert (median (t), 0.5, 1e-6); +%!assert (pdf (pd, [0:0.2:1]), [1, 1, 1, 1, 1, 1], 1e-4); +%!assert (pdf (t, [0:0.2:1]), [0, 1.6667, 1.6667, 1.6667, 1.6667, 0], 1e-4); +%!assert (pdf (pd, [-1, 1, NaN]), [0, 1, NaN], 1e-4); +%!assert (pdf (t, [-1, 1, NaN]), [0, 0, NaN], 1e-4); +%!assert (isequal (size (random (pd, 100, 50)), [100, 50])) +%!assert (any (random (t, 1000, 1) < 0.2), false); +%!assert (any (random (t, 1000, 1) > 0.8), false); +%!assert (std (pd), 0.2887, 1e-4); +%!assert (std (t), 0.1732, 1e-4); +%!assert (var (pd), 0.0833, 1e-4); +%!assert (var (t), 0.0300, 1e-4); + +## Test input validation +## 'BetaDistribution' constructor +%!error ... +%! BetaDistribution(0, 1) +%!error ... +%! BetaDistribution(Inf, 1) +%!error ... +%! BetaDistribution(i, 1) +%!error ... +%! BetaDistribution("a", 1) +%!error ... +%! BetaDistribution([1, 2], 1) +%!error ... +%! BetaDistribution(NaN, 1) +%!error ... +%! BetaDistribution(1, 0) +%!error ... +%! BetaDistribution(1, -1) +%!error ... +%! BetaDistribution(1, Inf) +%!error ... +%! BetaDistribution(1, i) +%!error ... +%! BetaDistribution(1, "a") +%!error ... +%! BetaDistribution(1, [1, 2]) +%!error ... +%! BetaDistribution(1, NaN) + +## 'cdf' method +%!error ... +%! cdf (BetaDistribution, 2, "uper") +%!error ... +%! cdf (BetaDistribution, 2, 3) + +## 'paramci' method +%!shared x +%! randg ("seed", 1); +%! x = betarnd (1, 1, [100, 1]); +%!error ... +%! paramci (BetaDistribution.fit (x), "alpha") +%!error ... +%! paramci (BetaDistribution.fit (x), "alpha", 0) +%!error ... +%! paramci (BetaDistribution.fit (x), "alpha", 1) +%!error ... +%! paramci (BetaDistribution.fit (x), "alpha", [0.5 2]) +%!error ... +%! paramci (BetaDistribution.fit (x), "alpha", "") +%!error ... +%! paramci (BetaDistribution.fit (x), "alpha", {0.05}) +%!error ... +%! paramci (BetaDistribution.fit (x), "parameter", "a", "alpha", {0.05}) +%!error ... +%! paramci (BetaDistribution.fit (x), "parameter", {"a", "b", "param"}) +%!error ... +%! paramci (BetaDistribution.fit (x), "alpha", 0.01, ... +%! "parameter", {"a", "b", "param"}) +%!error ... +%! paramci (BetaDistribution.fit (x), "parameter", "param") +%!error ... +%! paramci (BetaDistribution.fit (x), "alpha", 0.01, "parameter", "param") +%!error ... +%! paramci (BetaDistribution.fit (x), "NAME", "value") +%!error ... +%! paramci (BetaDistribution.fit (x), "alpha", 0.01, "NAME", "value") +%!error ... +%! paramci (BetaDistribution.fit (x), "alpha", 0.01, "parameter", "a", ... +%! "NAME", "value") + +## 'plot' method +%!error ... +%! plot (BetaDistribution, "Parent") +%!error ... +%! plot (BetaDistribution, "PlotType", 12) +%!error ... +%! plot (BetaDistribution, "PlotType", {"pdf", "cdf"}) +%!error ... +%! plot (BetaDistribution, "PlotType", "pdfcdf") +%!error ... +%! plot (BetaDistribution, "Discrete", "pdfcdf") +%!error ... +%! plot (BetaDistribution, "Discrete", [1, 0]) +%!error ... +%! plot (BetaDistribution, "Discrete", {true}) +%!error ... +%! plot (BetaDistribution, "Parent", 12) +%!error ... +%! plot (BetaDistribution, "Parent", "hax") + +## 'proflik' method +%!error ... +%! proflik (BetaDistribution, 2) +%!error ... +%! proflik (BetaDistribution.fit (x), 3) +%!error ... +%! proflik (BetaDistribution.fit (x), [1, 2]) +%!error ... +%! proflik (BetaDistribution.fit (x), {1}) +%!error ... +%! proflik (BetaDistribution.fit (x), 1, ones (2)) +%!error ... +%! proflik (BetaDistribution.fit (x), 1, "Display") +%!error ... +%! proflik (BetaDistribution.fit (x), 1, "Display", 1) +%!error ... +%! proflik (BetaDistribution.fit (x), 1, "Display", {1}) +%!error ... +%! proflik (BetaDistribution.fit (x), 1, "Display", {"on"}) +%!error ... +%! proflik (BetaDistribution.fit (x), 1, "Display", ["on"; "on"]) +%!error ... +%! proflik (BetaDistribution.fit (x), 1, "Display", "onnn") +%!error ... +%! proflik (BetaDistribution.fit (x), 1, "NAME", "on") +%!error ... +%! proflik (BetaDistribution.fit (x), 1, {"NAME"}, "on") +%!error ... +%! proflik (BetaDistribution.fit (x), 1, {[1 2 3 4]}, "Display", "on") + +## 'truncate' method +%!error ... +%! truncate (BetaDistribution) +%!error ... +%! truncate (BetaDistribution, 2) +%!error ... +%! truncate (BetaDistribution, 4, 2) + +## Catch errors when using array of probability objects with available methods +%!shared pd +%! pd = BetaDistribution(1, 1); +%! pd(2) = BetaDistribution(1, 3); +%!error cdf (pd, 1) +%!error icdf (pd, 0.5) +%!error iqr (pd) +%!error mean (pd) +%!error median (pd) +%!error negloglik (pd) +%!error paramci (pd) +%!error pdf (pd, 1) +%!error plot (pd) +%!error proflik (pd, 2) +%!error random (pd) +%!error std (pd) +%!error ... +%! truncate (pd, 2, 4) +%!error var (pd) diff --git a/samples/Octave/ClusterCriterion.m b/samples/Octave/ClusterCriterion.m new file mode 100644 index 0000000000..1dd31b243a --- /dev/null +++ b/samples/Octave/ClusterCriterion.m @@ -0,0 +1,245 @@ +## Copyright (C) 2021 Stefano Guidoni +## Copyright (C) 2024 Andreas Bertsatos +## +## This file is part of the statistics package for GNU Octave. +## +## 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 . + +classdef ClusterCriterion < handle + ## -*- texinfo -*- + ## @deftypefn {statistics} {@var{obj} =} ClusterCriterion (@var{x}, @var{clust}, @var{criterion}) + ## + ## A clustering evaluation object as created by @code{evalclusters}. + ## + ## @code{ClusterCriterion} is a superclass for clustering evaluation objects + ## as created by @code{evalclusters}. + ## + ## List of public properties: + ## @table @code + ## @item @qcode{ClusteringFunction} + ## a valid clustering funtion name or function handle. It can be empty if + ## the clustering solutions are passed as an input matric. + ## + ## @item @qcode{CriterionName} + ## a valid criterion name to evaluate the clustering solutions. + ## + ## @item @qcode{CriterionValues} + ## a vector of values as generated by the evaluation criterion for each + ## clustering solution. + ## + ## @item @qcode{InspectedK} + ## the list of proposed cluster numbers. + ## + ## @item @qcode{Missing} + ## a logical vector of missing observations. When there are @code{NaN} + ## values in the data matrix, the corresponding observation is excluded. + ## + ## @item @qcode{NumObservations} + ## the number of non-missing observations in the data matrix. + ## + ## @item @qcode{OptimalK} + ## the optimal number of clusters. + ## + ## @item @qcode{OptimalY} + ## the clustering solution corresponding to @code{OptimalK}. + ## + ## @item @qcode{X} + ## the data matrix. + ## + ## @end table + ## + ## List of public methods: + ## @table @code + ## @item @qcode{addK} + ## add a list of numbers of clusters to evaluate. + ## + ## @item @qcode{compact} + ## return a compact clustering evaluation object. Not implemented + ## + ## @item @qcode{plot} + ## plot the clustering evaluation values against the corresponding number of + ## clusters. + ## + ## @end table + ## + ## @seealso{evalclusters, CalinskiHarabaszEvaluation, DaviesBouldinEvaluation, + ## GapEvaluation, SilhouetteEvaluation} + ## @end deftypefn + + properties (Access = public) + ## public properties + endproperties + + properties (GetAccess = public, SetAccess = protected) + ClusteringFunction = ""; + CriterionName = ""; + CriterionValues = []; + InspectedK = []; + Missing = []; + NumObservations = 0; + OptimalK = 0; + OptimalY = []; + X = []; + endproperties + + properties (Access = protected) + N = 0; # number of observations + P = 0; # number of variables + ClusteringSolutions = []; # + OptimalIndex = 0; # index of the optimal K + endproperties + + methods (Access = public) + + ## constructor + function this = ClusterCriterion (x, clust, KList) + ## parsing input data + if ((! ismatrix (x)) || (! isnumeric (x))) + error ("ClusterCriterion: 'x' must be a numeric matrix"); + endif + this.X = x; + this.N = rows (this.X); + this.P = columns (this.X); + ## look for missing values + for iter = 1 : this.N + if (any (find (x(iter, :) == NaN))) + this.Missing(iter) = true; + else + this.Missing(iter) = false; + endif + endfor + ## number of usable observations + this.NumObservations = sum (this.Missing == false); + + ## parsing the clustering algorithm + if (ischar (clust)) + if (any (strcmpi (clust, {"kmeans", "linkage", "gmdistribution"}))) + this.ClusteringFunction = lower (clust); + else + error ("ClusterCriterion: unknown clustering algorithm '%s'", clust); + endif + elseif (isa (clust, "function_handle")) + this.ClusteringFunction = clust; + elseif (ismatrix (clust)) + if (isnumeric (clust) && (length (size (clust)) == 2) && ... + (rows (clust) == this.N)) + this.ClusteringFunction = ""; + this.ClusteringSolutions = clust(find (this.Missing == false), :); + else + error ("ClusterCriterion: invalid matrix of clustering solutions"); + endif + else + error ("ClusterCriterion: invalid argument"); + endif + + ## parsing the list of cluster sizes to inspect + this.InspectedK = parseKList (this, KList); + endfunction + + ## -*- texinfo -*- + ## @deftypefn {ClusterCriterion} {@var{obj} =} addK (@var{obj}, @var{K}) + ## + ## Add a new cluster array to inspect the ClusterCriterion object. + ## + ## @end deftypefn + function this = addK (this, k) + + ## if there is not a clustering function, then we are using a predefined + ## set of clustering solutions, hence we cannot redefine the number of + ## solutions + if (isempty (this.ClusteringFunction)) + warning (["ClusterCriterion.addK: cannot redefine the list of cluster"... + "numbers to evaluate when there is not a clustering function"]); + return; + endif + + ## otherwise go on + newList = this.parseKList ([this.InspectedK k]); + + ## check if the list has changed + if (length (newList) == length (this.InspectedK)) + warning ("ClusterCriterion.addK: the list has not changed"); + else + ## update ClusteringSolutions and CriterionValues + ClusteringSolutions_tmp = zeros (this.NumObservations, ... + length (newList)); + CriterionValues_tmp = zeros (length (newList), 1); + for iter = 1 : length (this.InspectedK) + idx = find (newList == this.InspectedK(iter)); + + if (! isempty (idx)) + ClusteringSolutions_tmp(:, idx) = this.ClusteringSolutions(:, iter); + CriterionValues_tmp(idx) = this.CriterionValues(iter); + endif + endfor + this.ClusteringSolutions = ClusteringSolutions_tmp; + this.CriterionValues = CriterionValues_tmp; + + ## reset the old results + this.OptimalK = 0; + this.OptimalY = []; + this.OptimalIndex = 0; + + ## update the list of cluster numbers to evaluate + this.InspectedK = newList; + endif + endfunction + + ## -*- texinfo -*- + ## @deftypefn {ClusterCriterion} {} plot (@var{obj}) + ## @deftypefnx {ClusterCriterion} {@var{h} =} plot (@var{obj}) + ## + ## Plot the evaluation results. + ## + ## Plot the CriterionValues against InspectedK from the ClusterCriterion, + ## @var{obj}, to the current plot. It can also return a handle to the + ## current plot. + ## + ## @end deftypefn + function h = plot (this) + yLabel = sprintf ("%s value", this.CriterionName); + h = gca (); + hold on; + plot (this.InspectedK, this.CriterionValues, "bo-"); + plot (this.OptimalK, this.CriterionValues(this.OptimalIndex), "b*"); + xlabel ("number of clusters"); + ylabel (yLabel); + hold off; + endfunction + + ## -*- texinfo -*- + ## @deftypefn {ClusterCriterion} {@var{obj} =} compact (@var{obj}) + ## + ## Return a compact ClusterCriterion object (not implemented yet). + ## + ## @end deftypefn + function this = compact (this) + warning ("ClusterCriterion.compact: this method is not yet implemented."); + endfunction + + endmethods + + methods (Access = private) + ## check if a list of cluster sizes is correct + function retList = parseKList (this, KList) + if (isnumeric (KList) && isvector (KList) && all (find (KList > 0)) && ... + all (floor (KList) == KList)) + retList = unique (KList); + else + error (["ClusterCriterion: the list of cluster sizes must be an " ... + "array of positive integer numbers"]); + endif + endfunction + endmethods +endclassdef diff --git a/samples/Octave/anova1.m b/samples/Octave/anova1.m new file mode 100644 index 0000000000..fe799c9339 --- /dev/null +++ b/samples/Octave/anova1.m @@ -0,0 +1,351 @@ +## Copyright (C) 2021-2022 Andreas Bertsatos +## Copyright (C) 2022 Andrew Penn +## +## This file is part of the statistics package for GNU Octave. +## +## 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 . + +## -*- texinfo -*- +## @deftypefn {statistics} {@var{p} =} anova1 (@var{x}) +## @deftypefnx {statistics} {@var{p} =} anova1 (@var{x}, @var{group}) +## @deftypefnx {statistics} {@var{p} =} anova1 (@var{x}, @var{group}, @var{displayopt}) +## @deftypefnx {statistics} {@var{p} =} anova1 (@var{x}, @var{group}, @var{displayopt}, @var{vartype}) +## @deftypefnx {statistics} {[@var{p}, @var{atab}] =} anova1 (@var{x}, @dots{}) +## @deftypefnx {statistics} {[@var{p}, @var{atab}, @var{stats}] =} anova1 (@var{x}, @dots{}) +## +## Perform a one-way analysis of variance (ANOVA) for comparing the means of two +## or more groups of data under the null hypothesis that the groups are drawn +## from distributions with the same mean. For planned contrasts and/or +## diagnostic plots, use @qcode{anovan} instead. +## +## anova1 can take up to three input arguments: +## +## @itemize +## @item +## @var{x} contains the data and it can either be a vector or matrix. +## If @var{x} is a matrix, then each column is treated as a separate group. +## If @var{x} is a vector, then the @var{group} argument is mandatory. +## +## @item +## @var{group} contains the names for each group. If @var{x} is a matrix, then +## @var{group} can either be a cell array of strings of a character array, with +## one row per column of @var{x}. If you want to omit this argument, enter an +## empty array ([]). If @var{x} is a vector, then @var{group} must be a vector +## of the same length, or a string array or cell array of strings with one row +## for each element of @var{x}. @var{x} values corresponding to the same value +## of @var{group} are placed in the same group. +## +## @item +## @var{displayopt} is an optional parameter for displaying the groups contained +## in the data in a boxplot. If omitted, it is 'on' by default. If group names +## are defined in @var{group}, these are used to identify the groups in the +## boxplot. Use 'off' to omit displaying this figure. +## +## @item +## @var{vartype} is an optional parameter to used to indicate whether the +## groups can be assumed to come from populations with equal variance. When +## @qcode{vartype} is @qcode{"equal"} the variances are assumed to be equal +## (this is the default). When @qcode{vartype} is @qcode{"unequal"} the +## population variances are not assumed to be equal and Welch's ANOVA test is +## used instead. +## @end itemize +## +## anova1 can return up to three output arguments: +## +## @itemize +## @item +## @var{p} is the p-value of the null hypothesis that all group means are equal. +## +## @item +## @var{atab} is a cell array containing the results in a standard ANOVA table. +## +## @item +## @var{stats} is a structure containing statistics useful for performing +## a multiple comparison of means with the MULTCOMPARE function. +## @end itemize +## +## If anova1 is called without any output arguments, then it prints the results +## in a one-way ANOVA table to the standard output. It is also printed when +## @var{displayopt} is 'on'. +## +## +## Examples: +## +## @example +## x = meshgrid (1:6); +## x = x + normrnd (0, 1, 6, 6); +## anova1 (x, [], 'off'); +## [p, atab] = anova1(x); +## @end example +## +## +## @example +## x = ones (50, 4) .* [-2, 0, 1, 5]; +## x = x + normrnd (0, 2, 50, 4); +## groups = @{"A", "B", "C", "D"@}; +## anova1 (x, groups); +## @end example +## +## @seealso{anova2, anovan, multcompare} +## @end deftypefn + +function [p, anovatab, stats] = anova1 (x, group, displayopt, vartype) + + ## Check for valid number of input arguments + if (nargin < 1 || nargin > 4) + error ("anova1: invalid number of input arguments."); + endif + ## Add defaults + if (nargin < 2) + group = []; + endif + if (nargin < 3) + displayopt = "on"; + endif + if (nargin < 4) + vartype = "equal"; + endif + plotdata = ! (strcmp (displayopt, 'off')); + + ## Convert group to cell array from character array, make it a column + if (! isempty (group) && ischar (group)) + group = cellstr (group); + endif + if (size (group, 1) == 1) + group = group'; + endif + + ## If x is a matrix, convert it to column vector and create a + ## corresponging column vector for groups + if (length (x) < prod (size (x))) + [n, m] = size (x); + x = x(:); + gi = reshape (repmat ((1:m), n, 1), n*m, 1); + if (length (group) == 0) ## no group names are provided + group = gi; + elseif (size (group, 1) == m) ## group names exist and match columns + group = group(gi,:); + else + error ("anova1: columns in X and GROUP length do not match."); + endif + endif + + ## Check that x and group are the same size + if (! all (numel (x) == numel (group))) + error ("anova1: GROUP must be a vector with the same number of rows as x."); + endif + + ## Identify NaN values (if any) and remove them from X along with + ## their corresponding values from group vector + nonan = ! isnan (x); + x = x(nonan); + group = group(nonan, :); + + ## Convert group to indices and separate names + [group_id, group_names] = grp2idx (group); + group_id = group_id(:); + named = 1; + + ## Center data to improve accuracy and keep uncentered data for ploting + xorig = x; + mu = mean(x); + x = x - mu; + xr = x; + + ## Get group size and mean for each group + groups = size (group_names, 1); + xs = zeros (1, groups); + xm = xs; + xv = xs; + for j = 1:groups + group_size = find (group_id == j); + xs(j) = length (group_size); + xm(j) = mean (xr(group_size)); + xv(j) = var (xr(group_size), 0); + endfor + + ## Calculate statistics + lx = length (xr); ## Number of samples in groups + gm = mean (xr); ## Grand mean of groups + dfm = length (xm) - 1; ## degrees of freedom for model + dfe = lx - dfm - 1; ## degrees of freedom for error + SSM = xs .* (xm - gm) * (xm - gm)'; ## Sum of Squares for Model + SST = (xr(:) - gm)' * (xr(:) - gm); ## Sum of Squares Total + SSE = SST - SSM; ## Sum of Squares Error + if (dfm > 0) + MSM = SSM / dfm; ## Mean Square for Model + else + MSM = NaN; + endif + if (dfe > 0) + MSE = SSE / dfe; ## Mean Square for Error + else + MSE = NaN; + endif + ## Calculate F statistic + if (SSE != 0) ## Regular Matrix case. + switch (lower (vartype)) + case "equal" + ## Assume equal variances (Fisher's One-way ANOVA) + F = (SSM / dfm) / MSE; + case "unequal" + ## Accomodate for unequal variances (Welch's One-way ANOVA) + ## Calculate the sampling variance for each group (i.e. the square of the SEM) + sv = xv ./ xs; + ## Calculate weights as the reciprocal of the sampling variance + w = 1 ./ sv; + ## Calculate the origin + ori = sum (w .* xm) ./ sum (w); + ## Calculate Welch's F statistic + F = (groups - 1)^-1 * sum (w .* (xm - ori).^2) /... + (1 + ((2 * (groups - 2)/(groups^2 - 1)) * ... + sum ((1 - w / sum (w)).^2 .* (xs - 1).^-1))); + ## Welch's test does not use a pooled error term + MSE = NaN; + ## Correct the error degrees of freedom + dfe = (3 /(groups^2 - 1) * sum ((1 - w / sum (w)).^2 .* (xs-1).^-1))^-1; + otherwise + error ("anova1: invalid fourth (vartype) argument to anova1."); + endswitch + p = 1 - fcdf (F, dfm, dfe); ## Probability of F given equal means. + elseif (SSM == 0) ## Constant Matrix case. + F = 0; + p = 1; + else ## Perfect fit case. + F = Inf; + p = 0; + end + + ## Create results table (if requested) + if (nargout > 1) + switch (lower (vartype)) + case "equal" + anovatab = {"Source", "SS", "df", "MS", "F", "Prob>F"; ... + "Groups", SSM, dfm, MSM, F, p; ... + "Error", SSE, dfe, MSE, "", ""; ... + "Total", SST, dfm + dfe, "", "", ""}; + case "unequal" + anovatab = {"Source", "F", "df", "dfe", "F", "Prob>F"; ... + "Groups", SSM, dfm, dfe, F, p}; + endswitch + endif + ## Create stats structure (if requested) for MULTCOMPARE + if (nargout > 2) + if (length (group_names) > 0) + stats.gnames = group_names; + else + stats.gnames = strjust (num2str ((1:length (xm))'), 'left'); + end + stats.n = xs; + stats.source = 'anova1'; + stats.vartype = vartype; + stats.means = xm + mu; + stats.vars = xv; + stats.df = dfe; + stats.s = sqrt (MSE); + endif + ## Print results table on screen if no output argument was requested + if (nargout == 0 || plotdata) + switch (lower (vartype)) + case "equal" + printf("\n ANOVA Table\n\n"); + printf("Source SS df MS F Prob>F\n"); + printf("------------------------------------------------------\n"); + printf("Groups %10.4f %5.0f %10.4f %8.2f %9.4f\n", SSM, dfm, MSM, F, p); + printf("Error %10.4f %5.0f %10.4f\n", SSE, dfe, MSE); + printf("Total %10.4f %5.0f\n\n", SST, dfm + dfe); + case "unequal" + printf("\n Welch's ANOVA Table\n\n"); + printf("Source F df dfe Prob>F\n"); + printf("-----------------------------------------\n"); + printf("Groups %8.2f %5.0f %7.2f %10.4f\n\n", F, dfm, dfe, p); + endswitch + endif + ## Plot data using BOXPLOT (unless opted out) + if (plotdata) + boxplot (x, group_id, "Notch", "on", "Labels", group_names); + endif +endfunction + + +%!demo +%! x = meshgrid (1:6); +%! randn ("seed", 15); # for reproducibility +%! x = x + normrnd (0, 1, 6, 6); +%! anova1 (x, [], 'off'); + +%!demo +%! x = meshgrid (1:6); +%! randn ("seed", 15); # for reproducibility +%! x = x + normrnd (0, 1, 6, 6); +%! [p, atab] = anova1(x); + +%!demo +%! x = ones (50, 4) .* [-2, 0, 1, 5]; +%! randn ("seed", 13); # for reproducibility +%! x = x + normrnd (0, 2, 50, 4); +%! groups = {"A", "B", "C", "D"}; +%! anova1 (x, groups); + +%!demo +%! y = [54 87 45; 23 98 39; 45 64 51; 54 77 49; 45 89 50; 47 NaN 55]; +%! g = [1 2 3 ; 1 2 3 ; 1 2 3 ; 1 2 3 ; 1 2 3 ; 1 2 3 ]; +%! anova1 (y(:), g(:), "on", "unequal"); + +## testing against GEAR.DAT data file and results for one-factor ANOVA from +## https://www.itl.nist.gov/div898/handbook/eda/section3/eda354.htm +%!test +%! data = [1.006, 0.996, 0.998, 1.000, 0.992, 0.993, 1.002, 0.999, 0.994, 1.000, ... +%! 0.998, 1.006, 1.000, 1.002, 0.997, 0.998, 0.996, 1.000, 1.006, 0.988, ... +%! 0.991, 0.987, 0.997, 0.999, 0.995, 0.994, 1.000, 0.999, 0.996, 0.996, ... +%! 1.005, 1.002, 0.994, 1.000, 0.995, 0.994, 0.998, 0.996, 1.002, 0.996, ... +%! 0.998, 0.998, 0.982, 0.990, 1.002, 0.984, 0.996, 0.993, 0.980, 0.996, ... +%! 1.009, 1.013, 1.009, 0.997, 0.988, 1.002, 0.995, 0.998, 0.981, 0.996, ... +%! 0.990, 1.004, 0.996, 1.001, 0.998, 1.000, 1.018, 1.010, 0.996, 1.002, ... +%! 0.998, 1.000, 1.006, 1.000, 1.002, 0.996, 0.998, 0.996, 1.002, 1.006, ... +%! 1.002, 0.998, 0.996, 0.995, 0.996, 1.004, 1.004, 0.998, 0.999, 0.991, ... +%! 0.991, 0.995, 0.984, 0.994, 0.997, 0.997, 0.991, 0.998, 1.004, 0.997]; +%! group = [1:10] .* ones (10,10); +%! group = group(:); +%! [p, tbl] = anova1 (data, group, "off"); +%! assert (p, 0.022661, 1e-6); +%! assert (tbl{2,5}, 2.2969, 1e-4); +%! assert (tbl{2,3}, 9, 0); +%! assert (tbl{4,2}, 0.003903, 1e-6); +%! data = reshape (data, 10, 10); +%! [p, tbl, stats] = anova1 (data, [], "off"); +%! assert (p, 0.022661, 1e-6); +%! assert (tbl{2,5}, 2.2969, 1e-4); +%! assert (tbl{2,3}, 9, 0); +%! assert (tbl{4,2}, 0.003903, 1e-6); +%! means = [0.998, 0.9991, 0.9954, 0.9982, 0.9919, 0.9988, 1.0015, 1.0004, 0.9983, 0.9948]; +%! N = 10 * ones (1, 10); +%! assert (stats.means, means, 1e-6); +%! assert (length (stats.gnames), 10, 0); +%! assert (stats.n, N, 0); + +## testing against one-way ANOVA example dataset from GraphPad Prism 8 +%!test +%! y = [54 87 45; 23 98 39; 45 64 51; 54 77 49; 45 89 50; 47 NaN 55]; +%! g = [1 2 3 ; 1 2 3 ; 1 2 3 ; 1 2 3 ; 1 2 3 ; 1 2 3 ]; +%! [p, tbl] = anova1 (y(:), g(:), "off", "equal"); +%! assert (p, 0.00004163, 1e-6); +%! assert (tbl{2,5}, 22.573418, 1e-6); +%! assert (tbl{2,3}, 2, 0); +%! assert (tbl{3,3}, 14, 0); +%! [p, tbl] = anova1 (y(:), g(:), "off", "unequal"); +%! assert (p, 0.00208877, 1e-8); +%! assert (tbl{2,5}, 15.523192, 1e-6); +%! assert (tbl{2,3}, 2, 0); +%! assert (tbl{2,4}, 7.5786897, 1e-6); diff --git a/samples/Octave/silhouette.m b/samples/Octave/silhouette.m new file mode 100644 index 0000000000..e3c15c4713 --- /dev/null +++ b/samples/Octave/silhouette.m @@ -0,0 +1,206 @@ +## Copyright (C) 2016 Nan Zhou +## Copyright (C) 2021 Stefano Guidoni +## +## This file is part of the statistics package for GNU Octave. +## +## 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 . + +## -*- texinfo -*- +## @deftypefn {statistics} {} silhouette (@var{X}, @var{clust}) +## @deftypefnx {statistics} {[@var{si}, @var{h}] =} silhouette (@var{X}, @var{clust}) +## @deftypefnx {statistics} {[@var{si}, @var{h}] =} silhouette (@dots{}, @var{Metric}, @var{MetricArg}) +## +## Compute the silhouette values of clustered data and show them on a plot. +## +## @var{X} is a n-by-p matrix of n data points in a p-dimensional space. Each +## datapoint is assigned to a cluster using @var{clust}, a vector of n elements, +## one cluster assignment for each data point. +## +## Each silhouette value of @var{si}, a vector of size n, is a measure of the +## likelihood that a data point is accurately classified to the right cluster. +## Defining "a" as the mean distance between a point and the other points from +## its cluster, and "b" as the mean distance between that point and the points +## from other clusters, the silhouette value of the i-th point is: +## +## @tex +## \def\frac#1#2{{\begingroup#1\endgroup\over#2}} +## $$ S_i = \frac{b_i - a_i}{max(a_1,b_i)} $$ +## @end tex +## @ifnottex +## @verbatim +## bi - ai +## Si = ------------ +## max(ai,bi) +## @end verbatim +## @end ifnottex +## +## Each element of @var{si} ranges from -1, minimum likelihood of a correct +## classification, to 1, maximum likelihood. +## +## Optional input value @var{Metric} is the metric used to compute the distances +## between data points. Since @code{silhouette} uses @code{pdist} to compute +## these distances, @var{Metric} is quite similar to the option @var{Metric} of +## pdist and it can be: +## @itemize @bullet +## @item A known distance metric defined as a string: @qcode{Euclidean}, +## @qcode{sqEuclidean} (default), @qcode{cityblock}, @qcode{cosine}, +## @qcode{correlation}, @qcode{Hamming}, @qcode{Jaccard}. +## +## @item A vector as those created by @code{pdist}. In this case @var{X} does +## nothing. +## +## @item A function handle that is passed to @code{pdist} with @var{MetricArg} +## as optional inputs. +## @end itemize +## +## Optional return value @var{h} is a handle to the silhouette plot. +## +## @strong{Reference} +## Peter J. Rousseeuw, Silhouettes: a Graphical Aid to the Interpretation and +## Validation of Cluster Analysis. 1987. doi:10.1016/0377-0427(87)90125-7 +## @end deftypefn +## +## @seealso{dendrogram, evalclusters, kmeans, linkage, pdist} + +function [si, h] = silhouette (X, clust, metric = "sqeuclidean", varargin) + ## check the input parameters + if (nargin < 2) + print_usage (); + endif + + n = size (clust, 1); + + ## check size + if (! isempty (X)) + if (size (X, 1) != n) + error ("First dimension of X <%d> doesn't match that of clust <%d>",... + size (X, 1), n); + endif + endif + + ## check metric + if (ischar (metric)) + metric = lower (metric); + switch (metric) + case "sqeuclidean" + metric = "squaredeuclidean"; + case { "euclidean", "cityblock", "cosine", ... + "correlation", "hamming", "jaccard" } + ; + otherwise + error ("silhouette: invalid metric '%s'", metric); + endswitch + elseif (isnumeric (metric) && isvector (metric)) + ## X can be omitted when using this + distMatrix = squareform (metric); + if (size (distMatrix, 1) != n) + error ("First dimension of X <%d> doesn't match that of clust <%d>",... + size (distMatrix, 1), n); + endif + endif + + ## main + si = zeros(n, 1); + clusterIDs = unique (clust); # eg [1; 2; 3; 4] + m = length (clusterIDs); + + ## if only one cluster is defined, the silhouette value is not defined + if (m == 1) + si = NaN * ones (n, 1); + return; + endif + + ## distance matrix showing the distance for any two rows of X + if (! exist ('distMatrix', 'var')) + distMatrix = squareform (pdist (X, metric, varargin{:})); + endif + + ## calculate values of si one by one + for iii = 1 : length (si) + + ## allocate values to clusters + groupedValues = {}; + for jjj = 1 : m + groupedValues{clusterIDs(jjj)} = [distMatrix(iii, ... + clust == clusterIDs(jjj))]; + endfor + ## end allocation + + ## calculate a(i) + ## average distance of iii to all other objects in the same cluster + if (length (groupedValues{clust(iii)}) == 1) + si(iii) = 1; + continue; + else + a_i = (sum (groupedValues{clust(iii)})) / ... + (length (groupedValues{clust(iii)}) - 1); + endif + ## end a(i) + + + ## calculate b(i) + clusterIDs_new = clusterIDs; + ## remove the cluster iii in + clusterIDs_new(find (clusterIDs_new == clust(iii))) = []; + ## average distance of iii to all objects of another cluster + a_iii_2others = zeros (length (clusterIDs_new), 1); + for jjj = 1 : length (clusterIDs_new) + a_iii_2others(jjj) = mean (groupedValues{clusterIDs_new(jjj)}); + endfor + b_i = min (a_iii_2others); + ## end b(i) + + + ## calculate s(i) + si(iii) = (b_i - a_i) / (max ([a_i; b_i])); + ## end s(i) + endfor + + ## plot + ## a poor man silhouette graph + vBarsc = zeros (m, 1); + vPadding = [0; 0; 0; 0]; + Bars = vPadding; + + for i = 1 : m + vBar = si(find (clust == clusterIDs(i))); + vBarsc(i) = length (Bars) + (length (vBar) / 2); + Bars = [Bars; (sort (vBar, "descend")); vPadding]; + endfor + + figure(); + h = barh (Bars, "hist", "facecolor", [0 0.4471 0.7412]); + + xlabel ("Silhouette Value"); + ylabel ("Cluster"); + set (gca, "ytick", vBarsc, "yticklabel", clusterIDs); + ylim ([0 (length (Bars))]); + axis ("ij"); +endfunction + + +%!demo +%! load fisheriris; +%! X = meas(:,3:4); +%! cidcs = kmeans (X, 3, "Replicates", 5); +%! silhouette (X, cidcs); +%! y_labels(cidcs([1 51 101])) = unique (species); +%! set (gca, "yticklabel", y_labels); +%! title ("Fisher's iris data"); + +## Test input validation +%!error silhouette (); +%!error silhouette ([1 2; 1 1]); +%!error silhouette ([1 2; 1 1], [1 2 3]'); +%!error silhouette ([1 2; 1 1], [1 2]', "xxx"); diff --git a/samples/Octave/var.m b/samples/Octave/var.m new file mode 100644 index 0000000000..cde793acb0 --- /dev/null +++ b/samples/Octave/var.m @@ -0,0 +1,956 @@ +######################################################################## +## +## Copyright (C) 1995-2024 The Octave Project Developers +## +## See the file COPYRIGHT.md in the top-level directory of this +## distribution or . +## +## This file is part of Octave. +## +## 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 3 of the License, or +## (at your option) any later version. +## +## Octave 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 +## . +## +######################################################################## + +## -*- texinfo -*- +## @deftypefn {} {@var{v} =} var (@var{x}) +## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}) +## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}, @var{dim}) +## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}, @var{vecdim}) +## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}, @qcode{"all"}) +## @deftypefnx {} {@var{v} =} var (@dots{}, @var{nanflag}) +## @deftypefnx {} {[@var{v}, @var{m}] =} var (@dots{}) +## Compute the variance of the elements of the vector @var{x}. +## +## The variance is defined as +## @tex +## $$ {\rm var}(x) = {1\over N-1} \sum_{i=1}^N (x_i - \bar x )^2 $$ +## where $\bar{x}$ is the mean value of @var{x} and $N$ is the number of +## elements of @var{x}. +## @end tex +## @ifnottex +## +## @example +## @group +## var (@var{x}) = (1 / (N-1)) * SUM_i ((@var{x}(i) - mean(@var{x}))^2) +## @end group +## @end example +## +## @noindent +## where @math{N} is the number of elements of @var{x}. +## @end ifnottex +## +## If @var{x} is an array, compute the variance along the first non-singleton +## dimensions of @var{x}. +## +## The optional argument @var{w} determines the weighting scheme to use. Valid +## values are: +## +## @table @asis +## @item 0 [default]: +## Normalize with @math{N-1} (population variance). This provides the square +## root of the best unbiased estimator of the variance. +## +## @item 1: +## Normalize with @math{N} (sample variance). This provides the square root of +## the second moment around the mean. +## +## @item a vector: +## Compute the weighted variance with non-negative weights. The length of +## @var{w} must equal the size of @var{x} in the operating dimension. NaN +## values are permitted in @var{w}, will be multiplied with the associated +## values in @var{x}, and can be excluded by the @var{nanflag} option. +## +## @item an array: +## Similar to vector weights, but @var{w} must be the same size as @var{x}. If +## the operating dimension is supplied as @var{vecdim} or @qcode{"all"} and +## @var{w} is not a scalar, @var{w} must be an same-sized array. +## @end table +## +## Note: @var{w} must always be specified before specifying any of the +## following dimension options. To use the default value for @var{w} you +## may pass an empty input argument []. +## +## The optional variable @var{dim} forces @code{var} to operate over the +## specified dimension, which must be a positive integer-valued number. +## Specifying any singleton dimension in @var{x}, including any dimension +## exceeding @code{ndims (@var{x})}, will result in a variance of 0. +## +## Specifying the dimensions as @var{vecdim}, a vector of non-repeating +## dimensions, will return the variance calculated over the array slice defined +## by @var{vecdim}. If @var{vecdim} indexes all dimensions of @var{x}, then it +## is equivalent to the option @qcode{"all"}. Any dimension in @var{vecdim} +## greater than @code{ndims (@var{x})} is ignored. +## +## Specifying the dimension as @qcode{"all"} will force @code{var} to +## operate on all elements of @var{x}, and is equivalent to @code{var +## (@var{x}(:))}. +## +## The optional variable @var{nanflag} specifies whether to include or exclude +## NaN values from the calculation using any of the previously specified input +## argument combinations. The default value for @var{nanflag} is +## @qcode{"includenan"} which keeps NaN values in the calculation. To +## exclude NaN values set the value of @var{nanflag} to @qcode{"omitnan"}. +## The output will still contain NaN values if @var{x} consists of all NaN +## values in the operating dimension. +## +## The optional second output variable @var{m} contains the mean of the +## elements of @var{x} used to calculate the variance. If @var{v} is the +## weighted variance, then @var{m} is also the weighted mean. +## +## @seealso{std, mean, cov, skewness, kurtosis, moment} +## @end deftypefn + +function [v, m] = var (x, varargin) + + if (nargin < 1 || nargin > 4) + print_usage (); + endif + + ## initialize variables + all_flag = false; + omitnan = false; + nvarg = numel (varargin); + varg_chars = cellfun ('ischar', varargin); + + ## Check all char arguments. + if (nvarg == 3 && ! varg_chars(3)) + print_usage (); + endif + + if (any (varg_chars)) + for argin = varargin(varg_chars) + switch (lower (argin{1})) + case "all" + all_flag = true; + case "omitnan" + omitnan = true; + case "includenan" + omitnan = false; + otherwise + print_usage (); + endswitch + endfor + varargin(varg_chars) = []; + nvarg = numel (varargin); + endif + + ## FIXME: When sparse can broadcast ops then remove sparse checks and hacks. + x_issparse = issparse (x); + w = 0; + weighted = false; # true if weight vector/array used + vecdim = []; + vecempty = true; + vecdim_scalar_vector = [false, false]; # [false, false] for empty vecdim + szx = size (x); + ndx = ndims (x); + + ## Check numeric arguments + if (! (isnumeric (x))) + error ("var: X must be a numeric vector or matrix"); + endif + if (isa (x, "single")) + outtype = "single"; + else + outtype = "double"; + endif + + if (nvarg > 0) + if (nvarg > 2 || any (! cellfun ('isnumeric', varargin))) + print_usage (); + endif + ## Process weight input + if (any (varargin{1} < 0)) + error ("var: weights must not contain any negative values"); + endif + if (isscalar (varargin{1})) + w = varargin{1}; + if (! (w == 0 || w == 1) && ! isscalar (x)) + error ("var: normalization scalar must be either 0 or 1"); + endif + elseif (numel (varargin{1}) > 1) + weights = varargin{1}; + weighted = true; + endif + if (nvarg > 1) + ## Process dimension input + vecdim = varargin{2}; + if (! (vecempty = isempty (vecdim))) + ## Check for empty vecdim, won't change vsv if nonzero size empty + vecdim_scalar_vector = [isscalar(vecdim), isvector(vecdim)]; + endif + if (! (vecdim_scalar_vector(2) && all (vecdim > 0)) + || any (rem (vecdim, 1))) + error ("var: DIM must be a positive integer scalar or vector"); + endif + if (vecdim_scalar_vector(1) && vecdim > ndx && ! isempty (x)) + ## Scalar dimension larger than ndims(x), variance of any single number + ## is zero, except for inf, NaN, and empty values of x. + v = zeros (szx, outtype); + vn = ! isfinite (x); + v(vn) = NaN; + m = x; + return; + endif + if (vecdim_scalar_vector == [0 1] && (! all (diff (sort (vecdim))))) + error ("var: VECDIM must contain non-repeating positive integers"); + endif + endif + endif + + ## Check for conflicting input arguments + if (all_flag && ! vecempty) + error ("var: 'all' flag cannot be used with DIM or VECDIM options"); + endif + if (weighted) + if (all_flag) + if (isvector (weights)) + if (numel (weights) != numel (x)) + error ("var: weight vector element count does not match X"); + endif + elseif (! (isequal (size (weights), szx))) + error ("var: weight matrix or array does not match X in size"); + endif + + elseif (vecempty) + ## Find the first non-singleton dimension. + (dim = find (szx > 1, 1)) || (dim = 1); + if (isvector (weights)) + if (numel (weights) != szx(dim)) + error ("var: weight vector length does not match operating dimension"); + endif + elseif (! isequal (size (weights), szx)) + error ("var: weight matrix or array does not match X in size"); + endif + elseif (vecdim_scalar_vector(1)) + if (isvector (weights)) + if (numel (weights) != szx(vecdim)) + error ("var: weight vector length does not match operating dimension"); + endif + elseif (! isequal (size (weights), szx)) + error ("var: weight matrix or array does not match X in size"); + endif + + elseif (vecdim_scalar_vector(2) && ! (isequal (size (weights), szx))) + error ("var: weight matrix or array does not match X in size"); + endif + endif + + ## Handle special cases of empty or scalar X and exit early. + if (isempty (x)) + if (vecempty && (ndx == 2 || all (szx == 0))) + v = NaN (outtype); + if (nargout > 1) + m = NaN (outtype); + endif + return; + endif + if (vecdim_scalar_vector(1)) + szx(vecdim) = 1; + v = NaN (szx, outtype); + if (nargout > 1) + m = NaN (szx, outtype); + endif + return; + endif + elseif (isscalar (x)) + if (isfinite (x)) + v = zeros (outtype); + else + v = NaN (outtype); + endif + if (nargout > 1) + m = x; + endif + return; + endif + + if (nvarg == 0) + ## Only numeric input argument, no dimensions or weights. + if (all_flag) + x = x(:); + if (omitnan) + x = x(! isnan (x)); + endif + n = length (x); + m = sum (x) ./ n; + v = sum (abs (x - m) .^ 2) ./ (n - 1 + w); + if (n == 1) + v = 0; + endif + else + ## Find the first non-singleton dimension. + (dim = find (szx > 1, 1)) || (dim = 1); + n = szx(dim); + if (omitnan) + n = sum (! isnan (x), dim); + xn = isnan (x); + x(xn) = 0; + endif + m = sum (x, dim) ./ n; + dims = ones (1, ndx); + dims(dim) = szx(dim); + if (x_issparse) + m_exp = repmat (m, dims); + else + m_exp = m .* ones (dims); + endif + if (omitnan) + x(xn) = m_exp(xn); + endif + v = sumsq (x - m_exp, dim) ./ (n - 1 + w); + if (numel (n) == 1) + divby0 = (n .* ones (size (v))) == 1; + else + divby0 = n == 1; + endif + v(divby0) = 0; + endif + + elseif (nvarg == 1) + ## Two numeric input arguments, w or weights given. + if (all_flag) + x = x(:); + if (weighted) + weights = weights(:); + wx = weights .* x; + else + weights = ones (length (x), 1); + wx = x; + endif + + if (omitnan) + xn = isnan (wx); + wx = wx(! xn); + weights = weights(! xn); + x = x(! xn); + endif + n = length (wx); + m = sum (wx) ./ sum (weights); + if (weighted) + v = sum (weights .* (abs (x - m) .^ 2)) ./ sum (weights); + else + v = sum (weights .* (abs (x - m) .^ 2)) ./ (n - 1 + w); + if (n == 1) + v = 0; + endif + endif + + else + ## Find the first non-singleton dimension. + (dim = find (szx > 1, 1)) || (dim = 1); + if (! weighted) + weights = ones (szx); + wx = x; + else + if (isvector (weights)) + dims = 1:ndx; + dims([1, dim]) = [dim, 1]; + weights = zeros (szx) + permute (weights(:), dims); + endif + wx = weights .* x; + endif + n = size (wx, dim); + if (omitnan) + xn = isnan (wx); + n = sum (! xn, dim); + wx(xn) = 0; + weights(xn) = 0; + endif + m = sum (wx, dim) ./ sum (weights, dim); + dims = ones (1, ndims (wx)); + dims(dim) = size (wx, dim); + if (x_issparse) + m_exp = repmat (m, dims); + else + m_exp = m .* ones (dims); + endif + if (omitnan) + x(xn) = m_exp(xn); + endif + if (weighted) + v = sum (weights .* ((x - m_exp) .^ 2), dim) ./ sum (weights, dim); + else + v = sumsq (x - m_exp, dim) ./ (n - 1 + w); + if (numel (n) == 1) + divby0 = (n .* ones (size (v))) == 1; + else + divby0 = n == 1; + endif + v(divby0) = 0; + endif + endif + + elseif (nvarg == 2) + ## Three numeric input arguments, both w or weights and dim or vecdim given. + if (vecdim_scalar_vector(1)) + if (! weighted) + weights = ones (szx); + wx = x; + else + if (isvector (weights)) + dims = 1:ndx; + dims([1, vecdim]) = [vecdim, 1]; + weights = zeros (szx) + permute (weights(:), dims); + endif + wx = weights .* x; + endif + n = size (wx, vecdim); + if (omitnan) + n = sum (! isnan (wx), vecdim); + xn = isnan (wx); + wx(xn) = 0; + weights(xn) = 0; + endif + m = sum (wx, vecdim) ./ sum (weights, vecdim); + dims = ones (1, ndims (wx)); + dims(vecdim) = size (wx, vecdim); + if (x_issparse) + m_exp = repmat (m, dims); + else + m_exp = m .* ones (dims); + endif + if (omitnan) + x(xn) = m_exp(xn); + endif + if (weighted) + v = sum (weights .* ((x - m_exp) .^ 2), vecdim) ... + ./ sum (weights, vecdim); + else + v = sumsq (x - m_exp, vecdim); + vn = isnan (v); + v = v ./ (n - 1 + w); + if (numel (n) == 1) + divby0 = (n .* ones (size (v))) == 1; + else + divby0 = n == 1; + endif + v(divby0) = 0; + v(vn) = NaN; + endif + + else + ## Weights and nonscalar vecdim specified + + ## Ignore dimensions in VECDIM larger than actual array. + remdims = 1 : ndx; # All dimensions + vecdim(vecdim > ndx) = []; + ## Calculate permutation vector + remdims(vecdim) = []; # Delete dimensions specified by vecdim + nremd = numel (remdims); + + ## If all dimensions are given, it is equivalent to the 'all' flag. + if (nremd == 0) + x = x(:); + if (weighted) + weights = weights(:); + wx = weights .* x; + else + weights = ones (length (x), 1); + wx = x; + endif + + if (omitnan) + xn = isnan (wx); + wx = wx(! xn); + weights = weights(! xn); + x = x(! xn); + endif + n = length (wx); + m = sum (wx) ./ sum (weights); + if (weighted) + v = sum (weights .* (abs (x - m) .^ 2)) ./ sum (weights); + else + v = sum (weights .* (abs (x - m) .^ 2)) ./ (n - 1 + w); + if (n == 1) + v = 0; + endif + endif + + else + + ## FIXME: much of the reshaping can be skipped once Octave's sum can + ## take a vecdim argument. + + ## Apply weights + if (weighted) + wx = weights .* x; + else + weights = ones (szx); + wx = x; + endif + + ## Permute to push vecdims to back + perm = [remdims, vecdim]; + wx = permute (wx, perm); + weights = permute (weights, perm); + x = permute (x, perm); + + ## Reshape to squash all vecdims in final dimension + szwx = size (wx); + sznew = [szwx(1:nremd), prod(szwx(nremd+1:end))]; + wx = reshape (wx, sznew); + weights = reshape (weights, sznew); + x = reshape (x, sznew); + + ## Calculate var on final squashed dimension + dim = nremd + 1; + n = size (wx, dim); + if (omitnan) + xn = isnan (wx); + n = sum (! xn, dim); + wx(xn) = 0; + weights(xn) = 0; + endif + m = sum (wx, dim) ./ sum (weights, dim); + m_exp = zeros (sznew) + m; + if (omitnan) + x(xn) = m_exp(xn); + endif + if (weighted) + v = sum (weights .* ((x - m_exp) .^ 2), dim) ./ sum (weights, dim); + else + v = sumsq (x - m_exp, dim) ./ (n - 1 + w); + if (numel (n) == 1) + divby0 = n .* ones (size (v)) == 1; + else + divby0 = n == 1; + endif + v(divby0) = 0; + endif + + ## Inverse permute back to correct dimensions + v = ipermute (v, perm); + if (nargout > 1) + m = ipermute (m, perm); + endif + endif + endif + endif + + ## Preserve class type + if (nargout < 2) + if (strcmp (outtype, "single")) + v = single (v); + else + v = double (v); + endif + else + if (strcmp (outtype, "single")) + v = single (v); + m = single (m); + else + v = double (v); + m = double (m); + endif + endif + +endfunction + + +%!assert (var (13), 0) +%!assert (var (single (13)), single (0)) +%!assert (var ([1,2,3]), 1) +%!assert (var ([1,2,3], 1), 2/3, eps) +%!assert (var ([1,2,3], [], 1), [0,0,0]) +%!assert (var ([1,2,3], [], 3), [0,0,0]) +%!assert (var (5, 99), 0) +%!assert (var (5, 99, 1), 0) +%!assert (var (5, 99, 2), 0) +%!assert (var ([5 3], [99 99], 2), 1) +%!assert (var ([1:7], [1:7]), 3) +%!assert (var ([eye(3)], [1:3]), [5/36, 2/9, 1/4], eps) +%!assert (var (ones (2,2,2), [1:2], 3), [(zeros (2,2))]) +%!assert (var ([1 2; 3 4], 0, 'all'), var ([1:4])) +%!assert (var (reshape ([1:8], 2, 2, 2), 0, [1 3]), [17/3 17/3], eps) +%!assert (var ([1 2 3;1 2 3], [], [1 2]), 0.8, eps) + +## Test single input and optional arguments "all", DIM, "omitnan") +%!test +%! x = [-10:10]; +%! y = [x;x+5;x-5]; +%! assert (var (x), 38.5); +%! assert (var (y, [], 2), [38.5; 38.5; 38.5]); +%! assert (var (y, 0, 2), [38.5; 38.5; 38.5]); +%! assert (var (y, 1, 2), ones (3,1) * 36.66666666666666, 1e-14); +%! assert (var (y, "all"), 54.19354838709678, 1e-14); +%! y(2,4) = NaN; +%! assert (var (y, "all"), NaN); +%! assert (var (y, "all", "includenan"), NaN); +%! assert (var (y, "all", "omitnan"), 55.01533580116342, 1e-14); +%! assert (var (y, 0, 2, "includenan"), [38.5; NaN; 38.5]); +%! assert (var (y, [], 2), [38.5; NaN; 38.5]); +%! assert (var (y, [], 2, "omitnan"), [38.5; 37.81842105263158; 38.5], 1e-14); + +## Tests for different weight and omitnan code paths +%!assert (var ([1 NaN 3], [1 2 3], "omitnan"), 0.75, eps) +%!assert (var ([1 2 3], [1 NaN 3], "omitnan"), 0.75, eps) +%!assert (var (magic(3), [1 NaN 3], "omitnan"), [3 12 3], eps) +%!assert (var ([1 NaN 3], [1 2 3], "omitnan", "all"), 0.75, eps) +%!assert (var ([1 NaN 3], [1 2 3], "all", "omitnan"), 0.75, eps) +%!assert (var ([1 2 3], [1 NaN 3], "omitnan", "all"), 0.75, eps) +%!assert (var ([1 NaN 3], [1 2 3], 2, "omitnan"), 0.75, eps) +%!assert (var ([1 2 3], [1 NaN 3], 2, "omitnan"), 0.75, eps) +%!assert (var (magic(3), [1 NaN 3], 1, "omitnan"), [3 12 3], eps) +%!assert (var (magic(3), [1 NaN 3], 2, "omitnan"), [0.75;3;0.75], eps) +%!assert (var ([4 4; 4 6; 6 6], [1 3], 2, 'omitnan'), [0;0.75;0], eps) +%!assert (var ([4 NaN; 4 6; 6 6], [1 2 3], 1, 'omitnan'), [1 0]) +%!assert (var ([4 NaN; 4 6; 6 6], [1 3], 2, 'omitnan'), [0;0.75;0], eps) +%!assert (var (3*reshape(1:18, [3 3 2]), [1 2 3], 1, 'omitnan'), ones(1,3,2)*5) +%!assert (var (reshape(1:18, [3 3 2]), [1 2 3], 2, 'omitnan'), 5*ones(3,1,2)) +%!assert (var (3*reshape(1:18, [3 3 2]), ones (3,3,2), [1 2], 'omitnan'), ... +%! 60 * ones(1,1,2)) +%!assert (var (3*reshape(1:18, [3 3 2]), ones (3,3,2), [1 4], 'omitnan'), ... +%! 6 * ones(1,3,2)) +%!assert (var (6*reshape(1:18, [3 3 2]), ones (3,3,2), [1:3], 'omitnan'), 969) +%!test +%! x = reshape(1:18, [3 3 2]); +%! x([2, 14]) = NaN; +%! w = ones (3,3,2); +%! assert (var (16*x, w, [1:3], 'omitnan'), 6519); +%!test +%! x = reshape(1:18, [3 3 2]); +%! w = ones (3,3,2); +%! w([2, 14]) = NaN; +%! assert (var (16*x, w, [1:3], 'omitnan'), 6519); + +## Test input case insensitivity +%!assert (var ([1 2 3], "aLl"), 1) +%!assert (var ([1 2 3], "OmitNan"), 1) +%!assert (var ([1 2 3], "IncludeNan"), 1) + +## Test dimension indexing with vecdim in n-dimensional arrays +%!test +%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); +%! assert (size (var (x, 0, [3 2])), [10, 1, 1, 3]); +%! assert (size (var (x, 1, [1 2])), [1, 1, 6, 3]); +%! assert (size (var (x, [], [1 2 4])), [1, 1, 6]); +%! assert (size (var (x, 0, [1 4 3])), [1, 40]); +%! assert (size (var (x, [], [1 2 3 4])), [1, 1]); + +## Test matrix with vecdim, weighted, matrix weights, omitnan +%!assert (var (3*magic(3)), [63 144 63]) +%!assert (var (3*magic(3), 'omitnan'), [63 144 63]) +%!assert (var (3*magic(3), 1), [42 96 42]) +%!assert (var (3*magic(3), 1, 'omitnan'), [42 96 42]) +%!assert (var (3*magic(3), ones(1,3), 1), [42 96 42]) +%!assert (var (3*magic(3), ones(1,3), 1, 'omitnan'), [42 96 42]) +%!assert (var (2*magic(3), [1 1 NaN], 1, 'omitnan'), [25 16 1]) +%!assert (var (3*magic(3), ones(3,3)), [42 96 42]) +%!assert (var (3*magic(3), ones(3,3), 'omitnan'), [42 96 42]) +%!assert (var (3*magic(3), [1 1 1; 1 1 1; 1 NaN 1], 'omitnan'), [42 36 42]) +%!assert (var (3*magic(3), ones(3,3), 1), [42 96 42]) +%!assert (var (3*magic(3), ones(3,3), 1, 'omitnan'), [42 96 42]) +%!assert (var (3*magic(3), [1 1 1; 1 1 1; 1 NaN 1], 1, 'omitnan'), [42 36 42]) +%!assert (var (3*magic(3), ones(3,3), [1 4]), [42 96 42]) +%!assert (var (3*magic(3), ones(3,3), [1 4], 'omitnan'), [42 96 42]) +%!assert (var (3*magic(3), [1 1 1; 1 1 1; 1 NaN 1],[1 4],'omitnan'), [42 36 42]) + +## Test results with vecdim in n-dimensional arrays and "omitnan" +%!test +%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); +%! v = repmat (33.38912133891213, [10, 1, 1, 3]); +%! assert (var (x, 0, [3, 2]), v, 1e-14); +%! v = repmat (33.250, [10, 1, 1, 3]); +%! assert (var (x, 1, [3, 2]), v, 1e-14); +%! x(2,5,6,3) = NaN; +%! v(2,1,1,3) = NaN; +%! assert (var (x, 1, [3, 2]), v, 4e-14); +%! v = repmat (33.38912133891213, [10 1 1 3]); +%! v(2,1,1,3) = NaN; +%! assert (var (x, [], [3, 2]), v, 4e-14); +%! v(2,1,1,3) = 33.40177912169048; +%! assert (var (x, [], [3, 2], "omitnan"), v, 4e-14); + +## Testing weights vector & arrays +%!assert (var (ones (2,2,2), [1:2], 3), [(zeros (2, 2))]) +%!assert (var (magic (3), [1:9], "all"), 6.666666666666667, 1e-14) + +## Test exceeding dimensions +%!assert (var (ones (2,2), [], 3), zeros (2,2)) +%!assert (var (ones (2,2,2), [], 99), zeros (2,2,2)) +%!assert (var (magic (3), [], 3), zeros (3,3)) +%!assert (var (magic (3), [], 1), [7, 16, 7]) +%!assert (var (magic (3), [], [1 3]), [7, 16, 7]) +%!assert (var (magic (3), [], [1 99]), [7, 16, 7]) + +## Test empty inputs +%!assert (var ([]), NaN) +%!assert (class (var (single ([]))), "single") +%!assert (var ([],[],1), NaN(1,0)) +%!assert (var ([],[],2), NaN(0,1)) +%!assert (var ([],[],3), []) +%!assert (class (var (single ([]), [], 1)), "single") +%!assert (var (ones (1,0)), NaN) +%!assert (var (ones (1,0), [], 1), NaN(1,0)) +%!assert (var (ones (1,0), [], 2), NaN) +%!assert (var (ones (1,0), [], 3), NaN(1,0)) +%!assert (class (var (ones (1, 0, "single"), [], 1)), "single") +%!assert (var (ones (0,1)), NaN) +%!assert (var (ones (0,1), [], 1), NaN) +%!assert (var (ones (0,1), [], 2), NaN(0,1)) +%!assert (var (ones (0,1), [], 3), NaN(0,1)) +%!assert (var (ones (1,3,0,2)), NaN(1,1,0,2)) +%!assert (var (ones (1,3,0,2), [], 1), NaN(1,3,0,2)) +%!assert (var (ones (1,3,0,2), [], 2), NaN(1,1,0,2)) +%!assert (var (ones (1,3,0,2), [], 3), NaN(1,3,1,2)) +%!assert (var (ones (1,3,0,2), [], 4), NaN(1,3,0)) +%!test +%! [~, m] = var ([]); +%! assert (m, NaN); + +## Test optional mean output +%!test <*62395> +%! [~, m] = var (13); +%! assert (m, 13); +%! [~, m] = var (single(13)); +%! assert (m, single(13)); +%! [~, m] = var ([1, 2, 3; 3 2 1], []); +%! assert (m, [2 2 2]); +%! [~, m] = var ([1, 2, 3; 3 2 1], [], 1); +%! assert (m, [2 2 2]); +%! [~, m] = var ([1, 2, 3; 3 2 1], [], 2); +%! assert (m, [2 2]'); +%! [~, m] = var ([1, 2, 3; 3 2 1], [], 3); +%! assert (m, [1 2 3; 3 2 1]); + +## Test mean output, weighted inputs, vector dims +%!test <*62395> +%! [~, m] = var (5,99); +%! assert (m, 5); +%! [~, m] = var ([1:7], [1:7]); +%! assert (m, 5); +%! [~, m] = var ([eye(3)], [1:3]); +%! assert (m, [1/6, 1/3, 0.5], eps); +%! [~, m] = var (ones (2,2,2), [1:2], 3); +%! assert (m, ones (2,2)); +%! [~, m] = var ([1 2; 3 4], 0, 'all'); +%! assert (m, 2.5, eps); +%! [~, m] = var (reshape ([1:8], 2, 2, 2), 0, [1 3]); +%! assert (m, [3.5, 5.5], eps); +%!test +%! [v, m] = var (4 * eye (2), [1, 3]); +%! assert (v, [3, 3]); +%! assert (m, [1, 3]); + +## Test mean output, empty inputs, omitnan +%!test <*62395> +%! [~, m] = var ([]); +%! assert (m, NaN); +#%! [~, m] = var ([],[],1); +#%! assert (m, NaN(1,0)); +#%! [~, m] = var ([],[],2); +#%! assert (m, NaN(0,1)); +#%! [~, m] = var ([],[],3); +#%! assert (m, []); +#%! [~, m] = var (ones (1,3,0,2)); +#%! assert (m, NaN(1,1,0,2)); + +## Test mean output, nD array +%!test <*62395> +%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); +%! [~, m] = var (x, 0, [3 2]); +%! assert (m, mean (x, [3 2])); +%! [~, m] = var (x, 0, [1 2]); +%! assert (m, mean (x, [1 2])); +%! [~, m] = var (x, 0, [1 3 4]); +%! assert (m, mean (x, [1 3 4])); +%!test +%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); +%! x(2,5,6,3) = NaN; +%! [~, m] = var (x, 0, [3 2], "omitnan"); +%! assert (m, mean (x, [3 2], "omitnan")); + +## Test Inf and NaN inputs +%!test <*63203> +%! [v, m] = var (Inf); +%! assert (v, NaN); +%! assert (m, Inf); +%!test <*63203> +%! [v, m] = var (NaN); +%! assert (v, NaN); +%! assert (m, NaN); +%!test <*63203> +%! [v, m] = var ([1, Inf, 3]); +%! assert (v, NaN); +%! assert (m, Inf); +%!test <*63203> +%! [v, m] = var ([1, Inf, 3]'); +%! assert (v, NaN); +%! assert (m, Inf); +%!test <*63203> +%! [v, m] = var ([1, NaN, 3]); +%! assert (v, NaN); +%! assert (m, NaN); +%!test <*63203> +%! [v, m] = var ([1, NaN, 3]'); +%! assert (v, NaN); +%! assert (m, NaN); +%!test <*63203> +%! [v, m] = var ([1, Inf, 3], [], 1); +%! assert (v, [0, NaN, 0]); +%! assert (m, [1, Inf, 3]); +%!test <*63203> +%! [v, m] = var ([1, Inf, 3], [], 2); +%! assert (v, NaN); +%! assert (m, Inf); +%!test <*63203> +%! [v, m] = var ([1, Inf, 3], [], 3); +%! assert (v, [0, NaN, 0]); +%! assert (m, [1, Inf, 3]); +%!test <*63203> +%! [v, m] = var ([1, NaN, 3], [], 1); +%! assert (v, [0, NaN, 0]); +%! assert (m, [1, NaN, 3]); +%!test <*63203> +%! [v, m] = var ([1, NaN, 3], [], 2); +%! assert (v, NaN); +%! assert (m, NaN); +%!test <*63203> +%! [v, m] = var ([1, NaN, 3], [], 3); +%! assert (v, [0, NaN, 0]); +%! assert (m, [1, NaN, 3]); +%!test <*63203> +%! [v, m] = var ([1, 2, 3; 3, Inf, 5]); +%! assert (v, [2, NaN, 2]); +%! assert (m, [2, Inf, 4]); +%!test <*63203> +%! [v, m] = var ([1, Inf, 3; 3, Inf, 5]); +%! assert (v, [2, NaN, 2]); +%! assert (m, [2, Inf, 4]); +%!test <*63203> +%! [v, m] = var ([1, 2, 3; 3, NaN, 5]); +%! assert (v, [2, NaN, 2]); +%! assert (m, [2, NaN, 4]); +%!test <*63203> +%! [v, m] = var ([1, NaN, 3; 3, NaN, 5]); +%! assert (v, [2, NaN, 2]); +%! assert (m, [2, NaN, 4]); +%!test <*63203> +%! [v, m] = var ([Inf, 2, NaN]); +%! assert (v, NaN); +%! assert (m, NaN); +%!test <*63203> +%! [v, m] = var ([Inf, 2, NaN]'); +%! assert (v, NaN); +%! assert (m, NaN); +%!test <*63203> +%! [v, m] = var ([NaN, 2, Inf]); +%! assert (v, NaN); +%! assert (m, NaN); +%!test <*63203> +%! [v, m] = var ([NaN, 2, Inf]'); +%! assert (v, NaN); +%! assert (m, NaN); +%!test <*63203> +%! [v, m] = var ([Inf, 2, NaN], [], 1); +%! assert (v, [NaN, 0, NaN]); +%! assert (m, [Inf, 2, NaN]); +%!test <*63203> +%! [v, m] = var ([Inf, 2, NaN], [], 2); +%! assert (v, NaN); +%! assert (m, NaN); +%!test <*63203> +%! [v, m] = var ([NaN, 2, Inf], [], 1); +%! assert (v, [NaN, 0, NaN]); +%! assert (m, [NaN, 2, Inf]); +%!test <*63203> +%! [v, m] = var ([NaN, 2, Inf], [], 2); +%! assert (v, NaN); +%! assert (m, NaN); +%!test <*63203> +%! [v, m] = var ([1, 3, NaN; 3, 5, Inf]); +%! assert (v, [2, 2, NaN]); +%! assert (m, [2, 4, NaN]); +%!test <*63203> +%! [v, m] = var ([1, 3, Inf; 3, 5, NaN]); +%! assert (v, [2, 2, NaN]); +%! assert (m, [2, 4, NaN]); + +## Test sparse/diagonal inputs +%!test <*63291> +%! [v, m] = var (2 * eye (2)); +%! assert (v, [2, 2]); +%! assert (m, [1, 1]); +%!test <*63291> +%! [v, m] = var (4 * eye (2), [1, 3]); +%! assert (v, [3, 3]); +%! assert (m, [1, 3]); +%!test <*63291> +%! [v, m] = var (sparse (2 * eye (2))); +%! assert (full (v), [2, 2]); +%! assert (full (m), [1, 1]); +%!test <*63291> +%! [v, m] = var (sparse (4 * eye (2)), [1, 3]); +%! assert (full (v), [3, 3]); +%! assert (full (m), [1, 3]); +%!test<*63291> +%! [v, m] = var (sparse (eye (2))); +%! assert (issparse (v)); +%! assert (issparse (m)); +%!test<*63291> +%! [v, m] = var (sparse (eye (2)), [1, 3]); +%! assert (issparse (v)); +%! assert (issparse (m)); + +## Test input validation +%!error var () +%!error var (1, 2, "omitnan", 3) +%!error var (1, 2, 3, 4) +%!error var (1, 2, 3, 4, 5) +%!error var (1, "foo") +%!error var (1, [], "foo") +%!error var ([1 2 3], 2) +%!error var ([1 2], 2, "all") +%!error var ([1 2],0.5, "all") +%!error var (1, -1) +%!error var (1, [1 -1]) +%!error ... +%! var ([1 2 3], [1 -1 0]) +%!error var ({1:5}) +%!error var ("char") +%!error var (['A'; 'B']) +%!error var (1, [], ones (2,2)) +%!error var (1, 0, 1.5) +%!error var (1, [], 0) +%!error var (1, [], 1.5) +%!error var ([1 2 3], [], [-1 1]) +%!error ... +%! var (repmat ([1:20;6:25], [5 2 6 3]), 0, [1 2 2 2]) +%!error ... +%! var ([1 2], eye (2)) +%!error ... +%! var ([1 2 3 4], [1 2; 3 4]) +%!error ... +%! var ([1 2 3 4], [1 2; 3 4], 1) +%!error ... +%! var ([1 2 3 4], [1 2; 3 4], [2 3]) +%!error ... +%! var (ones (2, 2), [1 2], [1 2]) +%!error ... +%! var ([1 2 3 4; 5 6 7 8], [1 2 1 2 1; 1 2 1 2 1], 1) +%!error ... +%! var (repmat ([1:20;6:25], [5 2 6 3]), repmat ([1:20;6:25], [5 2 3]), [2 3]) +%!error var ([1 2 3; 2 3 4], [1 3 4]) +%!error var ([1 2], [1 2 3]) +%!error var (1, [1 2]) +%!error var ([1 2 3; 2 3 4], [1 3 4], 1) +%!error var ([1 2 3; 2 3 4], [1 3], 2) +%!error var ([1 2], [1 2], 1) +%!error <'all' flag cannot be used with DIM or VECDIM options> ... +%! var (1, [], 1, "all") +%!error ... +%! var ([1 2 3; 2 3 4], [1 3], "all") +%!error ... +%! var (repmat ([1:20;6:25], [5 2 6 3]), repmat ([1:20;6:25], [5 2 3]), "all") From 163f347f4d275f8a0b2ae75d4c142389f55e9338 Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Tue, 9 Jul 2024 16:45:27 +0300 Subject: [PATCH 05/21] Update test_heuristics.rb for Octave language --- test/test_heuristics.rb | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_heuristics.rb b/test/test_heuristics.rb index 884a039d9b..bfb138e097 100755 --- a/test/test_heuristics.rb +++ b/test/test_heuristics.rb @@ -633,6 +633,7 @@ def test_m_by_heuristics "M" => all_fixtures("M", "MDB.m"), "Mathematica" => all_fixtures("Mathematica", "*.m") - all_fixtures("Mathematica", "Problem12.m"), "MATLAB" => all_fixtures("MATLAB", "create_ieee_paper_plots.m"), + "Octave" => all_fixtures("Octave", "accumarray.m"), "Limbo" => all_fixtures("Limbo", "*.m"), nil => ambiguous }) From 873a169a6a09565cd48ea2d56ae0f872debe6a8c Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Tue, 9 Jul 2024 16:46:40 +0300 Subject: [PATCH 06/21] Update test_language.rb for Octave language --- test/test_language.rb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_language.rb b/test/test_language.rb index 0f781beeb0..287f9465b4 100644 --- a/test/test_language.rb +++ b/test/test_language.rb @@ -166,7 +166,7 @@ def test_find_all_by_alias def test_find_by_extension assert_equal [], Language.find_by_extension('.factor-rc') - assert_equal [Language['Limbo'], Language['M'], Language['MATLAB'], Language['MUF'], Language['Mathematica'], Language['Mercury'], Language['Objective-C']], Language.find_by_extension('foo.m') + assert_equal [Language['Limbo'], Language['M'], Language['MATLAB'], Language['MUF'], Language['Mathematica'], Language['Mercury'], Language['Objective-C'], Language['Octave']], Language.find_by_extension('foo.m') assert_equal [Language['Ruby']], Language.find_by_extension('foo.rb') assert_equal [Language['Ruby']], Language.find_by_extension('foo/bar.rb') assert_equal [Language['Ruby']], Language.find_by_extension('PKGBUILD.rb') From 8942b54d72e8d10f0ef897159e9098236b27c05d Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Tue, 9 Jul 2024 16:52:28 +0300 Subject: [PATCH 07/21] Update heuristics.yml for Octave language --- lib/linguist/heuristics.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/linguist/heuristics.yml b/lib/linguist/heuristics.yml index 5cfa8197eb..870767fd94 100644 --- a/lib/linguist/heuristics.yml +++ b/lib/linguist/heuristics.yml @@ -418,7 +418,7 @@ disambiguations: - language: Objective-C named_pattern: objectivec - language: Octave - pattern: '^\s*##' + pattern: '^\s*#' - language: Mercury pattern: ':- module' - language: MUF From c8b5dd28d0a87ed1e61c04c9c28a536d6c299f99 Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Tue, 9 Jul 2024 17:05:56 +0300 Subject: [PATCH 08/21] Update heuristics.yml for Octave language --- lib/linguist/heuristics.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/linguist/heuristics.yml b/lib/linguist/heuristics.yml index 870767fd94..5cfa8197eb 100644 --- a/lib/linguist/heuristics.yml +++ b/lib/linguist/heuristics.yml @@ -418,7 +418,7 @@ disambiguations: - language: Objective-C named_pattern: objectivec - language: Octave - pattern: '^\s*#' + pattern: '^\s*##' - language: Mercury pattern: ':- module' - language: MUF From 6355389bd9cc30697ebd66282cff4112513effb4 Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Tue, 9 Jul 2024 19:00:31 +0300 Subject: [PATCH 09/21] Update heuristics.yml for Octave language --- lib/linguist/heuristics.yml | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/lib/linguist/heuristics.yml b/lib/linguist/heuristics.yml index 5cfa8197eb..271b107f65 100644 --- a/lib/linguist/heuristics.yml +++ b/lib/linguist/heuristics.yml @@ -418,7 +418,10 @@ disambiguations: - language: Objective-C named_pattern: objectivec - language: Octave - pattern: '^\s*##' + pattern: + - '^\s*##' + - '^\s*endfunction' + - '^\s*endclassdef' - language: Mercury pattern: ':- module' - language: MUF From 5f3c3fb61c1ceb87a878f0d568737756594144cf Mon Sep 17 00:00:00 2001 From: Andrew Janke Date: Tue, 9 Jul 2024 10:42:21 -0400 Subject: [PATCH 10/21] Add language: Octave, as distinct from Matlab Octave is a "largely Matlab-compatible" language, but it has several syntax elements distinct from Matlab, and which are in common use by Octave coders and encouraged as conventional Octave style. This grammar supports Octave with those variants, and without new Matlab syntax elements which are not present in Octave. --- .gitmodules | 3 +++ vendor/grammars/Octave-language-grammar | 1 + 2 files changed, 4 insertions(+) create mode 160000 vendor/grammars/Octave-language-grammar diff --git a/.gitmodules b/.gitmodules index 05a2491880..f07b3dfb30 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1401,3 +1401,6 @@ [submodule "vendor/grammars/zephir-sublime"] path = vendor/grammars/zephir-sublime url = https://github.com/phalcon/zephir-sublime +[submodule "vendor/grammars/Octave-language-grammar"] + path = vendor/grammars/Octave-language-grammar + url = https://github.com/gnu-octave/Octave-language-grammar diff --git a/vendor/grammars/Octave-language-grammar b/vendor/grammars/Octave-language-grammar new file mode 160000 index 0000000000..0f9f17f60e --- /dev/null +++ b/vendor/grammars/Octave-language-grammar @@ -0,0 +1 @@ +Subproject commit 0f9f17f60e4668c7aaaee6d9a6ad59ce6e773908 From 45e3c6234da8102f6df67f9240043730666658f9 Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Thu, 11 Jul 2024 20:03:44 +0300 Subject: [PATCH 11/21] Update patterns for Octave language in heuristics.yml --- lib/linguist/heuristics.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/lib/linguist/heuristics.yml b/lib/linguist/heuristics.yml index 271b107f65..ea2a069513 100644 --- a/lib/linguist/heuristics.yml +++ b/lib/linguist/heuristics.yml @@ -422,6 +422,9 @@ disambiguations: - '^\s*##' - '^\s*endfunction' - '^\s*endclassdef' + - '^\s*endif' + - '^\s*endfor' + - '^\s*endwhile' - language: Mercury pattern: ':- module' - language: MUF From 25b9b87c7f915d16c864877ca8062ab033384568 Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Fri, 12 Jul 2024 12:39:55 +0300 Subject: [PATCH 12/21] Delete samples/Octave/BetaDistribution.m --- samples/Octave/BetaDistribution.m | 824 ------------------------------ 1 file changed, 824 deletions(-) delete mode 100644 samples/Octave/BetaDistribution.m diff --git a/samples/Octave/BetaDistribution.m b/samples/Octave/BetaDistribution.m deleted file mode 100644 index c71509099d..0000000000 --- a/samples/Octave/BetaDistribution.m +++ /dev/null @@ -1,824 +0,0 @@ -## Copyright (C) 2024 Andreas Bertsatos -## -## This file is part of the statistics package for GNU Octave. -## -## 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 . - -classdef BetaDistribution - ## -*- texinfo -*- - ## @deftypefn {statistics} BetaDistribution - ## - ## Beta probability distribution object. - ## - ## A @code{BetaDistribution} object consists of parameters, a model - ## description, and sample data for a beta probability distribution. - ## - ## The beta distribution uses the following parameters. - ## - ## @multitable @columnfractions 0.25 0.48 0.27 - ## @headitem @var{Parameter} @tab @var{Description} @tab @var{Support} - ## - ## @item @qcode{a} @tab 1st Shape parameter @tab @math{α > 0} - ## @item @qcode{b} @tab 2nd Shape parameter @tab @math{β > 0} - ## @end multitable - ## - ## There are several ways to create a @code{BetaDistribution} object. - ## - ## @itemize - ## @item Fit a distribution to data using the @code{fitdist} function. - ## @item Create a distribution with specified parameter values using the - ## @code{makedist} function. - ## @item Use the constructor @qcode{BetaDistribution (@var{a}, @var{b})} - ## to create a beta distribution with specified parameter values. - ## @item Use the static method @qcode{BetaDistribution.fit (@var{x}, - ## @var{censor}, @var{freq}, @var{options})} to a distribution to data @var{x}. - ## @end itemize - ## - ## It is highly recommended to use @code{fitdist} and @code{makedist} - ## functions to create probability distribution objects, instead of the - ## constructor and the aforementioned static method. - ## - ## A @code{BetaDistribution} object contains the following properties, - ## which can be accessed using dot notation. - ## - ## @multitable @columnfractions 0.25 0.25 0.25 0.25 - ## @item @qcode{DistributionName} @tab @qcode{DistributionCode} @tab - ## @qcode{NumParameters} @tab @qcode{ParameterNames} - ## @item @qcode{ParameterDescription} @tab @qcode{ParameterValues} @tab - ## @qcode{ParameterValues} @tab @qcode{ParameterCI} - ## @item @qcode{ParameterIsFixed} @tab @qcode{Truncation} @tab - ## @qcode{IsTruncated} @tab @qcode{InputData} - ## @end multitable - ## - ## A @code{BetaDistribution} object contains the following methods: - ## @code{cdf}, @code{icdf}, @code{iqr}, @code{mean}, @code{median}, - ## @code{negloglik}, @code{paramci}, @code{pdf}, @code{plot}, @code{proflik}, - ## @code{random}, @code{std}, @code{truncate}, @code{var}. - ## - ## Further information about the Beta distribution can be found at - ## @url{https://en.wikipedia.org/wiki/Beta_distribution} - ## - ## @seealso{fitdist, makedist, betacdf, betainv, betapdf, betarnd, lognfit, - ## betalike, betastat} - ## @end deftypefn - - properties (Dependent = true) - a - b - endproperties - - properties (GetAccess = public, Constant = true) - CensoringAllowed = false; - DistributionName = "BetaDistribution"; - DistributionCode = "beta"; - NumParameters = 2; - ParameterNames = {"a", "b"}; - ParameterDescription = {"1st Shape", "2nd Shape"}; - endproperties - - properties (GetAccess = public, Constant = true) - ParameterRange = [realmin, realmin; Inf, Inf]; - ParameterLogCI = [true, true]; - endproperties - - properties (GetAccess = public , SetAccess = protected) - ParameterValues - ParameterCI - ParameterCovariance - ParameterIsFixed - Truncation - IsTruncated - InputData - endproperties - - methods (Hidden) - - function this = BetaDistribution (a, b) - if (nargin == 0) - a = 1; - b = 1; - endif - checkparams (a, b); - this.InputData = []; - this.IsTruncated = false; - this.ParameterValues = [a, b]; - this.ParameterIsFixed = [true, true]; - this.ParameterCovariance = zeros (this.NumParameters); - endfunction - - function display (this) - fprintf ("%s =\n", inputname(1)); - __disp__ (this, "beta distribution"); - endfunction - - function disp (this) - __disp__ (this, "beta distribution"); - endfunction - - function this = set.a (this, a) - checkparams (a, this.b); - this.InputData = []; - this.ParameterValues(1) = a; - this.ParameterIsFixed = [true, true]; - this.ParameterCovariance = zeros (this.NumParameters); - endfunction - - function a = get.a (this) - a = this.ParameterValues(1); - endfunction - - function this = set.b (this, b) - checkparams (this.a, b); - this.InputData = []; - this.ParameterValues(2) = b; - this.ParameterIsFixed = [true, true]; - this.ParameterCovariance = zeros (this.NumParameters); - endfunction - - function b = get.b (this) - b = this.ParameterValues(2); - endfunction - - endmethods - - methods (Access = public) - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{p} =} cdf (@var{pd}, @var{x}) - ## @deftypefnx {BetaDistribution} {@var{p} =} cdf (@var{pd}, @var{x}, @qcode{"upper"}) - ## - ## Compute the cumulative distribution function (CDF). - ## - ## @code{@var{p} = cdf (@var{pd}, @var{x})} computes the CDF of the - ## probability distribution object, @var{pd}, evaluated at the values in - ## @var{x}. - ## - ## @code{@var{p} = cdf (@dots{}, @qcode{"upper"})} returns the complement of - ## the CDF of the probability distribution object, @var{pd}, evaluated at - ## the values in @var{x}. - ## - ## @end deftypefn - function p = cdf (this, x, uflag) - if (! isscalar (this)) - error ("cdf: requires a scalar probability distribution."); - endif - ## Check for "upper" flag - if (nargin > 2 && strcmpi (uflag, "upper")) - utail = true; - elseif (nargin > 2 && ! strcmpi (uflag, "upper")) - error ("cdf: invalid argument for upper tail."); - else - utail = false; - endif - ## Do the computations - p = betacdf (x, this.a, this.b); - if (this.IsTruncated) - lx = this.Truncation(1); - lb = x < lx; - ux = this.Truncation(2); - ub = x > ux; - p(lb) = 0; - p(ub) = 1; - p(! (lb | ub)) -= betacdf (lx, this.a, this.b); - p(! (lb | ub)) /= diff (betacdf ([lx, ux], this.a, this.b)); - endif - ## Apply uflag - if (utail) - p = 1 - p; - endif - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{p} =} icdf (@var{pd}, @var{p}) - ## - ## Compute the cumulative distribution function (CDF). - ## - ## @code{@var{p} = icdf (@var{pd}, @var{x})} computes the quantile (the - ## inverse of the CDF) of the probability distribution object, @var{pd}, - ## evaluated at the values in @var{x}. - ## - ## @end deftypefn - function x = icdf (this, p) - if (! isscalar (this)) - error ("icdf: requires a scalar probability distribution."); - endif - if (this.IsTruncated) - lp = betacdf (this.Truncation(1), this.a, this.b); - up = betacdf (this.Truncation(2), this.a, this.b); - ## Adjust p values within range of p @ lower limit and p @ upper limit - is_nan = p < 0 | p > 1; - p(is_nan) = NaN; - np = lp + (up - lp) .* p; - x = betainv (np, this.a, this.b); - x(x < this.Truncation(1)) = this.Truncation(1); - x(x > this.Truncation(2)) = this.Truncation(2); - else - x = betainv (p, this.a, this.b); - endif - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{r} =} iqr (@var{pd}) - ## - ## Compute the interquartile range of a probability distribution. - ## - ## @code{@var{r} = iqr (@var{pd})} computes the interquartile range of the - ## probability distribution object, @var{pd}. - ## - ## @end deftypefn - function r = iqr (this) - if (! isscalar (this)) - error ("iqr: requires a scalar probability distribution."); - endif - r = diff (icdf (this, [0.25, 0.75])); - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{m} =} mean (@var{pd}) - ## - ## Compute the mean of a probability distribution. - ## - ## @code{@var{m} = mean (@var{pd})} computes the mean of the probability - ## distribution object, @var{pd}. - ## - ## @end deftypefn - function m = mean (this) - if (! isscalar (this)) - error ("mean: requires a scalar probability distribution."); - endif - if (this.IsTruncated) - fm = @(x) x .* pdf (this, x); - m = integral (fm, this.Truncation(1), this.Truncation(2)); - else - m = betastat (this.a, this.b); - endif - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{m} =} median (@var{pd}) - ## - ## Compute the median of a probability distribution. - ## - ## @code{@var{m} = median (@var{pd})} computes the median of the probability - ## distribution object, @var{pd}. - ## - ## @end deftypefn - function m = median (this) - if (! isscalar (this)) - error ("median: requires a scalar probability distribution."); - endif - if (this.IsTruncated) - lx = this.Truncation(1); - ux = this.Truncation(2); - Fa_b = betacdf ([lx, ux], this.a, this.b); - m = betainv (sum (Fa_b) / 2, this.a, this.b); - else - m = betainv (0.5, this.a, this.b); - endif - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{nlogL} =} negloglik (@var{pd}) - ## - ## Compute the negative loglikelihood of a probability distribution. - ## - ## @code{@var{m} = negloglik (@var{pd})} computes the negative loglikelihood - ## of the probability distribution object, @var{pd}. - ## - ## @end deftypefn - function nlogL = negloglik (this) - if (! isscalar (this)) - error ("negloglik: requires a scalar probability distribution."); - endif - if (isempty (this.InputData)) - nlogL = []; - return - endif - nlogL = - betalike ([this.a, this.b], this.InputData.data, ... - this.InputData.freq); - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{ci} =} paramci (@var{pd}) - ## @deftypefnx {BetaDistribution} {@var{ci} =} paramci (@var{pd}, @var{Name}, @var{Value}) - ## - ## Compute the confidence intervals for probability distribution parameters. - ## - ## @code{@var{ci} = paramci (@var{pd})} computes the lower and upper - ## boundaries of the 95% confidence interval for each parameter of the - ## probability distribution object, @var{pd}. - ## - ## @code{@var{ci} = paramci (@var{pd}, @var{Name}, @var{Value})} computes the - ## confidence intervals with additional options specified specified by - ## @qcode{Name-Value} pair arguments listed below. - ## - ## @multitable @columnfractions 0.18 0.02 0.8 - ## @headitem @var{Name} @tab @tab @var{Value} - ## - ## @item @qcode{"Alpha"} @tab @tab A scalar value in the range @math{(0,1)} - ## specifying the significance level for the confidence interval. The - ## default value 0.05 corresponds to a 95% confidence interval. - ## - ## @item @qcode{"Parameter"} @tab @tab A character vector or a cell array of - ## character vectors specifying the parameter names for which to compute - ## confidence intervals. By default, @code{paramci} computes confidence - ## intervals for all distribution parameters. - ## @end multitable - ## - ## @code{paramci} is meaningful only when @var{pd} is fitted to data, - ## otherwise an empty array, @qcode{[]}, is returned. - ## - ## @end deftypefn - function ci = paramci (this, varargin) - if (! isscalar (this)) - error ("paramci: requires a scalar probability distribution."); - endif - if (isempty (this.InputData)) - ci = [this.ParameterValues; this.ParameterValues]; - else - ci = __paramci__ (this, varargin{:}); - endif - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{y} =} pdf (@var{pd}, @var{x}) - ## - ## Compute the probability distribution function (PDF). - ## - ## @code{@var{y} = pdf (@var{pd}, @var{x})} computes the PDF of the - ## probability distribution object, @var{pd}, evaluated at the values in - ## @var{x}. - ## - ## @end deftypefn - function y = pdf (this, x) - if (! isscalar (this)) - error ("pdf: requires a scalar probability distribution."); - endif - y = betapdf (x, this.a, this.b); - if (this.IsTruncated) - lx = this.Truncation(1); - lb = x < lx; - ux = this.Truncation(2); - ub = x > ux; - y(lb | ub) = 0; - y(! (lb | ub)) /= diff (betacdf ([lx, ux], this.a, this.b)); - endif - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {} plot (@var{pd}) - ## @deftypefnx {BetaDistribution} {} plot (@var{pd}, @var{Name}, @var{Value}) - ## @deftypefnx {BetaDistribution} {@var{h} =} plot (@dots{}) - ## - ## Plot a probability distribution object. - ## - ## @code{plot (@var{pd}} plots a probability density function (PDF) of the - ## probability distribution object @var{pd}. If @var{pd} contains data, - ## which have been fitted by @code{fitdist}, the PDF is superimposed over a - ## histogram of the data. - ## - ## @code{plot (@var{pd}, @var{Name}, @var{Value})} specifies additional - ## options with the @qcode{Name-Value} pair arguments listed below. - ## - ## @multitable @columnfractions 0.18 0.02 0.8 - ## @headitem @tab @var{Name} @tab @var{Value} - ## - ## @item @qcode{"PlotType"} @tab @tab A character vector specifying the plot - ## type. @qcode{"pdf"} plots the probability density function (PDF). When - ## @var{pd} is fit to data, the PDF is superimposed on a histogram of the - ## data. @qcode{"cdf"} plots the cumulative density function (CDF). When - ## @var{pd} is fit to data, the CDF is superimposed over an empirical CDF. - ## @qcode{"probability"} plots a probability plot using a CDF of the data - ## and a CDF of the fitted probability distribution. This option is - ## available only when @var{pd} is fitted to data. - ## - ## @item @qcode{"Discrete"} @tab @tab A logical scalar to specify whether to - ## plot the PDF or CDF of a discrete distribution object as a line plot or a - ## stem plot, by specifying @qcode{false} or @qcode{true}, respectively. By - ## default, it is @qcode{true} for discrete distributions and @qcode{false} - ## for continuous distributions. When @var{pd} is a continuous distribution - ## object, option is ignored. - ## - ## @item @qcode{"Parent"} @tab @tab An axes graphics object for plot. If - ## not specified, the @code{plot} function plots into the current axes or - ## creates a new axes object if one does not exist. - ## @end multitable - ## - ## @code{@var{h} = plot (@dots{})} returns a graphics handle to the plotted - ## objects. - ## - ## @end deftypefn - function [varargout] = plot (this, varargin) - if (! isscalar (this)) - error ("plot: requires a scalar probability distribution."); - endif - h = __plot__ (this, false, varargin{:}); - if (nargout > 0) - varargout{1} = h; - endif - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {[@var{nlogL}, @var{param}] =} proflik (@var{pd}, @var{pnum}) - ## @deftypefnx {BetaDistribution} {[@var{nlogL}, @var{param}] =} proflik (@var{pd}, @var{pnum}, @qcode{"Display"}, @var{display}) - ## @deftypefnx {BetaDistribution} {[@var{nlogL}, @var{param}] =} proflik (@var{pd}, @var{pnum}, @var{setparam}) - ## @deftypefnx {BetaDistribution} {[@var{nlogL}, @var{param}] =} proflik (@var{pd}, @var{pnum}, @var{setparam}, @qcode{"Display"}, @var{display}) - ## - ## Profile likelihood function for a probability distribution object. - ## - ## @code{[@var{nlogL}, @var{param}] = proflik (@var{pd}, @var{pnum})} - ## returns a vector @var{nlogL} of negative loglikelihood values and a - ## vector @var{param} of corresponding parameter values for the parameter in - ## the position indicated by @var{pnum}. By default, @code{proflik} uses - ## the lower and upper bounds of the 95% confidence interval and computes - ## 100 equispaced values for the selected parameter. @var{pd} must be - ## fitted to data. - ## - ## @code{[@var{nlogL}, @var{param}] = proflik (@var{pd}, @var{pnum}, - ## @qcode{"Display"}, @qcode{"on"})} also plots the profile likelihood - ## against the default range of the selected parameter. - ## - ## @code{[@var{nlogL}, @var{param}] = proflik (@var{pd}, @var{pnum}, - ## @var{setparam})} defines a user-defined range of the selected parameter. - ## - ## @code{[@var{nlogL}, @var{param}] = proflik (@var{pd}, @var{pnum}, - ## @var{setparam}, @qcode{"Display"}, @qcode{"on"})} also plots the profile - ## likelihood against the user-defined range of the selected parameter. - ## - ## For the beta distribution, @qcode{@var{pnum} = 1} selects the parameter - ## @qcode{a} and @qcode{@var{pnum} = 2} selects the parameter @var{b}. - ## - ## When opted to display the profile likelihood plot, @code{proflik} also - ## plots the baseline loglikelihood computed at the lower bound of the 95% - ## confidence interval and estimated maximum likelihood. The latter might - ## not be observable if it is outside of the used-defined range of parameter - ## values. - ## - ## @end deftypefn - function [varargout] = proflik (this, pnum, varargin) - if (! isscalar (this)) - error ("proflik: requires a scalar probability distribution."); - endif - if (isempty (this.InputData)) - error ("proflik: no fitted data available."); - endif - [varargout{1:nargout}] = __proflik__ (this, pnum, varargin{:}); - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{y} =} random (@var{pd}) - ## @deftypefnx {BetaDistribution} {@var{y} =} random (@var{pd}, @var{rows}) - ## @deftypefnx {BetaDistribution} {@var{y} =} random (@var{pd}, @var{rows}, @var{cols}, @dots{}) - ## @deftypefnx {BetaDistribution} {@var{y} =} random (@var{pd}, [@var{sz}]) - ## - ## Generate random arrays from the probability distribution object. - ## - ## @code{@var{r} = random (@var{pd})} returns a random number from the - ## distribution object @var{pd}. - ## - ## When called with a single size argument, @code{betarnd} returns a square - ## matrix with the dimension specified. When called with more than one - ## scalar argument, the first two arguments are taken as the number of rows - ## and columns and any further arguments specify additional matrix - ## dimensions. The size may also be specified with a row vector of - ## dimensions, @var{sz}. - ## - ## @end deftypefn - function r = random (this, varargin) - if (! isscalar (this)) - error ("random: requires a scalar probability distribution."); - endif - if (this.IsTruncated) - sz = [varargin{:}]; - ps = prod (sz); - ## Get an estimate of how many more random numbers we need to randomly - ## pick the appropriate size from - lx = this.Truncation(1); - ux = this.Truncation(2); - ratio = 1 / diff (betacdf ([lx, ux], this.a, this.b)); - nsize = fix (2 * ratio * ps); # times 2 to be on the safe side - ## Generate the numbers and remove out-of-bound random samples - r = betarnd (this.a, this.b, nsize, 1); - r(r < lx | r > ux) = []; - ## Randomly select the required size and reshape to requested dimensions - idx = randperm (numel (r), ps); - r = reshape (r(idx), sz); - else - r = betarnd (this.a, this.b, varargin{:}); - endif - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{s} =} std (@var{pd}) - ## - ## Compute the standard deviation of a probability distribution. - ## - ## @code{@var{s} = std (@var{pd})} computes the standard deviation of the - ## probability distribution object, @var{pd}. - ## - ## @end deftypefn - function s = std (this) - if (! isscalar (this)) - error ("std: requires a scalar probability distribution."); - endif - v = var (this); - s = sqrt (v); - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{t} =} truncate (@var{pd}, @var{lower}, @var{upper}) - ## - ## Truncate a probability distribution. - ## - ## @code{@var{t} = truncate (@var{pd})} returns a probability distribution - ## @var{t}, which is the probability distribution @var{pd} truncated to the - ## specified interval with lower limit, @var{lower}, and upper limit, - ## @var{upper}. If @var{pd} is fitted to data with @code{fitdist}, the - ## returned probability distribution @var{t} is not fitted, does not contain - ## any data or estimated values, and it is as it has been created with the - ## @var{makedist} function, but it includes the truncation interval. - ## - ## @end deftypefn - function this = truncate (this, lower, upper) - if (! isscalar (this)) - error ("truncate: requires a scalar probability distribution."); - endif - if (nargin < 3) - error ("truncate: missing input argument."); - elseif (lower >= upper) - error ("truncate: invalid lower upper limits."); - endif - this.Truncation = [lower, upper]; - this.IsTruncated = true; - this.InputData = []; - this.ParameterIsFixed = [true, true]; - this.ParameterCovariance = zeros (this.NumParameters); - endfunction - - ## -*- texinfo -*- - ## @deftypefn {BetaDistribution} {@var{v} =} var (@var{pd}) - ## - ## Compute the variance of a probability distribution. - ## - ## @code{@var{v} = var (@var{pd})} computes the standard deviation of the - ## probability distribution object, @var{pd}. - ## - ## @end deftypefn - function v = var (this) - if (! isscalar (this)) - error ("var: requires a scalar probability distribution."); - endif - if (this.IsTruncated) - fm = @(x) x .* pdf (this, x); - m = integral (fm, this.Truncation(1), this.Truncation(2)); - fv = @(x) ((x - m) .^ 2) .* pdf (this, x); - v = integral (fv, this.Truncation(1), this.Truncation(2)); - else - [~, v] = betastat (this.a, this.b); - endif - endfunction - - endmethods - - methods (Static, Hidden) - - function pd = fit (x, varargin) - ## Check input arguments - if (nargin < 2) - alpha = 0.05; - else - alpha = varargin{1}; - endif - if (nargin < 3) - freq = []; - else - freq = varargin{2}; - endif - if (nargin < 4) - options.Display = "off"; - options.MaxFunEvals = 400; - options.MaxIter = 200; - options.TolX = 1e-6; - else - options = varargin{3}; - endif - ## Fit data - [phat, pci] = betafit (x, alpha, freq, options); - [~, acov] = betalike (phat, x, freq); - ## Create fitted distribution object - pd = BetaDistribution.makeFitted (phat, pci, acov, x, freq); - endfunction - - function pd = makeFitted (phat, pci, acov, x, freq) - a = phat(1); - b = phat(2); - pd = BetaDistribution (a, b); - pd.ParameterCI = pci; - pd.ParameterIsFixed = [false, false]; - pd.ParameterCovariance = acov; - pd.InputData = struct ("data", x, "cens", [], "freq", freq); - endfunction - - endmethods - -endclassdef - -function checkparams (a, b) - if (! (isscalar (a) && isnumeric (a) && isreal (a) && isfinite (a) - && a > 0)) - error ("BetaDistribution: A must be a positive real scalar.") - endif - if (! (isscalar (b) && isnumeric (b) && isreal (b) - && isfinite (b) && b > 0)) - error ("BetaDistribution: B must be a positive real scalar.") - endif -endfunction - -## Test output -%!shared pd, t -%! pd = BetaDistribution; -%! t = truncate (pd, 0.2, 0.8); -%!assert (cdf (pd, [0:0.2:1]), [0, 0.2, 0.4, 0.6, 0.8, 1], 1e-4); -%!assert (cdf (t, [0:0.2:1]), [0, 0, 0.3333, 0.6667, 1, 1], 1e-4); -%!assert (cdf (pd, [-1, 1, NaN]), [0, 1, NaN], 1e-4); -%!assert (cdf (t, [-1, 1, NaN]), [0, 1, NaN], 1e-4); -%!assert (icdf (pd, [0:0.2:1]), [0, 0.2, 0.4, 0.6, 0.8, 1], 1e-4); -%!assert (icdf (t, [0:0.2:1]), [0.2, 0.32, 0.44, 0.56, 0.68, 0.8], 1e-4); -%!assert (icdf (pd, [-1, 0.4:0.2:1, NaN]), [NaN, 0.4, 0.6, 0.8, 1, NaN], 1e-4); -%!assert (icdf (t, [-1, 0.4:0.2:1, NaN]), [NaN, 0.44, 0.56, 0.68, 0.8, NaN], 1e-4); -%!assert (iqr (pd), 0.5, 1e-4); -%!assert (iqr (t), 0.3, 1e-4); -%!assert (mean (pd), 0.5); -%!assert (mean (t), 0.5, 1e-6); -%!assert (median (pd), 0.5); -%!assert (median (t), 0.5, 1e-6); -%!assert (pdf (pd, [0:0.2:1]), [1, 1, 1, 1, 1, 1], 1e-4); -%!assert (pdf (t, [0:0.2:1]), [0, 1.6667, 1.6667, 1.6667, 1.6667, 0], 1e-4); -%!assert (pdf (pd, [-1, 1, NaN]), [0, 1, NaN], 1e-4); -%!assert (pdf (t, [-1, 1, NaN]), [0, 0, NaN], 1e-4); -%!assert (isequal (size (random (pd, 100, 50)), [100, 50])) -%!assert (any (random (t, 1000, 1) < 0.2), false); -%!assert (any (random (t, 1000, 1) > 0.8), false); -%!assert (std (pd), 0.2887, 1e-4); -%!assert (std (t), 0.1732, 1e-4); -%!assert (var (pd), 0.0833, 1e-4); -%!assert (var (t), 0.0300, 1e-4); - -## Test input validation -## 'BetaDistribution' constructor -%!error ... -%! BetaDistribution(0, 1) -%!error ... -%! BetaDistribution(Inf, 1) -%!error ... -%! BetaDistribution(i, 1) -%!error ... -%! BetaDistribution("a", 1) -%!error ... -%! BetaDistribution([1, 2], 1) -%!error ... -%! BetaDistribution(NaN, 1) -%!error ... -%! BetaDistribution(1, 0) -%!error ... -%! BetaDistribution(1, -1) -%!error ... -%! BetaDistribution(1, Inf) -%!error ... -%! BetaDistribution(1, i) -%!error ... -%! BetaDistribution(1, "a") -%!error ... -%! BetaDistribution(1, [1, 2]) -%!error ... -%! BetaDistribution(1, NaN) - -## 'cdf' method -%!error ... -%! cdf (BetaDistribution, 2, "uper") -%!error ... -%! cdf (BetaDistribution, 2, 3) - -## 'paramci' method -%!shared x -%! randg ("seed", 1); -%! x = betarnd (1, 1, [100, 1]); -%!error ... -%! paramci (BetaDistribution.fit (x), "alpha") -%!error ... -%! paramci (BetaDistribution.fit (x), "alpha", 0) -%!error ... -%! paramci (BetaDistribution.fit (x), "alpha", 1) -%!error ... -%! paramci (BetaDistribution.fit (x), "alpha", [0.5 2]) -%!error ... -%! paramci (BetaDistribution.fit (x), "alpha", "") -%!error ... -%! paramci (BetaDistribution.fit (x), "alpha", {0.05}) -%!error ... -%! paramci (BetaDistribution.fit (x), "parameter", "a", "alpha", {0.05}) -%!error ... -%! paramci (BetaDistribution.fit (x), "parameter", {"a", "b", "param"}) -%!error ... -%! paramci (BetaDistribution.fit (x), "alpha", 0.01, ... -%! "parameter", {"a", "b", "param"}) -%!error ... -%! paramci (BetaDistribution.fit (x), "parameter", "param") -%!error ... -%! paramci (BetaDistribution.fit (x), "alpha", 0.01, "parameter", "param") -%!error ... -%! paramci (BetaDistribution.fit (x), "NAME", "value") -%!error ... -%! paramci (BetaDistribution.fit (x), "alpha", 0.01, "NAME", "value") -%!error ... -%! paramci (BetaDistribution.fit (x), "alpha", 0.01, "parameter", "a", ... -%! "NAME", "value") - -## 'plot' method -%!error ... -%! plot (BetaDistribution, "Parent") -%!error ... -%! plot (BetaDistribution, "PlotType", 12) -%!error ... -%! plot (BetaDistribution, "PlotType", {"pdf", "cdf"}) -%!error ... -%! plot (BetaDistribution, "PlotType", "pdfcdf") -%!error ... -%! plot (BetaDistribution, "Discrete", "pdfcdf") -%!error ... -%! plot (BetaDistribution, "Discrete", [1, 0]) -%!error ... -%! plot (BetaDistribution, "Discrete", {true}) -%!error ... -%! plot (BetaDistribution, "Parent", 12) -%!error ... -%! plot (BetaDistribution, "Parent", "hax") - -## 'proflik' method -%!error ... -%! proflik (BetaDistribution, 2) -%!error ... -%! proflik (BetaDistribution.fit (x), 3) -%!error ... -%! proflik (BetaDistribution.fit (x), [1, 2]) -%!error ... -%! proflik (BetaDistribution.fit (x), {1}) -%!error ... -%! proflik (BetaDistribution.fit (x), 1, ones (2)) -%!error ... -%! proflik (BetaDistribution.fit (x), 1, "Display") -%!error ... -%! proflik (BetaDistribution.fit (x), 1, "Display", 1) -%!error ... -%! proflik (BetaDistribution.fit (x), 1, "Display", {1}) -%!error ... -%! proflik (BetaDistribution.fit (x), 1, "Display", {"on"}) -%!error ... -%! proflik (BetaDistribution.fit (x), 1, "Display", ["on"; "on"]) -%!error ... -%! proflik (BetaDistribution.fit (x), 1, "Display", "onnn") -%!error ... -%! proflik (BetaDistribution.fit (x), 1, "NAME", "on") -%!error ... -%! proflik (BetaDistribution.fit (x), 1, {"NAME"}, "on") -%!error ... -%! proflik (BetaDistribution.fit (x), 1, {[1 2 3 4]}, "Display", "on") - -## 'truncate' method -%!error ... -%! truncate (BetaDistribution) -%!error ... -%! truncate (BetaDistribution, 2) -%!error ... -%! truncate (BetaDistribution, 4, 2) - -## Catch errors when using array of probability objects with available methods -%!shared pd -%! pd = BetaDistribution(1, 1); -%! pd(2) = BetaDistribution(1, 3); -%!error cdf (pd, 1) -%!error icdf (pd, 0.5) -%!error iqr (pd) -%!error mean (pd) -%!error median (pd) -%!error negloglik (pd) -%!error paramci (pd) -%!error pdf (pd, 1) -%!error plot (pd) -%!error proflik (pd, 2) -%!error random (pd) -%!error std (pd) -%!error ... -%! truncate (pd, 2, 4) -%!error var (pd) From 0aa5feda5ddc93299c574fdbebd497c5d78dae79 Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Fri, 12 Jul 2024 12:41:11 +0300 Subject: [PATCH 13/21] Delete samples/Octave/accumarray.m --- samples/Octave/accumarray.m | 466 ------------------------------------ 1 file changed, 466 deletions(-) delete mode 100644 samples/Octave/accumarray.m diff --git a/samples/Octave/accumarray.m b/samples/Octave/accumarray.m deleted file mode 100644 index 4f850c12ea..0000000000 --- a/samples/Octave/accumarray.m +++ /dev/null @@ -1,466 +0,0 @@ -######################################################################## -## -## Copyright (C) 2007-2024 The Octave Project Developers -## -## See the file COPYRIGHT.md in the top-level directory of this -## distribution or . -## -## This file is part of Octave. -## -## 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 3 of the License, or -## (at your option) any later version. -## -## Octave 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 -## . -## -######################################################################## - -## -*- texinfo -*- -## @deftypefn {} {@var{A} =} accumarray (@var{subs}, @var{vals}) -## @deftypefnx {} {@var{A} =} accumarray (@var{subs}, @var{vals}, @var{sz}) -## @deftypefnx {} {@var{A} =} accumarray (@var{subs}, @var{vals}, @var{sz}, @var{fcn}) -## @deftypefnx {} {@var{A} =} accumarray (@var{subs}, @var{vals}, @var{sz}, @var{fcn}, @var{fillval}) -## @deftypefnx {} {@var{A} =} accumarray (@var{subs}, @var{vals}, @var{sz}, @var{fcn}, @var{fillval}, @var{issparse}) -## -## Create an array by accumulating the elements of a vector into the -## positions defined by their subscripts. -## -## The subscripts are defined by the rows of the matrix @var{subs} and the -## values by @var{vals}. Each row of @var{subs} corresponds to one of the -## values in @var{vals}. If @var{vals} is a scalar, it will be used for each -## of the row of @var{subs}. If @var{subs} is a cell array of vectors, all -## vectors must be of the same length, and the subscripts in the @var{k}th -## vector must correspond to the @var{k}th dimension of the result. -## -## The size of the matrix will be determined by the subscripts -## themselves. However, if @var{sz} is defined it determines the matrix -## size. The length of @var{sz} must correspond to the number of columns -## in @var{subs}. An exception is if @var{subs} has only one column, in -## which case @var{sz} may be the dimensions of a vector and the -## subscripts of @var{subs} are taken as the indices into it. -## -## The default action of @code{accumarray} is to sum the elements with -## the same subscripts. This behavior can be modified by defining the -## @var{fcn} function. This should be a function or function handle -## that accepts a column vector and returns a scalar. The result of the -## function should not depend on the order of the subscripts. -## -## The elements of the returned array that have no subscripts associated -## with them are set to zero. Defining @var{fillval} to some other value -## allows these values to be defined. This behavior changes, however, -## for certain values of @var{fcn}. If @var{fcn} is @code{@@min} -## (respectively, @code{@@max}) then the result will be filled with the -## minimum (respectively, maximum) integer if @var{vals} is of integral -## type, logical false (respectively, logical true) if @var{vals} is of -## logical type, zero if @var{fillval} is zero and all values are -## non-positive (respectively, non-negative), and NaN otherwise. -## -## By default @code{accumarray} returns a full matrix. If -## @var{issparse} is logically true, then a sparse matrix is returned -## instead. -## -## The following @code{accumarray} example constructs a frequency table -## that in the first column counts how many occurrences each number in -## the second column has, taken from the vector @var{x}. Note the usage -## of @code{unique} for assigning to all repeated elements of @var{x} -## the same index (@pxref{XREFunique,,@code{unique}}). -## -## @example -## @group -## @var{x} = [91, 92, 90, 92, 90, 89, 91, 89, 90, 100, 100, 100]; -## [@var{u}, ~, @var{j}] = unique (@var{x}); -## [accumarray(@var{j}', 1), @var{u}'] -## @result{} 2 89 -## 3 90 -## 2 91 -## 2 92 -## 3 100 -## @end group -## @end example -## -## Another example, where the result is a multi-dimensional 3-D array and -## the default value (zero) appears in the output: -## -## @example -## @group -## accumarray ([1, 1, 1; -## 2, 1, 2; -## 2, 3, 2; -## 2, 1, 2; -## 2, 3, 2], 101:105) -## @result{} ans(:,:,1) = [101, 0, 0; 0, 0, 0] -## @result{} ans(:,:,2) = [0, 0, 0; 206, 0, 208] -## @end group -## @end example -## -## The sparse option can be used as an alternative to the @code{sparse} -## constructor (@pxref{XREFsparse,,@code{sparse}}). Thus -## -## @example -## sparse (@var{i}, @var{j}, @var{sv}) -## @end example -## -## @noindent -## can be written with @code{accumarray} as -## -## @example -## accumarray ([@var{i}, @var{j}], @var{sv}', [], [], 0, true) -## @end example -## -## @noindent -## For repeated indices, @code{sparse} adds the corresponding value. To -## take the minimum instead, use @code{min} as an accumulator function: -## -## @example -## accumarray ([@var{i}, @var{j}], @var{sv}', [], @@min, 0, true) -## @end example -## -## The complexity of accumarray in general for the non-sparse case is -## generally O(M+N), where N is the number of subscripts and M is the -## maximum subscript (linearized in multi-dimensional case). If -## @var{fcn} is one of @code{@@sum} (default), @code{@@max}, -## @code{@@min} or @code{@@(x) @{x@}}, an optimized code path is used. -## Note that for general reduction function the interpreter overhead can -## play a major part and it may be more efficient to do multiple -## accumarray calls and compute the results in a vectorized manner. -## -## @seealso{accumdim, unique, sparse} -## @end deftypefn - -function A = accumarray (subs, vals, sz = [], fcn = [], fillval = [], issparse = []) - - if (nargin < 2) - print_usage (); - endif - - lenvals = length (vals); - - if (iscell (subs)) - subs = cellfun (@vec, subs, "uniformoutput", false); - ndims = numel (subs); - if (ndims == 1) - subs = subs{1}; - endif - - lensubs = cellfun (@length, subs); - - if (any (lensubs != lensubs(1)) || (lenvals > 1 && lenvals != lensubs(1))) - error ("accumarray: dimension mismatch"); - endif - - else - ndims = columns (subs); - if (lenvals > 1 && lenvals != rows (subs)) - error ("accumarray: dimension mismatch"); - endif - endif - - if (isempty (fcn)) - fcn = @sum; - elseif (! is_function_handle (fcn)) - error ("accumarray: FCN must be a function handle"); - endif - - if (isempty (fillval)) - fillval = 0; - endif - - if (isempty (issparse)) - issparse = false; - endif - - if (issparse) - - ## Sparse case. - ## Avoid linearizing the subscripts, because it could overflow. - - if (fillval != 0) - error ("accumarray: FILLVAL must be zero in the sparse case"); - endif - - ## Ensure subscripts are a two-column matrix. - if (iscell (subs)) - subs = [subs{:}]; - endif - - ## Validate dimensions. - if (ndims == 1) - subs(:,2) = 1; - elseif (ndims != 2) - error ("accumarray: in the sparse case, needs 1 or 2 subscripts"); - endif - - if (isnumeric (vals) || islogical (vals)) - vals = double (vals); - else - error ("accumarray: in the sparse case, values must be numeric or logical"); - endif - - if (fcn != @sum) - - ## Reduce values. This is not needed if we're about to sum them, - ## because "sparse" can do that. - - ## Sort indices. - [subs, idx] = sortrows (subs); - n = rows (subs); - ## Identify runs. - jdx = find (any (diff (subs, 1, 1), 2)); - jdx = [jdx; n]; - - vals = cellfun (fcn, mat2cell (vals(:)(idx), diff ([0; jdx]))); - subs = subs(jdx, :); - mode = "unique"; - else - mode = "sum"; - endif - - ## Form the sparse matrix. - if (isempty (sz)) - A = sparse (subs(:,1), subs(:,2), vals, mode); - elseif (length (sz) == 2) - - ## Row vector case - if (sz(1) == 1) - [i, j] = deal (subs(:,2), subs(:,1)); - else - [i, j] = deal (subs(:,1), subs(:,2)); - endif - A = sparse (i, j, vals, sz(1), sz(2), mode); - else - error ("accumarray: dimensions mismatch"); - endif - - else - - ## Linearize subscripts. - if (ndims > 1) - if (isempty (sz)) - if (iscell (subs)) - sz = cellfun ("max", subs); - else - sz = max (subs, [], 1); - endif - elseif (ndims != length (sz)) - error ("accumarray: dimensions mismatch"); - endif - - ## Convert multidimensional subscripts. - if (isnumeric (subs)) - subs = num2cell (subs, 1); - endif - subs = sub2ind (sz, subs{:}); # creates index cache - elseif (! isempty (sz) && length (sz) < 2) - error ("accumarray: needs at least 2 dimensions"); - elseif (! isindex (subs)) # creates index cache - error ("accumarray: indices must be positive integers"); - endif - - - ## Some built-in reductions handled efficiently. - - if (fcn == @sum) - ## Fast summation. - if (isempty (sz)) - A = __accumarray_sum__ (subs, vals); - else - A = __accumarray_sum__ (subs, vals, prod (sz)); - ## set proper shape. - A = reshape (A, sz); - endif - - ## we fill in nonzero fill value. - if (fillval != 0) - mask = true (size (A)); - mask(subs) = false; - A(mask) = fillval; - endif - elseif (fcn == @max) - ## Fast maximization. - - if (isinteger (vals)) - zero = intmin (vals); - elseif (islogical (vals)) - zero = false; - elseif (fillval == 0 && all (vals(:) >= 0)) - ## This is a common case - fillval is zero, all numbers - ## nonegative. - zero = 0; - else - zero = NaN; # Neutral value. - endif - - if (isempty (sz)) - A = __accumarray_max__ (subs, vals, zero); - else - A = __accumarray_max__ (subs, vals, zero, prod (sz)); - A = reshape (A, sz); - endif - - if (fillval != zero && ! (isnan (fillval) || isnan (zero))) - mask = true (size (A)); - mask(subs) = false; - A(mask) = fillval; - endif - elseif (fcn == @min) - ## Fast minimization. - - if (isinteger (vals)) - zero = intmax (vals); - elseif (islogical (vals)) - zero = true; - elseif (fillval == 0 && all (vals(:) <= 0)) - ## This is a common case - fillval is zero, all numbers - ## non-positive. - zero = 0; - else - zero = NaN; # Neutral value. - endif - - if (isempty (sz)) - A = __accumarray_min__ (subs, vals, zero); - else - A = __accumarray_min__ (subs, vals, zero, prod (sz)); - A = reshape (A, sz); - endif - - if (fillval != zero && ! (isnan (fillval) || isnan (zero))) - mask = true (size (A)); - mask(subs) = false; - A(mask) = fillval; - endif - else - - ## The general case. Reduce values. - n = rows (subs); - if (numel (vals) == 1) - vals = vals(ones (1, n), 1); - else - vals = vals(:); - endif - - ## Sort indices. - [subs, idx] = sort (subs); - ## Identify runs. - jdx = find (subs(1:n-1) != subs(2:n)); - if (n != 0) # bug #47287 - jdx = [jdx; n]; - endif - vals = mat2cell (vals(idx), diff ([0; jdx])); - ## Optimize the case when function is @(x) {x}, i.e., we just want - ## to collect the values to cells. - persistent simple_cell_str = func2str (@(x) {x}); - if (! strcmp (func2str (fcn), simple_cell_str)) - vals = cellfun (fcn, vals); - endif - - subs = subs(jdx); - - if (isempty (sz)) - sz = max (subs); - ## If subs is empty, sz will be too, and length will be 0, hence "<= 1" - if (length (sz) <= 1) - sz(2) = 1; - endif - endif - - ## Construct matrix of fillvals. - if (iscell (vals)) - A = cell (sz); - elseif (fillval == 0) - A = zeros (sz, class (vals)); - else - A = repmat (fillval, sz); - endif - - ## Set the reduced values. - A(subs) = vals; - endif - endif - -endfunction - - -%!assert (accumarray ([1; 2; 4; 2; 4], 101:105), [101; 206; 0; 208]) -%!assert (accumarray ([1 1 1; 2 1 2; 2 3 2; 2 1 2; 2 3 2], 101:105), -%! cat (3, [101 0 0; 0 0 0], [0 0 0; 206 0 208])) - -%!assert (accumarray ([1 1 1; 2 1 2; 2 3 2; 2 1 2; 2 3 2], 101:105, [], @(x) sin (sum (x))), -%! sin (cat (3, [101,0,0;0,0,0],[0,0,0;206,0,208]))) - -%!assert (accumarray ({[1 3 3 2 3 1 2 2 3 3 1 2], [3 4 2 1 4 3 4 2 2 4 3 4], [1 1 2 2 1 1 2 1 1 1 2 2]}, 101:112), -%! cat (3, [0 0 207 0; 0 108 0 0; 0 109 0 317], [0 0 111 0; 104 0 0 219; 0 103 0 0])) - -%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 101:105, [2 4], @max, NaN), -%! [101 NaN NaN NaN; 104 NaN 105 NaN]) - -%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 101:105, [], @prod), -%! [101 0 0; 10608 0 10815]) -%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 101:105, [2 4], @prod, 0, true), -%! sparse ([1 2 2], [1 1 3], [101 10608 10815], 2, 4)) -%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 1, [2 4]), [1 0 0 0; 2 0 2 0]) -%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 101:105, [2 4], @(x) length (x) > 1), -%! [false false false false; true false true false]) - -%!assert (accumarray ([1; 2], [3; 4], [2, 1], @min, [], 0), [3; 4]) -%!assert (accumarray ([1; 2], [3; 4], [2, 1], @min, [], 1), sparse ([3; 4])) -%!assert (accumarray ([1; 2], [3; 4], [1, 2], @min, [], 0), [3, 4]) -%!assert (accumarray ([1; 2], [3; 4], [1, 2], @min, [], 1), sparse ([3, 4])) - -%!test -%! A = accumarray ([1 1; 2 1; 2 3; 2 1; 2 3], 101:105, [2,4], @(x) {x}); -%! assert (A{2},[102; 104]); - -%!test -%! subs = ceil (rand (2000, 3)*10); -%! vals = rand (2000, 1); -%! assert (accumarray (subs, vals, [], @max), -%! accumarray (subs, vals, [], @(x) max (x))); - -%!test -%! subs = ceil (rand (2000, 1)*100); -%! vals = rand (2000, 1); -%! assert (accumarray (subs, vals, [100, 1], @min, NaN), -%! accumarray (subs, vals, [100, 1], @(x) min (x), NaN)); - -%!test -%! subs = ceil (rand (2000, 2)*30); -%! subsc = num2cell (subs, 1); -%! vals = rand (2000, 1); -%! assert (accumarray (subsc, vals, [], [], 0, true), -%! accumarray (subs, vals, [], [], 0, true)); - -%!test -%! subs = ceil (rand (2000, 3)*10); -%! subsc = num2cell (subs, 1); -%! vals = rand (2000, 1); -%! assert (accumarray (subsc, vals, [], @max), -%! accumarray (subs, vals, [], @max)); - -%!error accumarray (1:5) -%!error accumarray ([1,2,3],1:2) - -## Handle empty arrays -%!test <*47287> -%! ## min, max, and sum are special cases within accumarray so test them. -%! funcs = {@(x) length (x) > 1, @min, @max, @sum}; -%! for idx = 1:numel (funcs) -%! assert (accumarray (zeros (0, 1), [], [0 1] , funcs{idx}), zeros (0, 1)); -%! assert (accumarray (zeros (0, 1), [], [1 0] , funcs{idx}), zeros (1, 0)); -%! assert (accumarray (zeros (0, 1), [], [] , funcs{idx}), zeros (0, 1)); -%! endfor - -## Matlab returns an array of doubles even though FCN returns cells. In -## Octave, we do not have that bug, at least for this case. -%!assert (accumarray (zeros (0, 1), [], [0 1] , @(x) {x}), cell (0, 1)) - -%!error -%! accumarray ([1; 2; 3], [1; 2; 3], [3 1], '@(x) {x}') From 58eb7de5bbdf99a69552663ec870df04a44eea71 Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Fri, 12 Jul 2024 12:41:57 +0300 Subject: [PATCH 14/21] Delete samples/Octave/var.m --- samples/Octave/var.m | 956 ------------------------------------------- 1 file changed, 956 deletions(-) delete mode 100644 samples/Octave/var.m diff --git a/samples/Octave/var.m b/samples/Octave/var.m deleted file mode 100644 index cde793acb0..0000000000 --- a/samples/Octave/var.m +++ /dev/null @@ -1,956 +0,0 @@ -######################################################################## -## -## Copyright (C) 1995-2024 The Octave Project Developers -## -## See the file COPYRIGHT.md in the top-level directory of this -## distribution or . -## -## This file is part of Octave. -## -## 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 3 of the License, or -## (at your option) any later version. -## -## Octave 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 -## . -## -######################################################################## - -## -*- texinfo -*- -## @deftypefn {} {@var{v} =} var (@var{x}) -## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}) -## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}, @var{dim}) -## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}, @var{vecdim}) -## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}, @qcode{"all"}) -## @deftypefnx {} {@var{v} =} var (@dots{}, @var{nanflag}) -## @deftypefnx {} {[@var{v}, @var{m}] =} var (@dots{}) -## Compute the variance of the elements of the vector @var{x}. -## -## The variance is defined as -## @tex -## $$ {\rm var}(x) = {1\over N-1} \sum_{i=1}^N (x_i - \bar x )^2 $$ -## where $\bar{x}$ is the mean value of @var{x} and $N$ is the number of -## elements of @var{x}. -## @end tex -## @ifnottex -## -## @example -## @group -## var (@var{x}) = (1 / (N-1)) * SUM_i ((@var{x}(i) - mean(@var{x}))^2) -## @end group -## @end example -## -## @noindent -## where @math{N} is the number of elements of @var{x}. -## @end ifnottex -## -## If @var{x} is an array, compute the variance along the first non-singleton -## dimensions of @var{x}. -## -## The optional argument @var{w} determines the weighting scheme to use. Valid -## values are: -## -## @table @asis -## @item 0 [default]: -## Normalize with @math{N-1} (population variance). This provides the square -## root of the best unbiased estimator of the variance. -## -## @item 1: -## Normalize with @math{N} (sample variance). This provides the square root of -## the second moment around the mean. -## -## @item a vector: -## Compute the weighted variance with non-negative weights. The length of -## @var{w} must equal the size of @var{x} in the operating dimension. NaN -## values are permitted in @var{w}, will be multiplied with the associated -## values in @var{x}, and can be excluded by the @var{nanflag} option. -## -## @item an array: -## Similar to vector weights, but @var{w} must be the same size as @var{x}. If -## the operating dimension is supplied as @var{vecdim} or @qcode{"all"} and -## @var{w} is not a scalar, @var{w} must be an same-sized array. -## @end table -## -## Note: @var{w} must always be specified before specifying any of the -## following dimension options. To use the default value for @var{w} you -## may pass an empty input argument []. -## -## The optional variable @var{dim} forces @code{var} to operate over the -## specified dimension, which must be a positive integer-valued number. -## Specifying any singleton dimension in @var{x}, including any dimension -## exceeding @code{ndims (@var{x})}, will result in a variance of 0. -## -## Specifying the dimensions as @var{vecdim}, a vector of non-repeating -## dimensions, will return the variance calculated over the array slice defined -## by @var{vecdim}. If @var{vecdim} indexes all dimensions of @var{x}, then it -## is equivalent to the option @qcode{"all"}. Any dimension in @var{vecdim} -## greater than @code{ndims (@var{x})} is ignored. -## -## Specifying the dimension as @qcode{"all"} will force @code{var} to -## operate on all elements of @var{x}, and is equivalent to @code{var -## (@var{x}(:))}. -## -## The optional variable @var{nanflag} specifies whether to include or exclude -## NaN values from the calculation using any of the previously specified input -## argument combinations. The default value for @var{nanflag} is -## @qcode{"includenan"} which keeps NaN values in the calculation. To -## exclude NaN values set the value of @var{nanflag} to @qcode{"omitnan"}. -## The output will still contain NaN values if @var{x} consists of all NaN -## values in the operating dimension. -## -## The optional second output variable @var{m} contains the mean of the -## elements of @var{x} used to calculate the variance. If @var{v} is the -## weighted variance, then @var{m} is also the weighted mean. -## -## @seealso{std, mean, cov, skewness, kurtosis, moment} -## @end deftypefn - -function [v, m] = var (x, varargin) - - if (nargin < 1 || nargin > 4) - print_usage (); - endif - - ## initialize variables - all_flag = false; - omitnan = false; - nvarg = numel (varargin); - varg_chars = cellfun ('ischar', varargin); - - ## Check all char arguments. - if (nvarg == 3 && ! varg_chars(3)) - print_usage (); - endif - - if (any (varg_chars)) - for argin = varargin(varg_chars) - switch (lower (argin{1})) - case "all" - all_flag = true; - case "omitnan" - omitnan = true; - case "includenan" - omitnan = false; - otherwise - print_usage (); - endswitch - endfor - varargin(varg_chars) = []; - nvarg = numel (varargin); - endif - - ## FIXME: When sparse can broadcast ops then remove sparse checks and hacks. - x_issparse = issparse (x); - w = 0; - weighted = false; # true if weight vector/array used - vecdim = []; - vecempty = true; - vecdim_scalar_vector = [false, false]; # [false, false] for empty vecdim - szx = size (x); - ndx = ndims (x); - - ## Check numeric arguments - if (! (isnumeric (x))) - error ("var: X must be a numeric vector or matrix"); - endif - if (isa (x, "single")) - outtype = "single"; - else - outtype = "double"; - endif - - if (nvarg > 0) - if (nvarg > 2 || any (! cellfun ('isnumeric', varargin))) - print_usage (); - endif - ## Process weight input - if (any (varargin{1} < 0)) - error ("var: weights must not contain any negative values"); - endif - if (isscalar (varargin{1})) - w = varargin{1}; - if (! (w == 0 || w == 1) && ! isscalar (x)) - error ("var: normalization scalar must be either 0 or 1"); - endif - elseif (numel (varargin{1}) > 1) - weights = varargin{1}; - weighted = true; - endif - if (nvarg > 1) - ## Process dimension input - vecdim = varargin{2}; - if (! (vecempty = isempty (vecdim))) - ## Check for empty vecdim, won't change vsv if nonzero size empty - vecdim_scalar_vector = [isscalar(vecdim), isvector(vecdim)]; - endif - if (! (vecdim_scalar_vector(2) && all (vecdim > 0)) - || any (rem (vecdim, 1))) - error ("var: DIM must be a positive integer scalar or vector"); - endif - if (vecdim_scalar_vector(1) && vecdim > ndx && ! isempty (x)) - ## Scalar dimension larger than ndims(x), variance of any single number - ## is zero, except for inf, NaN, and empty values of x. - v = zeros (szx, outtype); - vn = ! isfinite (x); - v(vn) = NaN; - m = x; - return; - endif - if (vecdim_scalar_vector == [0 1] && (! all (diff (sort (vecdim))))) - error ("var: VECDIM must contain non-repeating positive integers"); - endif - endif - endif - - ## Check for conflicting input arguments - if (all_flag && ! vecempty) - error ("var: 'all' flag cannot be used with DIM or VECDIM options"); - endif - if (weighted) - if (all_flag) - if (isvector (weights)) - if (numel (weights) != numel (x)) - error ("var: weight vector element count does not match X"); - endif - elseif (! (isequal (size (weights), szx))) - error ("var: weight matrix or array does not match X in size"); - endif - - elseif (vecempty) - ## Find the first non-singleton dimension. - (dim = find (szx > 1, 1)) || (dim = 1); - if (isvector (weights)) - if (numel (weights) != szx(dim)) - error ("var: weight vector length does not match operating dimension"); - endif - elseif (! isequal (size (weights), szx)) - error ("var: weight matrix or array does not match X in size"); - endif - elseif (vecdim_scalar_vector(1)) - if (isvector (weights)) - if (numel (weights) != szx(vecdim)) - error ("var: weight vector length does not match operating dimension"); - endif - elseif (! isequal (size (weights), szx)) - error ("var: weight matrix or array does not match X in size"); - endif - - elseif (vecdim_scalar_vector(2) && ! (isequal (size (weights), szx))) - error ("var: weight matrix or array does not match X in size"); - endif - endif - - ## Handle special cases of empty or scalar X and exit early. - if (isempty (x)) - if (vecempty && (ndx == 2 || all (szx == 0))) - v = NaN (outtype); - if (nargout > 1) - m = NaN (outtype); - endif - return; - endif - if (vecdim_scalar_vector(1)) - szx(vecdim) = 1; - v = NaN (szx, outtype); - if (nargout > 1) - m = NaN (szx, outtype); - endif - return; - endif - elseif (isscalar (x)) - if (isfinite (x)) - v = zeros (outtype); - else - v = NaN (outtype); - endif - if (nargout > 1) - m = x; - endif - return; - endif - - if (nvarg == 0) - ## Only numeric input argument, no dimensions or weights. - if (all_flag) - x = x(:); - if (omitnan) - x = x(! isnan (x)); - endif - n = length (x); - m = sum (x) ./ n; - v = sum (abs (x - m) .^ 2) ./ (n - 1 + w); - if (n == 1) - v = 0; - endif - else - ## Find the first non-singleton dimension. - (dim = find (szx > 1, 1)) || (dim = 1); - n = szx(dim); - if (omitnan) - n = sum (! isnan (x), dim); - xn = isnan (x); - x(xn) = 0; - endif - m = sum (x, dim) ./ n; - dims = ones (1, ndx); - dims(dim) = szx(dim); - if (x_issparse) - m_exp = repmat (m, dims); - else - m_exp = m .* ones (dims); - endif - if (omitnan) - x(xn) = m_exp(xn); - endif - v = sumsq (x - m_exp, dim) ./ (n - 1 + w); - if (numel (n) == 1) - divby0 = (n .* ones (size (v))) == 1; - else - divby0 = n == 1; - endif - v(divby0) = 0; - endif - - elseif (nvarg == 1) - ## Two numeric input arguments, w or weights given. - if (all_flag) - x = x(:); - if (weighted) - weights = weights(:); - wx = weights .* x; - else - weights = ones (length (x), 1); - wx = x; - endif - - if (omitnan) - xn = isnan (wx); - wx = wx(! xn); - weights = weights(! xn); - x = x(! xn); - endif - n = length (wx); - m = sum (wx) ./ sum (weights); - if (weighted) - v = sum (weights .* (abs (x - m) .^ 2)) ./ sum (weights); - else - v = sum (weights .* (abs (x - m) .^ 2)) ./ (n - 1 + w); - if (n == 1) - v = 0; - endif - endif - - else - ## Find the first non-singleton dimension. - (dim = find (szx > 1, 1)) || (dim = 1); - if (! weighted) - weights = ones (szx); - wx = x; - else - if (isvector (weights)) - dims = 1:ndx; - dims([1, dim]) = [dim, 1]; - weights = zeros (szx) + permute (weights(:), dims); - endif - wx = weights .* x; - endif - n = size (wx, dim); - if (omitnan) - xn = isnan (wx); - n = sum (! xn, dim); - wx(xn) = 0; - weights(xn) = 0; - endif - m = sum (wx, dim) ./ sum (weights, dim); - dims = ones (1, ndims (wx)); - dims(dim) = size (wx, dim); - if (x_issparse) - m_exp = repmat (m, dims); - else - m_exp = m .* ones (dims); - endif - if (omitnan) - x(xn) = m_exp(xn); - endif - if (weighted) - v = sum (weights .* ((x - m_exp) .^ 2), dim) ./ sum (weights, dim); - else - v = sumsq (x - m_exp, dim) ./ (n - 1 + w); - if (numel (n) == 1) - divby0 = (n .* ones (size (v))) == 1; - else - divby0 = n == 1; - endif - v(divby0) = 0; - endif - endif - - elseif (nvarg == 2) - ## Three numeric input arguments, both w or weights and dim or vecdim given. - if (vecdim_scalar_vector(1)) - if (! weighted) - weights = ones (szx); - wx = x; - else - if (isvector (weights)) - dims = 1:ndx; - dims([1, vecdim]) = [vecdim, 1]; - weights = zeros (szx) + permute (weights(:), dims); - endif - wx = weights .* x; - endif - n = size (wx, vecdim); - if (omitnan) - n = sum (! isnan (wx), vecdim); - xn = isnan (wx); - wx(xn) = 0; - weights(xn) = 0; - endif - m = sum (wx, vecdim) ./ sum (weights, vecdim); - dims = ones (1, ndims (wx)); - dims(vecdim) = size (wx, vecdim); - if (x_issparse) - m_exp = repmat (m, dims); - else - m_exp = m .* ones (dims); - endif - if (omitnan) - x(xn) = m_exp(xn); - endif - if (weighted) - v = sum (weights .* ((x - m_exp) .^ 2), vecdim) ... - ./ sum (weights, vecdim); - else - v = sumsq (x - m_exp, vecdim); - vn = isnan (v); - v = v ./ (n - 1 + w); - if (numel (n) == 1) - divby0 = (n .* ones (size (v))) == 1; - else - divby0 = n == 1; - endif - v(divby0) = 0; - v(vn) = NaN; - endif - - else - ## Weights and nonscalar vecdim specified - - ## Ignore dimensions in VECDIM larger than actual array. - remdims = 1 : ndx; # All dimensions - vecdim(vecdim > ndx) = []; - ## Calculate permutation vector - remdims(vecdim) = []; # Delete dimensions specified by vecdim - nremd = numel (remdims); - - ## If all dimensions are given, it is equivalent to the 'all' flag. - if (nremd == 0) - x = x(:); - if (weighted) - weights = weights(:); - wx = weights .* x; - else - weights = ones (length (x), 1); - wx = x; - endif - - if (omitnan) - xn = isnan (wx); - wx = wx(! xn); - weights = weights(! xn); - x = x(! xn); - endif - n = length (wx); - m = sum (wx) ./ sum (weights); - if (weighted) - v = sum (weights .* (abs (x - m) .^ 2)) ./ sum (weights); - else - v = sum (weights .* (abs (x - m) .^ 2)) ./ (n - 1 + w); - if (n == 1) - v = 0; - endif - endif - - else - - ## FIXME: much of the reshaping can be skipped once Octave's sum can - ## take a vecdim argument. - - ## Apply weights - if (weighted) - wx = weights .* x; - else - weights = ones (szx); - wx = x; - endif - - ## Permute to push vecdims to back - perm = [remdims, vecdim]; - wx = permute (wx, perm); - weights = permute (weights, perm); - x = permute (x, perm); - - ## Reshape to squash all vecdims in final dimension - szwx = size (wx); - sznew = [szwx(1:nremd), prod(szwx(nremd+1:end))]; - wx = reshape (wx, sznew); - weights = reshape (weights, sznew); - x = reshape (x, sznew); - - ## Calculate var on final squashed dimension - dim = nremd + 1; - n = size (wx, dim); - if (omitnan) - xn = isnan (wx); - n = sum (! xn, dim); - wx(xn) = 0; - weights(xn) = 0; - endif - m = sum (wx, dim) ./ sum (weights, dim); - m_exp = zeros (sznew) + m; - if (omitnan) - x(xn) = m_exp(xn); - endif - if (weighted) - v = sum (weights .* ((x - m_exp) .^ 2), dim) ./ sum (weights, dim); - else - v = sumsq (x - m_exp, dim) ./ (n - 1 + w); - if (numel (n) == 1) - divby0 = n .* ones (size (v)) == 1; - else - divby0 = n == 1; - endif - v(divby0) = 0; - endif - - ## Inverse permute back to correct dimensions - v = ipermute (v, perm); - if (nargout > 1) - m = ipermute (m, perm); - endif - endif - endif - endif - - ## Preserve class type - if (nargout < 2) - if (strcmp (outtype, "single")) - v = single (v); - else - v = double (v); - endif - else - if (strcmp (outtype, "single")) - v = single (v); - m = single (m); - else - v = double (v); - m = double (m); - endif - endif - -endfunction - - -%!assert (var (13), 0) -%!assert (var (single (13)), single (0)) -%!assert (var ([1,2,3]), 1) -%!assert (var ([1,2,3], 1), 2/3, eps) -%!assert (var ([1,2,3], [], 1), [0,0,0]) -%!assert (var ([1,2,3], [], 3), [0,0,0]) -%!assert (var (5, 99), 0) -%!assert (var (5, 99, 1), 0) -%!assert (var (5, 99, 2), 0) -%!assert (var ([5 3], [99 99], 2), 1) -%!assert (var ([1:7], [1:7]), 3) -%!assert (var ([eye(3)], [1:3]), [5/36, 2/9, 1/4], eps) -%!assert (var (ones (2,2,2), [1:2], 3), [(zeros (2,2))]) -%!assert (var ([1 2; 3 4], 0, 'all'), var ([1:4])) -%!assert (var (reshape ([1:8], 2, 2, 2), 0, [1 3]), [17/3 17/3], eps) -%!assert (var ([1 2 3;1 2 3], [], [1 2]), 0.8, eps) - -## Test single input and optional arguments "all", DIM, "omitnan") -%!test -%! x = [-10:10]; -%! y = [x;x+5;x-5]; -%! assert (var (x), 38.5); -%! assert (var (y, [], 2), [38.5; 38.5; 38.5]); -%! assert (var (y, 0, 2), [38.5; 38.5; 38.5]); -%! assert (var (y, 1, 2), ones (3,1) * 36.66666666666666, 1e-14); -%! assert (var (y, "all"), 54.19354838709678, 1e-14); -%! y(2,4) = NaN; -%! assert (var (y, "all"), NaN); -%! assert (var (y, "all", "includenan"), NaN); -%! assert (var (y, "all", "omitnan"), 55.01533580116342, 1e-14); -%! assert (var (y, 0, 2, "includenan"), [38.5; NaN; 38.5]); -%! assert (var (y, [], 2), [38.5; NaN; 38.5]); -%! assert (var (y, [], 2, "omitnan"), [38.5; 37.81842105263158; 38.5], 1e-14); - -## Tests for different weight and omitnan code paths -%!assert (var ([1 NaN 3], [1 2 3], "omitnan"), 0.75, eps) -%!assert (var ([1 2 3], [1 NaN 3], "omitnan"), 0.75, eps) -%!assert (var (magic(3), [1 NaN 3], "omitnan"), [3 12 3], eps) -%!assert (var ([1 NaN 3], [1 2 3], "omitnan", "all"), 0.75, eps) -%!assert (var ([1 NaN 3], [1 2 3], "all", "omitnan"), 0.75, eps) -%!assert (var ([1 2 3], [1 NaN 3], "omitnan", "all"), 0.75, eps) -%!assert (var ([1 NaN 3], [1 2 3], 2, "omitnan"), 0.75, eps) -%!assert (var ([1 2 3], [1 NaN 3], 2, "omitnan"), 0.75, eps) -%!assert (var (magic(3), [1 NaN 3], 1, "omitnan"), [3 12 3], eps) -%!assert (var (magic(3), [1 NaN 3], 2, "omitnan"), [0.75;3;0.75], eps) -%!assert (var ([4 4; 4 6; 6 6], [1 3], 2, 'omitnan'), [0;0.75;0], eps) -%!assert (var ([4 NaN; 4 6; 6 6], [1 2 3], 1, 'omitnan'), [1 0]) -%!assert (var ([4 NaN; 4 6; 6 6], [1 3], 2, 'omitnan'), [0;0.75;0], eps) -%!assert (var (3*reshape(1:18, [3 3 2]), [1 2 3], 1, 'omitnan'), ones(1,3,2)*5) -%!assert (var (reshape(1:18, [3 3 2]), [1 2 3], 2, 'omitnan'), 5*ones(3,1,2)) -%!assert (var (3*reshape(1:18, [3 3 2]), ones (3,3,2), [1 2], 'omitnan'), ... -%! 60 * ones(1,1,2)) -%!assert (var (3*reshape(1:18, [3 3 2]), ones (3,3,2), [1 4], 'omitnan'), ... -%! 6 * ones(1,3,2)) -%!assert (var (6*reshape(1:18, [3 3 2]), ones (3,3,2), [1:3], 'omitnan'), 969) -%!test -%! x = reshape(1:18, [3 3 2]); -%! x([2, 14]) = NaN; -%! w = ones (3,3,2); -%! assert (var (16*x, w, [1:3], 'omitnan'), 6519); -%!test -%! x = reshape(1:18, [3 3 2]); -%! w = ones (3,3,2); -%! w([2, 14]) = NaN; -%! assert (var (16*x, w, [1:3], 'omitnan'), 6519); - -## Test input case insensitivity -%!assert (var ([1 2 3], "aLl"), 1) -%!assert (var ([1 2 3], "OmitNan"), 1) -%!assert (var ([1 2 3], "IncludeNan"), 1) - -## Test dimension indexing with vecdim in n-dimensional arrays -%!test -%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); -%! assert (size (var (x, 0, [3 2])), [10, 1, 1, 3]); -%! assert (size (var (x, 1, [1 2])), [1, 1, 6, 3]); -%! assert (size (var (x, [], [1 2 4])), [1, 1, 6]); -%! assert (size (var (x, 0, [1 4 3])), [1, 40]); -%! assert (size (var (x, [], [1 2 3 4])), [1, 1]); - -## Test matrix with vecdim, weighted, matrix weights, omitnan -%!assert (var (3*magic(3)), [63 144 63]) -%!assert (var (3*magic(3), 'omitnan'), [63 144 63]) -%!assert (var (3*magic(3), 1), [42 96 42]) -%!assert (var (3*magic(3), 1, 'omitnan'), [42 96 42]) -%!assert (var (3*magic(3), ones(1,3), 1), [42 96 42]) -%!assert (var (3*magic(3), ones(1,3), 1, 'omitnan'), [42 96 42]) -%!assert (var (2*magic(3), [1 1 NaN], 1, 'omitnan'), [25 16 1]) -%!assert (var (3*magic(3), ones(3,3)), [42 96 42]) -%!assert (var (3*magic(3), ones(3,3), 'omitnan'), [42 96 42]) -%!assert (var (3*magic(3), [1 1 1; 1 1 1; 1 NaN 1], 'omitnan'), [42 36 42]) -%!assert (var (3*magic(3), ones(3,3), 1), [42 96 42]) -%!assert (var (3*magic(3), ones(3,3), 1, 'omitnan'), [42 96 42]) -%!assert (var (3*magic(3), [1 1 1; 1 1 1; 1 NaN 1], 1, 'omitnan'), [42 36 42]) -%!assert (var (3*magic(3), ones(3,3), [1 4]), [42 96 42]) -%!assert (var (3*magic(3), ones(3,3), [1 4], 'omitnan'), [42 96 42]) -%!assert (var (3*magic(3), [1 1 1; 1 1 1; 1 NaN 1],[1 4],'omitnan'), [42 36 42]) - -## Test results with vecdim in n-dimensional arrays and "omitnan" -%!test -%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); -%! v = repmat (33.38912133891213, [10, 1, 1, 3]); -%! assert (var (x, 0, [3, 2]), v, 1e-14); -%! v = repmat (33.250, [10, 1, 1, 3]); -%! assert (var (x, 1, [3, 2]), v, 1e-14); -%! x(2,5,6,3) = NaN; -%! v(2,1,1,3) = NaN; -%! assert (var (x, 1, [3, 2]), v, 4e-14); -%! v = repmat (33.38912133891213, [10 1 1 3]); -%! v(2,1,1,3) = NaN; -%! assert (var (x, [], [3, 2]), v, 4e-14); -%! v(2,1,1,3) = 33.40177912169048; -%! assert (var (x, [], [3, 2], "omitnan"), v, 4e-14); - -## Testing weights vector & arrays -%!assert (var (ones (2,2,2), [1:2], 3), [(zeros (2, 2))]) -%!assert (var (magic (3), [1:9], "all"), 6.666666666666667, 1e-14) - -## Test exceeding dimensions -%!assert (var (ones (2,2), [], 3), zeros (2,2)) -%!assert (var (ones (2,2,2), [], 99), zeros (2,2,2)) -%!assert (var (magic (3), [], 3), zeros (3,3)) -%!assert (var (magic (3), [], 1), [7, 16, 7]) -%!assert (var (magic (3), [], [1 3]), [7, 16, 7]) -%!assert (var (magic (3), [], [1 99]), [7, 16, 7]) - -## Test empty inputs -%!assert (var ([]), NaN) -%!assert (class (var (single ([]))), "single") -%!assert (var ([],[],1), NaN(1,0)) -%!assert (var ([],[],2), NaN(0,1)) -%!assert (var ([],[],3), []) -%!assert (class (var (single ([]), [], 1)), "single") -%!assert (var (ones (1,0)), NaN) -%!assert (var (ones (1,0), [], 1), NaN(1,0)) -%!assert (var (ones (1,0), [], 2), NaN) -%!assert (var (ones (1,0), [], 3), NaN(1,0)) -%!assert (class (var (ones (1, 0, "single"), [], 1)), "single") -%!assert (var (ones (0,1)), NaN) -%!assert (var (ones (0,1), [], 1), NaN) -%!assert (var (ones (0,1), [], 2), NaN(0,1)) -%!assert (var (ones (0,1), [], 3), NaN(0,1)) -%!assert (var (ones (1,3,0,2)), NaN(1,1,0,2)) -%!assert (var (ones (1,3,0,2), [], 1), NaN(1,3,0,2)) -%!assert (var (ones (1,3,0,2), [], 2), NaN(1,1,0,2)) -%!assert (var (ones (1,3,0,2), [], 3), NaN(1,3,1,2)) -%!assert (var (ones (1,3,0,2), [], 4), NaN(1,3,0)) -%!test -%! [~, m] = var ([]); -%! assert (m, NaN); - -## Test optional mean output -%!test <*62395> -%! [~, m] = var (13); -%! assert (m, 13); -%! [~, m] = var (single(13)); -%! assert (m, single(13)); -%! [~, m] = var ([1, 2, 3; 3 2 1], []); -%! assert (m, [2 2 2]); -%! [~, m] = var ([1, 2, 3; 3 2 1], [], 1); -%! assert (m, [2 2 2]); -%! [~, m] = var ([1, 2, 3; 3 2 1], [], 2); -%! assert (m, [2 2]'); -%! [~, m] = var ([1, 2, 3; 3 2 1], [], 3); -%! assert (m, [1 2 3; 3 2 1]); - -## Test mean output, weighted inputs, vector dims -%!test <*62395> -%! [~, m] = var (5,99); -%! assert (m, 5); -%! [~, m] = var ([1:7], [1:7]); -%! assert (m, 5); -%! [~, m] = var ([eye(3)], [1:3]); -%! assert (m, [1/6, 1/3, 0.5], eps); -%! [~, m] = var (ones (2,2,2), [1:2], 3); -%! assert (m, ones (2,2)); -%! [~, m] = var ([1 2; 3 4], 0, 'all'); -%! assert (m, 2.5, eps); -%! [~, m] = var (reshape ([1:8], 2, 2, 2), 0, [1 3]); -%! assert (m, [3.5, 5.5], eps); -%!test -%! [v, m] = var (4 * eye (2), [1, 3]); -%! assert (v, [3, 3]); -%! assert (m, [1, 3]); - -## Test mean output, empty inputs, omitnan -%!test <*62395> -%! [~, m] = var ([]); -%! assert (m, NaN); -#%! [~, m] = var ([],[],1); -#%! assert (m, NaN(1,0)); -#%! [~, m] = var ([],[],2); -#%! assert (m, NaN(0,1)); -#%! [~, m] = var ([],[],3); -#%! assert (m, []); -#%! [~, m] = var (ones (1,3,0,2)); -#%! assert (m, NaN(1,1,0,2)); - -## Test mean output, nD array -%!test <*62395> -%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); -%! [~, m] = var (x, 0, [3 2]); -%! assert (m, mean (x, [3 2])); -%! [~, m] = var (x, 0, [1 2]); -%! assert (m, mean (x, [1 2])); -%! [~, m] = var (x, 0, [1 3 4]); -%! assert (m, mean (x, [1 3 4])); -%!test -%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); -%! x(2,5,6,3) = NaN; -%! [~, m] = var (x, 0, [3 2], "omitnan"); -%! assert (m, mean (x, [3 2], "omitnan")); - -## Test Inf and NaN inputs -%!test <*63203> -%! [v, m] = var (Inf); -%! assert (v, NaN); -%! assert (m, Inf); -%!test <*63203> -%! [v, m] = var (NaN); -%! assert (v, NaN); -%! assert (m, NaN); -%!test <*63203> -%! [v, m] = var ([1, Inf, 3]); -%! assert (v, NaN); -%! assert (m, Inf); -%!test <*63203> -%! [v, m] = var ([1, Inf, 3]'); -%! assert (v, NaN); -%! assert (m, Inf); -%!test <*63203> -%! [v, m] = var ([1, NaN, 3]); -%! assert (v, NaN); -%! assert (m, NaN); -%!test <*63203> -%! [v, m] = var ([1, NaN, 3]'); -%! assert (v, NaN); -%! assert (m, NaN); -%!test <*63203> -%! [v, m] = var ([1, Inf, 3], [], 1); -%! assert (v, [0, NaN, 0]); -%! assert (m, [1, Inf, 3]); -%!test <*63203> -%! [v, m] = var ([1, Inf, 3], [], 2); -%! assert (v, NaN); -%! assert (m, Inf); -%!test <*63203> -%! [v, m] = var ([1, Inf, 3], [], 3); -%! assert (v, [0, NaN, 0]); -%! assert (m, [1, Inf, 3]); -%!test <*63203> -%! [v, m] = var ([1, NaN, 3], [], 1); -%! assert (v, [0, NaN, 0]); -%! assert (m, [1, NaN, 3]); -%!test <*63203> -%! [v, m] = var ([1, NaN, 3], [], 2); -%! assert (v, NaN); -%! assert (m, NaN); -%!test <*63203> -%! [v, m] = var ([1, NaN, 3], [], 3); -%! assert (v, [0, NaN, 0]); -%! assert (m, [1, NaN, 3]); -%!test <*63203> -%! [v, m] = var ([1, 2, 3; 3, Inf, 5]); -%! assert (v, [2, NaN, 2]); -%! assert (m, [2, Inf, 4]); -%!test <*63203> -%! [v, m] = var ([1, Inf, 3; 3, Inf, 5]); -%! assert (v, [2, NaN, 2]); -%! assert (m, [2, Inf, 4]); -%!test <*63203> -%! [v, m] = var ([1, 2, 3; 3, NaN, 5]); -%! assert (v, [2, NaN, 2]); -%! assert (m, [2, NaN, 4]); -%!test <*63203> -%! [v, m] = var ([1, NaN, 3; 3, NaN, 5]); -%! assert (v, [2, NaN, 2]); -%! assert (m, [2, NaN, 4]); -%!test <*63203> -%! [v, m] = var ([Inf, 2, NaN]); -%! assert (v, NaN); -%! assert (m, NaN); -%!test <*63203> -%! [v, m] = var ([Inf, 2, NaN]'); -%! assert (v, NaN); -%! assert (m, NaN); -%!test <*63203> -%! [v, m] = var ([NaN, 2, Inf]); -%! assert (v, NaN); -%! assert (m, NaN); -%!test <*63203> -%! [v, m] = var ([NaN, 2, Inf]'); -%! assert (v, NaN); -%! assert (m, NaN); -%!test <*63203> -%! [v, m] = var ([Inf, 2, NaN], [], 1); -%! assert (v, [NaN, 0, NaN]); -%! assert (m, [Inf, 2, NaN]); -%!test <*63203> -%! [v, m] = var ([Inf, 2, NaN], [], 2); -%! assert (v, NaN); -%! assert (m, NaN); -%!test <*63203> -%! [v, m] = var ([NaN, 2, Inf], [], 1); -%! assert (v, [NaN, 0, NaN]); -%! assert (m, [NaN, 2, Inf]); -%!test <*63203> -%! [v, m] = var ([NaN, 2, Inf], [], 2); -%! assert (v, NaN); -%! assert (m, NaN); -%!test <*63203> -%! [v, m] = var ([1, 3, NaN; 3, 5, Inf]); -%! assert (v, [2, 2, NaN]); -%! assert (m, [2, 4, NaN]); -%!test <*63203> -%! [v, m] = var ([1, 3, Inf; 3, 5, NaN]); -%! assert (v, [2, 2, NaN]); -%! assert (m, [2, 4, NaN]); - -## Test sparse/diagonal inputs -%!test <*63291> -%! [v, m] = var (2 * eye (2)); -%! assert (v, [2, 2]); -%! assert (m, [1, 1]); -%!test <*63291> -%! [v, m] = var (4 * eye (2), [1, 3]); -%! assert (v, [3, 3]); -%! assert (m, [1, 3]); -%!test <*63291> -%! [v, m] = var (sparse (2 * eye (2))); -%! assert (full (v), [2, 2]); -%! assert (full (m), [1, 1]); -%!test <*63291> -%! [v, m] = var (sparse (4 * eye (2)), [1, 3]); -%! assert (full (v), [3, 3]); -%! assert (full (m), [1, 3]); -%!test<*63291> -%! [v, m] = var (sparse (eye (2))); -%! assert (issparse (v)); -%! assert (issparse (m)); -%!test<*63291> -%! [v, m] = var (sparse (eye (2)), [1, 3]); -%! assert (issparse (v)); -%! assert (issparse (m)); - -## Test input validation -%!error var () -%!error var (1, 2, "omitnan", 3) -%!error var (1, 2, 3, 4) -%!error var (1, 2, 3, 4, 5) -%!error var (1, "foo") -%!error var (1, [], "foo") -%!error var ([1 2 3], 2) -%!error var ([1 2], 2, "all") -%!error var ([1 2],0.5, "all") -%!error var (1, -1) -%!error var (1, [1 -1]) -%!error ... -%! var ([1 2 3], [1 -1 0]) -%!error var ({1:5}) -%!error var ("char") -%!error var (['A'; 'B']) -%!error var (1, [], ones (2,2)) -%!error var (1, 0, 1.5) -%!error var (1, [], 0) -%!error var (1, [], 1.5) -%!error var ([1 2 3], [], [-1 1]) -%!error ... -%! var (repmat ([1:20;6:25], [5 2 6 3]), 0, [1 2 2 2]) -%!error ... -%! var ([1 2], eye (2)) -%!error ... -%! var ([1 2 3 4], [1 2; 3 4]) -%!error ... -%! var ([1 2 3 4], [1 2; 3 4], 1) -%!error ... -%! var ([1 2 3 4], [1 2; 3 4], [2 3]) -%!error ... -%! var (ones (2, 2), [1 2], [1 2]) -%!error ... -%! var ([1 2 3 4; 5 6 7 8], [1 2 1 2 1; 1 2 1 2 1], 1) -%!error ... -%! var (repmat ([1:20;6:25], [5 2 6 3]), repmat ([1:20;6:25], [5 2 3]), [2 3]) -%!error var ([1 2 3; 2 3 4], [1 3 4]) -%!error var ([1 2], [1 2 3]) -%!error var (1, [1 2]) -%!error var ([1 2 3; 2 3 4], [1 3 4], 1) -%!error var ([1 2 3; 2 3 4], [1 3], 2) -%!error var ([1 2], [1 2], 1) -%!error <'all' flag cannot be used with DIM or VECDIM options> ... -%! var (1, [], 1, "all") -%!error ... -%! var ([1 2 3; 2 3 4], [1 3], "all") -%!error ... -%! var (repmat ([1:20;6:25], [5 2 6 3]), repmat ([1:20;6:25], [5 2 3]), "all") From 29ce7d454fd38b23352f69b138d7e4b78b0faeb8 Mon Sep 17 00:00:00 2001 From: Andreas Bertsatos Date: Fri, 12 Jul 2024 12:42:45 +0300 Subject: [PATCH 15/21] Update test_heuristics.rb --- test/test_heuristics.rb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_heuristics.rb b/test/test_heuristics.rb index bfb138e097..41c879ee66 100755 --- a/test/test_heuristics.rb +++ b/test/test_heuristics.rb @@ -633,7 +633,7 @@ def test_m_by_heuristics "M" => all_fixtures("M", "MDB.m"), "Mathematica" => all_fixtures("Mathematica", "*.m") - all_fixtures("Mathematica", "Problem12.m"), "MATLAB" => all_fixtures("MATLAB", "create_ieee_paper_plots.m"), - "Octave" => all_fixtures("Octave", "accumarray.m"), + "Octave" => all_fixtures("Octave", "*.m"), "Limbo" => all_fixtures("Limbo", "*.m"), nil => ambiguous }) From b15123b669607d57225f2a9672c5a5f19a531958 Mon Sep 17 00:00:00 2001 From: pr0m1th3as Date: Fri, 12 Jul 2024 15:38:58 +0300 Subject: [PATCH 16/21] add the cached license file for Octave --- .../Octave-language-grammar.dep.yml | 31 +++++++++++++++++++ 1 file changed, 31 insertions(+) create mode 100644 vendor/licenses/git_submodule/Octave-language-grammar.dep.yml diff --git a/vendor/licenses/git_submodule/Octave-language-grammar.dep.yml b/vendor/licenses/git_submodule/Octave-language-grammar.dep.yml new file mode 100644 index 0000000000..fe96bdce46 --- /dev/null +++ b/vendor/licenses/git_submodule/Octave-language-grammar.dep.yml @@ -0,0 +1,31 @@ +--- +name: Octave-language-grammar +version: 0f9f17f60e4668c7aaaee6d9a6ad59ce6e773908 +type: git_submodule +homepage: https://github.com/gnu-octave/Octave-language-grammar +license: mit +licenses: +- sources: LICENSE + text: | + MIT License + + Copyright (C) 2024 The Octave Project Developers + + Permission is hereby granted, free of charge, to any person obtaining a copy + of this software and associated documentation files (the "Software"), to deal + in the Software without restriction, including without limitation the rights + to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + copies of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included in all + copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +notices: [] From be1b7e132a569837cf554a43ae94f49ff78638bb Mon Sep 17 00:00:00 2001 From: pr0m1th3as Date: Fri, 12 Jul 2024 15:41:42 +0300 Subject: [PATCH 17/21] update Octave language_id as generated by script/update-ids script --- lib/linguist/languages.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/lib/linguist/languages.yml b/lib/linguist/languages.yml index 962d75754c..0d2df2834b 100644 --- a/lib/linguist/languages.yml +++ b/lib/linguist/languages.yml @@ -943,11 +943,11 @@ Caddyfile: type: data color: "#22b638" aliases: - - Caddy + - Caddy extensions: - - ".caddyfile" + - ".caddyfile" filenames: - - Caddyfile + - Caddyfile ace_mode: text tm_scope: source.Caddyfile language_id: 615465151 @@ -4887,7 +4887,7 @@ Octave: - ".m" tm_scope: source.octave ace_mode: octave - language_id: 312 + language_id: 854759768 Odin: type: programming color: "#60AFFE" @@ -5363,7 +5363,7 @@ Pkl: extensions: - ".pkl" interpreters: - - pkl + - pkl tm_scope: source.pkl ace_mode: text language_id: 288822799 From 9dd005d75dc3880803b4dc03eb045077a89d8330 Mon Sep 17 00:00:00 2001 From: pr0m1th3as Date: Fri, 12 Jul 2024 19:25:40 +0300 Subject: [PATCH 18/21] add Octave to grammars.yml --- grammars.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/grammars.yml b/grammars.yml index 23c95ad2f7..2f123c29b3 100644 --- a/grammars.yml +++ b/grammars.yml @@ -87,6 +87,8 @@ vendor/grammars/ObjectScript.tmBundle: - source.objectscript_class - source.objectscript_csp - source.objectscript_macros +vendor/grammars/Octave-language-grammar: +- source.octave vendor/grammars/PHP-Twig.tmbundle: - text.html.twig vendor/grammars/PogoScript.tmbundle: From 303b6b94bf89d36c3d4a7b0d390626a65478da54 Mon Sep 17 00:00:00 2001 From: pr0m1th3as Date: Fri, 12 Jul 2024 20:30:59 +0300 Subject: [PATCH 19/21] change Octave's ace_mode to text --- lib/linguist/languages.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/linguist/languages.yml b/lib/linguist/languages.yml index 0d2df2834b..48408f3849 100644 --- a/lib/linguist/languages.yml +++ b/lib/linguist/languages.yml @@ -4886,7 +4886,7 @@ Octave: extensions: - ".m" tm_scope: source.octave - ace_mode: octave + ace_mode: text language_id: 854759768 Odin: type: programming From 709c48a6e7be561cae001c534729f2eda14dae46 Mon Sep 17 00:00:00 2001 From: pr0m1th3as Date: Tue, 23 Jul 2024 16:37:01 +0300 Subject: [PATCH 20/21] add codemirror code and mime_type for MATLAB and Octave in the languages.yml --- lib/linguist/languages.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lib/linguist/languages.yml b/lib/linguist/languages.yml index 48408f3849..bf2f0db356 100644 --- a/lib/linguist/languages.yml +++ b/lib/linguist/languages.yml @@ -4025,6 +4025,8 @@ MATLAB: - ".m" tm_scope: source.matlab ace_mode: matlab + codemirror_code: octave + codemirror_mime_type: text/x-octave language_id: 225 MAXScript: type: programming @@ -4887,6 +4889,8 @@ Octave: - ".m" tm_scope: source.octave ace_mode: text + codemirror_mode: octave + codemirror_mime_type: text/x-octave language_id: 854759768 Odin: type: programming From 485c9783f017cda93c35baacd86658a1fb56d915 Mon Sep 17 00:00:00 2001 From: Colin Seymour Date: Mon, 5 Aug 2024 18:07:29 +0100 Subject: [PATCH 21/21] Update lib/linguist/languages.yml --- lib/linguist/languages.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/linguist/languages.yml b/lib/linguist/languages.yml index bf2f0db356..db6949e46c 100644 --- a/lib/linguist/languages.yml +++ b/lib/linguist/languages.yml @@ -4025,7 +4025,7 @@ MATLAB: - ".m" tm_scope: source.matlab ace_mode: matlab - codemirror_code: octave + codemirror_mode: octave codemirror_mime_type: text/x-octave language_id: 225 MAXScript: