Mercurial > octave-nkf
changeset 3191:e4f4b2d26ee9
[project @ 1998-10-23 05:43:59 by jwe]
line wrap: on
line diff
--- a/scripts/ChangeLog Mon Oct 19 22:37:57 1998 +0000 +++ b/scripts/ChangeLog Fri Oct 23 05:44:01 1998 +0000 @@ -1,3 +1,230 @@ +Fri Oct 23 00:21:55 1998 John W. Eaton <jwe@bevo.che.wisc.edu> + + * configure.in: Add finance/Makefile, statistics/base/Makefile, + statistics/distributions/Makefile, statistics/models/Makefile, + and statistics/tests/Makefile to the list of files to create. + + * finance/Makefile.in, statistics/base/Makefile.in, + statistics/distributions/Makefile.in, statistics/models/Makefile.in, + statistics/tests/Makefile.in: New files. + * statistics/Makefile.in: Delete file lists. Now only handle + subdirectories. + * Makefile.in (SUBDIRS): Add finance. + + * Move the following files from statistics to statistics/base: + + corrcoef.m + kurtosis.m + mahalanobis.m + median.m + ols.m + skewness.m + std.m + + New files, from Kurt Hornik's octave-ci package: + + * finance (new directory): + + fv.m + fvl.m + irr.m + nper.m + npv.m + pmt.m + pv.m + pvl.m + rate.m + vol.m + + * linear-algebra/dmult.m + + * signal: + + arch_fit.m + arch_rnd.m + arch_test.m + arma_rnd.m + autocor.m + autocov.m + autoreg_matrix.m + bartlett.m + blackman.m + diffpara.m + durbinlevinson.m + fractdiff.m + hamming.m + hanning.m + hurst.m + periodogram.m + rectangle_lw.m + rectangle_sw.m + sinetone.m + sinewave.m + spectral_adf.m + spectral_xdf.m + spencer.m + stft.m + synthesis.m + triangle_lw.m + triangle_sw.m + yulewalker.m + + * statistics/base (new directory): + + center.m + cloglog.m + cor.m + cov.m + cut.m + iqr.m + kendall.m + logit.m + mean.m + meansq.m + moment.m + ppplot.m + probit.m + qqplot.m + range.m + ranks.m + run_count.m + spearman.m + statistics.m + studentize.m + table.m + values.m + var.m + + (Replaces cov.m and mean.m with new versions.) + + * statistics/distributions (new directory): + + beta_cdf.m + beta_inv.m + beta_pdf.m + beta_rnd.m + binomial_cdf.m + binomial_inv.m + binomial_pdf.m + binomial_rnd.m + cauchy_cdf.m + cauchy_inv.m + cauchy_pdf.m + cauchy_rnd.m + chisquare_cdf.m + chisquare_inv.m + chisquare_pdf.m + chisquare_rnd.m + discrete_cdf.m + discrete_inv.m + discrete_pdf.m + discrete_rnd.m + empirical_cdf.m + empirical_inv.m + empirical_pdf.m + empirical_rnd.m + exponential_cdf.m + exponential_inv.m + exponential_pdf.m + exponential_rnd.m + f_cdf.m + f_inv.m + f_pdf.m + f_rnd.m + gamma_cdf.m + gamma_inv.m + gamma_pdf.m + gamma_rnd.m + geometric_cdf.m + geometric_inv.m + geometric_pdf.m + geometric_rnd.m + hypergeometric_cdf.m + hypergeometric_inv.m + hypergeometric_pdf.m + hypergeometric_rnd.m + kolmogorov_smirnov_cdf.m + laplace_cdf.m + laplace_inv.m + laplace_pdf.m + laplace_rnd.m + logistic_cdf.m + logistic_inv.m + logistic_pdf.m + logistic_rnd.m + lognormal_cdf.m + lognormal_inv.m + lognormal_pdf.m + lognormal_rnd.m + normal_cdf.m + normal_inv.m + normal_pdf.m + normal_rnd.m + pascal_cdf.m + pascal_inv.m + pascal_pdf.m + pascal_rnd.m + poisson_cdf.m + poisson_inv.m + poisson_pdf.m + poisson_rnd.m + stdnormal_cdf.m + stdnormal_inv.m + stdnormal_pdf.m + stdnormal_rnd.m + t_cdf.m + t_inv.m + t_pdf.m + t_rnd.m + uniform_cdf.m + uniform_inv.m + uniform_pdf.m + uniform_rnd.m + weibull_cdf.m + weibull_inv.m + weibull_pdf.m + weibull_rnd.m + wiener_rnd.m + + * statistics/models (new directory): + + logistic_regression.m + logistic_regression_derivatives.m + logistic_regression_likelihood.m + + * statistics/tests (new directory): + + anova.m + bartlett_test.m + chisquare_test_homogeneity.m + chisquare_test_independence.m + cor_test.m + f_test_regression.m + hotelling_test.m + hotelling_test_2.m + kolmogorov_smirnov_test.m + kolmogorov_smirnov_test_2.m + kruskal_wallis_test.m + manova.m + mcnemar_test.m + prop_test_2.m + run_test.m + sign_test.m + t_test.m + t_test_2.m + t_test_regression.m + u_test.m + var_test.m + welch_test.m + wilcoxon_test.m + z_test.m + z_test_2.m + +Thu Oct 22 12:25:55 1998 John W. Eaton <jwe@bevo.che.wisc.edu> + + * time/date.m: Use %Y, not %y in format string, for Matlab 5 + compatibility and to avoid Y2K problems. + Mon Oct 19 17:26:35 1998 John W. Eaton <jwe@bevo.che.wisc.edu> * polynomial/polyfit.m: Just use the \ operator to handle the
--- a/scripts/configure.in Mon Oct 19 22:37:57 1998 +0000 +++ b/scripts/configure.in Fri Oct 23 05:44:01 1998 +0000 @@ -27,8 +27,10 @@ AC_PROG_INSTALL AC_OUTPUT(Makefile audio/Makefile control/Makefile elfun/Makefile - general/Makefile image/Makefile io/Makefile + finance/Makefile general/Makefile image/Makefile io/Makefile linear-algebra/Makefile miscellaneous/Makefile plot/Makefile polynomial/Makefile set/Makefile signal/Makefile specfun/Makefile special-matrix/Makefile startup/Makefile - statistics/Makefile strings/Makefile time/Makefile) + statistics/Makefile statistics/base/Makefile + statistics/distributions/Makefile statistics/models/Makefile + statistics/tests/Makefile strings/Makefile time/Makefile)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/finance/fv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,76 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: fv (r, n, p [, l] [, method]) +## +## Returns the future value at the end of period n of an investment +## which consisting of n payments of p in each period, assuming an +## interest rate r. +## +## With the optional scalar argument l, one can specify an additional +## lump-sum payment. With the optional argument `method', one can +## specify whether the payments are made at the end ("e", default) or at +## the beginning ("b") of each period. +## +## Note that the rate r is not specified in percent, i.e., one has to +## write 0.05 rather than 5 %. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Future value of an investment + +function v = fv (r, n, p, l, m) + + if ((nargin < 3) || (nargin > 5)) + usage ("fv (r, n, p [, l] [, method])"); + endif + + if !(is_scalar (r) && (r > -1)) + error ("fv: r must be a scalar > -1"); + elseif !(is_scalar (n) && (n > 0)) + error ("fv: n must be a positive scalar"); + elseif !is_scalar (p) + error ("fv: p must be a scalar."); + endif + + if (r != 0) + v = p * ((1 + r)^n - 1) / r; + else + v = p * n; + endif + + if (nargin > 3) + if (nargin == 5) + if !isstr (m) + error ("fv: `method' must be a string"); + endif + elseif isstr (l) + m = l; + l = 0; + else + m = "e"; + endif + if strcmp (m, "b") + v = v * (1 + r); + endif + if is_scalar (l) + v = v + fvl (r, n, l); + else + error ("fv: l must be a scalar"); + endif + endif + +endfunction +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/finance/fvl.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,41 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: fvl (r, n, l) +## +## Returns the future value at the end of n periods of an initial lump +## sum investment l, given a per-period interest rate r. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Future value of an initial lump sum investment + +function v = fvl (r, n, l) + + if (nargin != 3) + usage ("fvl (r, n, l)"); + endif + + if !(is_scalar (r) && (r > -1)) + error ("fvl: r has to be a scalar > -1"); + elseif !(is_scalar (n) && (n > 0)) + error ("fvl: n has to be a positive scalar"); + elseif !is_scalar (l) + error ("fvl: l has to be a scalar"); + endif + + v = l * (1 + r)^n; + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/finance/irr.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,56 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: irr (p [, i]) +## +## Computes the internal rate of return of a series of payments p from +## an initial investment i, i.e., the solution of npv (r, p) = i. If the +## second argument is omitted, i = 0 is used. +## +## See also: npv; pv, rate. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Internal rate of return of an investment + +function r = irr (p, i) + + if (nargin == 1) + i = 0; + elseif !(nargin == 2) + usage ("irr (p [, i])"); + endif + + tmp = output_precision; + output_precision = 15; + if !(is_vector (p)) + error ("irr: p must be a vector"); + else + p_string = type p; + endif + + if !is_scalar (i) + error ("irr: i must be a scalar"); + endif + + string = ["function delta = f (r) ", ... + "delta = npv (r, %s) - %g; end"]; + eval (sprintf (string, p_string, i)); + + r = fsolve ("f", 0.01); + + output_precision = tmp; + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/finance/nper.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,79 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: nper (r, p, a [, l] [, method]) +## +## Computes the number of regular payments of p necessary to amortize a +## loan of amount a and interest r. +## +## With the optional scalar argument l, one can specify an additional +## lump-sum payment of l made at the end of the amortization time. With +## the optional string argument `method', one can specify whether +## payments are made at the end ("e", default) or at the beginning ("b") +## of each period. +## +## Note that the rate r is not specified in percent, i.e., one has to +## write 0.05 rather than 5 %. +## +## See also: pv, pmt, rate; npv. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Number of payments needed for amortizing a loan + +function n = nper (r, p, a, l, m) + + if ((nargin < 3) || (nargin > 5)) + usage ("nper (r, p, a [, l] [, method])"); + endif + + if !(is_scalar (r) && (r > -1)) + error ("nper: r must be a scalar > -1"); + elseif !is_scalar (p) + error ("nper: p must be a scalar"); + elseif !is_scalar (a) + error ("nper: a must be a scalar"); + endif + + if (nargin == 5) + if !isstr (m) + error ("nper: `method' must be a string"); + endif + elseif (nargin == 4) + if isstr (l) + m = l; + l = 0; + else + m = "e"; + endif + else + m = "e"; + l = 0; + endif + + if strcmp (m, "b") + p = p * (1 + r); + endif + + q = (p - r * a) / (p - r * l); + + if (q > 0) + n = - log (q) / log (1 + r); + else + n = Inf; + endif + +endfunction +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/finance/npv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,72 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: npv (r, p [, i]) +## +## Returns the net present value of a series of irregular (i.e., not +## necessarily identical) payments p which occur at the ends of n +## consecutive periods. r specifies the one-period interest rates and +## can either be a scalar (constant rates) or a vector of the same +## length as p. +## +## With the optional scalar argument i, one can specify an initial +## investment. +## +## Note that rates are not specified in percent, i.e., one has to write +## 0.05 rather than 5 %. +## +## See also: irr; pv. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Net present value of a series of payments + +function v = npv (r, p, i) + + if ((nargin < 2) || (nargin > 3)) + usage ("npv (r, p [, i]"); + endif + + if !(is_vector (p)) + error ("npv: p has to be a vector"); + else + n = length (p); + p = reshape (p, 1, n); + endif + + if any (any (r <= -1)) + error ("npv: all interest rates must be > -1"); + endif + if is_scalar (r) + d = 1 ./ (1 + r) .^ (0 : n); + elseif (is_vector (r) && (length (r) == n)) + d = [1, 1 ./ cumprod (reshape (1 + r, 1, n))]; + else + error (["npv: r has be a scalar ", ... + "or a vector of the same length as p"]); + endif + + if (nargin == 3) + if !is_scalar (i) + error ("npv: I_0 must be a scalar"); + endif + else + i = 0; + endif + + p = [i, p]; + v = sum (d .* p); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/finance/pmt.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,73 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: pmt (r, n, a [, l] [, method]) +## +## Compute the amount of periodic payment necessary to amortize a loan +## of amount a with interest rate r in n periods. +## +## With the optional scalar argument l, one can specify an initial +## lump-sum payment. With the optional string argument `method', one can +## specify whether payments are made at the end ("e", default) or at the +## beginning ("b") of each period. +## +## See also: pv, nper, rate + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Amount of periodic payment needed to amortize a loan + +function p = pmt (r, n, a, l, m) + + if ((nargin < 3) || (nargin > 5)) + usage ("pmt (r, n, a [, l] [, method])"); + endif + + if !(is_scalar (r) && (r > -1)) + error ("pmt: rate must be a scalar > -1"); + elseif !(is_scalar (n) && (n > 0)) + error ("pmt: n must be a positive scalar"); + elseif !(is_scalar (a) && (a > 0)) + error ("pmt: a must be a positive scalar."); + endif + + if (nargin == 5) + if !isstr (m) + error ("pmt: `method' must be a string"); + endif + elseif (nargin == 4) + if isstr (l) + m = l; + l = 0; + else + m = "e"; + endif + else + l = 0; + m = "e"; + endif + + p = r * (a - l * (1 + r)^(-n)) / (1 - (1 + r)^(-n)); + + if strcmp (m, "b") + p = p / (1 + r); + endif + + +endfunction + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/finance/pv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,79 @@ +## Copyright (C) 1995, 1996 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: pv (r, n, p [, l] [, method]) +## +## Returns the present value of an investment that will pay off p for n +## consecutive periods, assuming an interest r. +## +## With the optional scalar argument l, one can specify an additional +## lump-sum payment made at the end of n periods. +## +## With the optional string argument `method', one can specify whether +## payments are made at the end ("e", default) or at the beginning ("b") +## of each period. +## +## Note that the rate r is not specified in percent, i.e., one has to +## write 0.05 rather than 5 %. +## +## See also: pmt, nper, rate; npv. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Present value of an investment + +function v = pv (r, n, p, l, m) + + if ((nargin < 3) || (nargin > 5)) + usage ("pv (r, n, p [, l] [, method])"); + endif + + if !(is_scalar (r) && (r > -1)) + error ("pv: r must be a scalar > -1"); + elseif !(is_scalar (n) && (n > 0)) + error ("pv: n must be a positive scalar"); + elseif !is_scalar (p) + error ("pv: p must be a scalar."); + endif + + if (r != 0) + v = p * (1 - (1 + r)^(-n)) / r; + else + v = p * n; + endif + + if (nargin > 3) + if (nargin == 5) + if !isstr (m) + error ("pv: `method' must be a string"); + endif + elseif isstr (l) + m = l; + l = 0; + else + m = "e"; + endif + if strcmp (m, "b") + v = v * (1 + r); + endif + if is_scalar (l) + v = v + pvl (r, n, l); + else + error ("pv: l must be a scalar"); + endif + endif + +endfunction +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/finance/pvl.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,41 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: pvl (r, n, p) +## +## Returns the present value of an investment that will pay off p in one +## lump sum at the end of n periods, given the interest rate r. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Present value of an investment that pays off at the end + +function v = pvl (r, n, p) + + if (nargin != 3) + usage ("pvl (r, n, p)"); + endif + + if !(is_scalar (r) && (r > -1)) + error ("pvl: r has to be a scalar > -1"); + elseif !(is_scalar (n) && (n > 0)) + error ("pvl: n has to be a positive scalar"); + elseif !is_scalar (p) + error ("pvl: p has to be a scalar"); + endif + + v = p / (1 + r)^n; + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/finance/rate.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,72 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: rate (n, p, v [, l] [, method]) +## +## Computes the rate of return on an investment of present value v which +## pays p in n consecutive periods. +## +## With the optional scalar argument l, one can specify an additional +## lump-sum payment made at the end of n periods. With the optional +## string argument `method', one can specify whether payments are made +## at the end ("e", default) or at the beginning ("b") of each period. +## +## See also: pv, pmt, nper; npv. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Rate of return of an investment + +function r = rate (n, p, v, l, m) + + if ((nargin < 3) || (nargin > 5)) + usage ("rate (n, p, v [, l] [, method])"); + endif + + if !(is_scalar (n) && (n > 0)) + error ("rate: n must be a positive scalar"); + elseif !is_scalar (p) + error ("rate: p must be a scalar"); + elseif !is_scalar (v) + error ("rate: p must be a scalar"); + endif + + if (nargin == 5) + if !isstr (m) + error ("rate: `method' must be a string"); + endif + elseif (nargin == 4) + if isstr (l) + m = l; + l = 0; + else + m = "e"; + endif + else + l = 0; + m = "e"; + endif + + if !is_scalar (l) + error ("rate: l must be a scalar"); + endif + + string = ["function delta = f (r) ", ... + "delta = pv (r, %g, %g, %g, \"%s\") - %g; end"]; + eval (sprintf (string, n, p, l, m, v)); + + [r, info] = fsolve ("f", 0); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/finance/vol.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,70 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: vol (X, m [, n]) +## +## vol returns the volatility of each column of the input matrix X. m is +## the number of data sets per period (e.g. the number of data per year +## if you want to compute the volatility per year). The optional +## parameter n gives the number of past periods used for computation, if +## n is omitted, n=1 is used. If T is the number of rows of X, vol +## returns the volatility from n*m to T. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Volatility of financial time series data + +function retval = vol (X, m, n) + + if (nargin < 2) + usage ("vol (X, m [, n])"); + endif + + [xr, xc] = size (X); + + if (nargin > 2) + if (n * m > xr) + error ("vol: I need more data!"); + endif + else + n = 1; + if (n * m > xr) + error ("vol: I need more data!"); + endif + endif + + U = zeros (xr - 1, xc); + + if all (X) + U = X ((2 : xr), :) ./ X((1 : (xr-1)), :); + else + error ("vol: zero element in X"); + endif + + U = log(U); + U = U - ones (xr - 1, 1) * sum (U) / (xr - 1); + + retval = zeros (xr - n * m, xc); + + retval(1, :) = sumsq (U((1 : n*m), :)); + for i = 2 : (xr - n * m) + retval(i, :) = retval(i - 1, :) ... + - U(i - 1, :).^2 + U(i + n * m - 1, :).^2; + endfor + + retval = sqrt (retval * m / (n * m - 1)); + +endfunction +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/linear-algebra/dmult.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,38 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: dmult (a, B) +## +## If a is a vector of length rows (B), return diag (a) * B (but +## computed much more efficiently). + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Rescale the rows of a matrix + +function M = dmult (a, B) + + if (nargin != 2) + usage ("dmult (a, B)"); + endif + + s = size (a); + if ((min (s) > 1) || (max (s) != rows (B))) + error ("dmult: a must be a vector of length rows (B)"); + endif + + M = (reshape (a, max (s), 1) * ones (1, columns (B))) .* B; + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/arch_fit.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,108 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: [a, b] = arch_fit (y, X, p [, ITER [, gamma [, a0, b0]]]) +## +## Fits an ARCH regression model to the time series y using the scoring +## algorithm in Engle's original ARCH paper. The model is +## y(t) = b(1) * x(t,1) + ... + b(k) * x(t,k) + e(t), +## h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(p+1) * e(t-p)^2, +## where e(t) is N(0, h(t)), given y up to time t-1 and X up to t. +## +## If invoked as arch_fit (y, k, p) with a positive integer k, fit an +## ARCH(k,p) process, i.e., do the above with the t-th row of X given by +## [1, y(t-1), ..., y(t-k)]. +## +## Optionally, one can specify the number of iterations ITER, the +## updating factor gamma, and initial values a0 and b0 for the scoring +## algorithm. +## +## The input parameters are: +## y ... time series (vector) +## X ... matrix of (ordinary) regressors or order of +## autoregression +## p ... order of the regression of the residual variance + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Fit an ARCH regression model + +function [a, b] = arch_fit (y, X, p, ITER, gamma, a0, b0) + + if ((nargin < 3) || (nargin == 6) || (nargin > 7)) + usage ("arch_fit (y, X, p [, ITER [, gamma [, a0, b0]]])"); + endif + + if !(is_vector (y)) + error ("arch_test: y must be a vector"); + endif + + T = length (y); + y = reshape (y, T, 1); + [rx, cx] = size (X); + if ((rx == 1) && (cx == 1)) + X = autoreg_matrix (y, X); + elseif !(rx == T) + error (["arch_test: ", ... + "either rows (X) == length (y), or X is a scalar"]); + endif + + [T, k] = size (X); + + if (nargin == 7) + a = a0; + b = b0; + e = y - X * b; + else + [b, v_b, e] = ols (y, X); + a = [v_b, zeros (1,p)]'; + if (nargin < 5) + gamma = 0.1; + if (nargin < 4) + ITER = 50; + endif + endif + endif + + esq = e.^2; + Z = autoreg_matrix (esq, p); + + for i = 1 : ITER; + h = Z * a; + tmp = esq ./ h.^2 - 1 ./ h; + s = 1 ./ h(1:T-p); + for j = 1 : p; + s = s - a(j+1) * tmp(j+1:T-p+j); + endfor + r = 1 ./ h(1:T-p); + for j=1:p; + r = r + 2 * h(j+1:T-p+j).^2 .* esq(1:T-p); + endfor + r = sqrt (r); + X_tilde = X(1:T-p, :) .* (r * ones (1,k)); + e_tilde = e(1:T-p) .*s ./ r; + delta_b = inv (X_tilde' * X_tilde) * X_tilde' * e_tilde; + b = b + gamma * delta_b; + e = y - X * b; + esq = e .^ 2; + Z = autoreg_matrix (esq, p); + h = Z * a; + f = esq ./ h - ones(T,1); + Z_tilde = Z ./ (h * ones (1, p+1)); + delta_a = inv (Z_tilde' * Z_tilde) * Z_tilde' * f; + a = a + gamma * delta_a; + endfor + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/arch_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,87 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: y = arch_rnd (a, b, T) +## +## Simulates an ARCH sequence y of length T with AR coefficients b and +## CH coefficients a. +## I.e., y follows the model +## y(t) = b(1) + b(2) * y(t-1) + ... + b(lb) * y(t-lb+1) + e(t), +## where e(t), given y up to time t-1, is N(0, h(t)), with +## h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(la) * e(t-la+1)^2. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Simulate an ARCH process + +function y = arch_rnd (a, b, T) + + if (nargin != 3) + usage ("arch_rnd (a, b, T)"); + endif + + if !( (min (size (a)) == 1) && (min (size (b)) == 1) ) + error ("arch_rnd: a and b must both be scalars or vectors"); + endif + if !( is_scalar (T) && (T > 0) && (rem (T, 1) == 0) ) + error ("arch_rnd: T must be a positive integer"); + endif + + if !(a(1) > 0) + error ("arch_rnd: a(1) must be positive"); + endif + ## perhaps add a test for the roots of a(z) here ... + + la = length (a); + a = reshape (a, 1, la); + if (la == 1) + a = [a, 0]; + la = la + 1; + endif + lb = length (b); + b = reshape (b, 1, lb); + if (lb == 1) + b = [b, 0]; + lb = lb + 1; + endif + M = max([la lb]); + + e = zeros (T, 1); + h = zeros (T, 1); + y = zeros (T, 1); + + h(1) = a(1); + e(1) = sqrt (h(1)) * randn; + y(1) = b(1) + e(1); + + for t= 2 : M; + ta = min ([t la]); + h(t) = a(1) + a(2:ta) * e(t-1:t-ta+1).^2; + e(t) = sqrt (h(t)) * randn; + tb = min ([t lb]); + y(t) = b(1) + b(2:tb) * y(t-1:t-tb+1) + e(t); + endfor + if (T > M) + for t = M+1 : T; + h(t) = a(1) + a(2:la) * e(t-1:t-la+1).^2; + e(t) = sqrt (h(t)) * randn; + y(t) = b(1) + b(2:lb) * y(t-1:t-tb+1) + e(t); + endfor + endif + + y = y(1:T); + +endfunction + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/arch_test.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,70 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: [pval, lm] = arch_test (y, X, p) +## [pval, lm] = arch_test (y, k, p) +## +## arch_test (y, X, p) performs a Lagrange Multiplier (LM) test of the +## null hypothesis of no conditional heteroscedascity in the linear +## regression model y = X * b + e against the alternative of CH(p). +## I.e., the model is +## y(t) = b(1) * x(t,1) + ... + b(k) * x(t,k) + e(t), +## where given y up to t-1 and x up to t, e(t) is N(0, h(t)) with +## h(t) = v + a(1) * e(t-1)^2 + ... + a(p) * e(t-p)^2, +## and the null is a(1) == ... == a(p) == 0. +## +## arch_test (y, k, p) does the same in a linear autoregression model of +## order k, i.e., with [1, y(t-1), ..., y(t-k)] as the t-th row of X. +## +## Under the null, lm approximately has a chisquare distribution with p +## degrees of freedom. pval is the p-value (1 minus the CDF of this +## distribution at lm) of the test. +## +## If no output argument is given, the p-value is displayed. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Test for conditional heteroscedascity + +function [pval, lm] = arch_test (y, X, p) + + if (nargin != 3) + error ("arch_test needs 3 input arguments"); + endif + + if !(is_vector (y)) + error ("arch_test: y must be a vector"); + endif + T = length (y); + y = reshape (y, T, 1); + [rx, cx] = size (X); + if ((rx == 1) && (cx == 1)) + X = autoreg_matrix (y, X); + elseif !(rx == T) + error (["arch_test: ", ... + "either rows(X) == length(y), or X is a scalar"]); + endif + if !(is_scalar(p) && (rem(p, 1) == 0) && (p > 0)) + error ("arch_test: p must be a positive integer."); + endif + + [b, v_b, e] = ols (y, X); + Z = autoreg_matrix (e.^2, p); + f = e.^2 / v_b - ones (T, 1); + f = Z' * f; + lm = f' * inv (Z'*Z) * f / 2; + pval = 1 - chisquare_cdf (lm, p); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/arma_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,79 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: arma_rnd (a, b, v, t [, n]) +## +## Returns a simulation of the ARMA model +## x(n) = a(1) * x(n-1) + ... + a(k) * x(n-k) + +## + e(n) + b(1) * e(n-1) + ... + b(l) * e(n-l), +## where k is the length of vector a, l is the length of vector b and e +## is gaussian white noise with variance v. The function returns a +## vector of length t. +## +## The optional parameter n gives the number of dummy x(i) used for +## initialization, i.e., a sequence of length t+n is generated and +## x(n+1:t+n) is returned. If n is omitted, n=100 is used. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Simulate an ARMA process + +function x = arma_rnd (a, b, v, t, n) + + unwind_protect + orig_listelemok = empty_list_elements_ok; + empty_list_elements_ok = "true"; + + if (nargin == 4) + n = 100; + elseif (nargin == 5) + if (!is_scalar (t)) + error ("arma_rnd: n must be a scalar"); + endif + else + usage ("arma_rnd (a, b, v, t [, n])"); + endif + + if ( (min (size (a)) > 1) || (min (size (b)) > 1) ) + error ("arma_rnd: a and b must not be matrices"); + endif + + if (!is_scalar (t)) + error ("arma_rnd: t must be a scalar"); + endif + + ar = length (a); + br = length (b); + + a = reshape (a, ar, 1); + b = reshape (b, br, 1); + + a = [1; -a]; # apply our notational convention + b = [1; b]; + + n = min (n, ar + br); + + e = sqrt (v) * randn (t + n, 1); + + x = filter (b, a, e); + x = x(n + 1 : t + n); + + unwind_protect_cleanup + + empty_list_elements_ok = orig_listelemok; + + end_unwind_protect + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/autocor.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,44 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: autocor (X [, h]) +## +## returns the autocorrelations from lag 0 to h of vector X. +## If h is omitted, all autocorrelations are computed. +## If X is a matrix, the autocorrelations of every single column are +## computed. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Compute autocorrelations + +function retval = autocor (X, h) + + if (nargin == 1) + retval = autocov (X); + elseif (nargin == 2) + retval = autocov (X, h); + else + usage ("autocor (X [, h])"); + endif + + if (min (retval (1,:)) != 0) + retval = retval ./ ( ones (rows (retval), 1) * retval(1, :) ); + endif + +endfunction + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/autocov.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,60 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: autocov (X [, h]) +## +## returns the autocovariances from lag 0 to h of vector X. +## If h is omitted, all autocovariances are computed. +## If X is a matrix, the autocovariances of every single column are +## computed. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Compute autocovariances + +function retval = autocov (X, h) + + [n c] = size (X); + + if (is_vector (X)) + n = length (X); + c = 1; + X = reshape (X, n, 1); + endif + + X = center (X); + + if (nargin == 1) + h = n - 1; + endif + + retval = zeros (h + 1, c); + + unwind_protect + + oldpcv = prefer_column_vectors; + prefer_column_vectors = "false"; + + for i = 0 : h + retval(i+1, :) = diag (X(i+1:n, :).' * conj (X(1:n-i, :))) / n; + endfor + + unwind_protect_cleanup + + prefer_column_vectors = oldpcv; + + end_unwind_protect + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/autoreg_matrix.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,44 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: X = autoreg_matrix (y, k) +## +## Given a time series (vector) y, returns a matrix X with ones in the +## first column and the first k lagged values of y in the other columns. +## I.e., for t > k, [1, y(t-1), ..., y(t-k)] is the t-th row of X. X can +## be used as regressor matrix in autoregressions. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Design matrix for autoregressions + +function X = autoreg_matrix (y, k) + + if (nargin != 2) + usage ("autoreg_matrix (y, k)"); + endif + + if !(is_vector (y)) + error ("autoreg_matrix: y must be a vector"); + endif + + T = length (y); + y = reshape (y, T, 1); + X = ones (T, k+1); + for j = 1 : k; + X(:, j+1) = [zeros (j, 1); y(1:T-j)]; + endfor + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/bartlett.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,49 @@ +## Copyright (C) 1995, 1996, 1997 Andreas Weingessel +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: bartlett (m) +## +## Returns the filter coefficients of a Bartlett (triangular) window of +## length m. +## +## For a definition of the Bartlett window, see e.g. A. V. Oppenheim & +## R. W. Schafer, "Discrete-Time Signal Processing". + +## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at> +## Description: Coefficients of the Bartlett (triangular) window + +function c = bartlett (m) + + if (nargin != 1) + usage ("bartlett (m)"); + endif + + if !(is_scalar (m) && (m == round (m)) && (m > 0)) + error ("bartlett: m has to be an integer > 0"); + endif + + if (m == 1) + c = 1; + else + m = m - 1; + n = fix (m / 2); + c (1 : n+1) = 2 * (0 : n)' / m; + c (n+2 : m+1) = 2 - 2 * (n+1 : m)'/m; + endif + + c = c'; + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/blackman.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,45 @@ +## Copyright (C) 1995, 1996, 1997 Andreas Weingessel +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: blackman (m) +## +## Returns the filter coefficients of a Blackman window of length m. +## +## For a definition of the Blackman window, see e.g. A. V. Oppenheim & +## R. W. Schafer, "Discrete-Time Signal Processing". + +## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at> +## Description: Coefficients of the Blackman window + +function c = blackman (m) + + if (nargin != 1) + usage ("blackman (m)"); + endif + + if !(is_scalar (m) && (m == round (m)) && (m > 0)) + error ("blackman: m has to be an integer > 0"); + endif + + if (m == 1) + c = 1; + else + m = m - 1; + k = (0 : m)' / m; + c = 0.42 - 0.5 * cos (2 * pi * k) + 0.08 * cos (4 * pi * k); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/diffpara.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,83 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: [d, D] = diffpara (X [, a [, b]]) +## +## Returns the estimator d for the differencing parameter of an +## integrated time series. +## +## The frequencies from [2*pi*a/T, 2*pi*b/T] are used for the +## estimation. If b is omitted, the interval [2*pi/T, 2*pi*a/T] is used, +## if both b and a are omitted then a = 0.5 * sqrt(T) and b = 1.5 * +## sqrt(T) is used, where T is the sample size. If X is a matrix, the +## differencing parameter of every single column is estimated. +## +## D contains the estimators for all frequencies in the intervals +## described above, d is simply mean(D). +## +## Reference: Brockwell, Peter J. & Davis, Richard A. Time Series: +## Theory and Methods Springer 1987 + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Estimate the fractional differencing parameter + +function [d, D] = diffpara (X, a, b) + + if ((nargin < 1) || (nargin > 3)) + usage ("[d [, D]] = diffpara (X [, a [, b]])"); + else + if is_vector (X) + n = length (X); + k = 1; + X = reshape (X, n, 1); + else + [n k] = size(X); + endif + if (nargin == 1) + a = 0.5 * sqrt (n); + b = 1.5 * sqrt (n); + elseif (nargin == 2) + b = a; + a = 1; + endif + endif + + if !(is_scalar (a) && is_scalar (b)) + error ("diffpara: a and b must be scalars"); + endif + + D = zeros (b - a + 1, k); + + for l = 1:k + + w = 2 * pi * (1 : n-1) / n; + + x = 2 * log (abs( 1 - exp (-i*w))); + y = log (periodogram (X(2:n,l))); + + x = center (x); + y = center (y); + + for m = a:b + D(m-a+1) = - x(1:m) * y(1:m) / sumsq (x(1:m)); + endfor + + endfor + + d = mean (D); + +endfunction +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/durbinlevinson.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,89 @@ +## Copyright (C) 1995 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: durbinlevinson (c, [oldphi, oldv]) +## +## Performs one step of the Durbin-Levinson algorithm. +## +## The vector c_t = [gamma_0, ..., gamma_t] contains the autocovariances +## from lag 0 to t, oldphi the coefficients based on c_(t-1) and oldv +## the corresponding error. +## +## If oldphi is omitted, all steps from 1 to t of the algorithm are +## performed. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Perform one step of the Durbin-Levinson algorithm + +function [newphi, newv] = durbinlevinson (c, oldphi, oldv) + + if( !((nargin == 1) || (nargin == 3)) ) + usage ("durbinlevinson (c, [oldphi, oldv])"); + endif + + if( columns (c) > 1 ) + c=c'; + endif + + newphi = 0; + newv = 0; + + if (nargin == 3) + + t = length (oldphi) + 1; + + if (length (c) < t+1) + error ("durbilevinson: c too small"); + endif + + if (oldv == 0) + error ("durbinlevinson: oldv = 0"); + endif + + if (rows (oldphi) > 1 ) + oldphi = oldphi'; + endif + + newphi = zeros (1, t); + newphi(1) = ( c(t+1) - oldphi * c(2:t) ) / oldv; + for i = 2 : t + newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1); + endfor + newv = ( 1 - newphi(1)^2 ) * oldv; + + elseif(nargin == 1) + + tt = length (c)-1; + oldphi = c(2) / c(1); + oldv = ( 1 - oldphi^2 ) * c(1); + + for t = 2 : tt + + newphi = zeros (1, t); + newphi(1) = ( c(t+1) - oldphi * c(2:t) ) / oldv; + for i = 2 : t + newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1); + endfor + newv = ( 1 - newphi(1)^2 ) * oldv; + + oldv = newv; + oldphi = newphi; + + endfor + + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/fftshift.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,55 @@ +## Copyright (C) 1997 by Vincent Cautaerts +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: y = fftshift(W) +## +## Performs a shift of the vector V, for use with the "fft" and "ifft" +## functions, in order the move the frequency 0 to the centre of +## the vector or matrix. +## +## If V is a vector of E datas corresponding to E time samples spaced +## of Dt each, then fftshift (fft (V)) correspond to frequencies +## +## f = linspace (-E/(4*DT), (E/2-1)/(2*DT), E) +## +## If V is a matrix, does the same holds for rows and columns. + +## Author: Vincent Cautaerts <vincent@comf5.comm.eng.osaka-u.ac.jp> +## Created: July 1997 +## Adapted-By: jwe + +function retval = fftshift (V) + + retval = 0; + + if (nargin != 1) + usage ("usage: fftshift (X)"); + endif + + if (is_vector (V)) + x = length (V); + xx = ceil (x/2); + retval = V([xx+1:x, 1:xx]); + elseif (is_matrix (V)) + [x, y] = size (V); + xx = ceil (x/2); + yy = ceil (y/2); + retval = V([xx+1:x, 1:xx], [yy+1:y, 1:yy]); + else + error ("fftshift: expecting vector or matrix argument"); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/fractdiff.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,62 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: fractdiff(x, d) +## +## Computes the fractional differences (1-L)^d x where L denotes the +## lag-operator and d > -1. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Compute fractional differences + +function retval = fractdiff (x, d) + + N = 100; + + if !is_vector (x) + error ("fractdiff: x must be a vector") + endif + + if !is_scalar (d) + error ("fractdiff: d must be a scalar") + endif + + + if (d >= 1) + for k = 1 : d + x = x(2 : length (x)) - x(1 : length (x) - 1); + endfor + endif + + if (d > -1) + + d = rem (d, 1); + + if (d != 0) + n = (0 : N)'; + w = real (gamma (-d+n) ./ gamma (-d) ./ gamma (n+1)); + retval = fftfilt (w, x); + retval = retval(1 : length (x)); + else + retval = x; + endif + + else + error ("fractdiff: d must be > -1"); + + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/hamming.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,44 @@ +## Copyright (C) 1995, 1996, 1997 Andreas Weingessel +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: hamming (m) +## +## Returns the filter coefficients of a Hamming window of length m. +## +## For a definition of the Hamming window, see e.g. A. V. Oppenheim & +## R. W. Schafer, "Discrete-Time Signal Processing". + +## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at> +## Description: Coefficients of the Hamming window + +function c = hamming (m) + + if (nargin != 1) + usage ("hamming (m)"); + endif + + if !(is_scalar (m) && (m == round (m)) && (m > 0)) + error ("hamming: m has to be an integer > 0"); + endif + + if (m == 1) + c = 1; + else + m = m - 1; + c = 0.54 - 0.46 * cos (2 * pi * (0:m)' / m); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/hanning.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,44 @@ +## Copyright (C) 1995, 1996, 1997 Andreas Weingessel +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: hanning (m) +## +## Returns the filter coefficients of a Hanning window of length m. +## +## For a definition of this window type, see e.g. A. V. Oppenheim & +## R. W. Schafer, "Discrete-Time Signal Processing". + +## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at> +## Description: Coefficients of the Hanning window + +function c = hanning (m) + + if (nargin != 1) + usage ("hanning (m)"); + endif + + if !(is_scalar (m) && (m == round (m)) && (m > 0)) + error ("hanning: m has to be an integer > 0"); + endif + + if (m == 1) + c = 1; + else + m = m - 1; + c = 0.5 - 0.5 * cos (2 * pi * (0 : m)' / m); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/hurst.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,45 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: H = hurst (x) +## +## Estimates the Hurst parameter of sample x via the rescaled range +## statistic. If x is a matrix, the parameter is estimated for every +## single column. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Estimate the Hurst parameter + +function H = hurst (x) + + if (nargin != 1) + usage ("hurst (x)"); + endif + + if (is_scalar (x)) + error ("hurst: x must not be a scalar") + elseif is_vector (x) + x = reshape (x, length (x), 1); + end + + [xr xc] = size (x); + + s = std (x); + w = cumsum (x - mean (x)); + RS = (max(w) - min(w)) ./ s; + H = log (RS) / log (xr); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/periodogram.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,41 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: periodogram (x) +## +## For a data matrix x from a sample of size n, return the periodogram. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Compute the periodogram + +function retval = periodogram (x) + + [r c] = size(x); + + if (r == 1) + r = c; + endif + + retval = (abs (fft (x - mean (x)))) .^ 2 / r; + +endfunction + + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/rectangle_lw.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,32 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: retval = rectangle_lw (n, b) +## +## Rectangular lag window. Subfunction used for spectral density +## estimation. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Rectangular lag window + +function retval = rectangle_lw (n, b) + + retval = zeros (n, 1); + t = floor (1 / b); + + retval (1:t, 1) = ones (t, 1); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/rectangle_sw.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: retval = rectangle_sw (n, b) +## +## Rectangular spectral window. Subfunction used for spectral density +## estimation. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Rectangular spectral window + +function retval = rectangle_sw (n, b) + + retval = zeros (n, 1); + retval(1) = 2 / b + 1; + + l = (2:n)' - 1; + l = 2 * pi * l / n; + + retval(2:n) = sin( (2/b + 1) * l / 2 ) ./ sin (l / 2); + +endfunction + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/sinetone.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,60 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: sinetone (freq [, rate [, sec [, ampl]]]) +## +## Returns a sinetone of frequency freq with length of sec seconds at +## sampling rate rate and with amplitude ampl. freq and ampl may be +## vectors of common size. +## +## Defaults are rate = 8000, sec = 1 and ampl = 64. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Compute a sine tone + +function retval = sinetone (f, r, s, a) + + if (nargin == 1) + r = 8000; + s = 1; + a = 64; + elseif (nargin == 2) + s = 1; + a = 64; + elseif (nargin == 3) + a = 64; + elseif ((nargin < 1) || (nargin > 4)) + usage ("sinetone (freq [, rate [, sec [, ampl]]])"); + endif + + [err, f, a] = common_size (f, a); + if (err || ! is_vector (f)) + error ("sinetone: freq and ampl must be vectors of common size"); + endif + + if !(is_scalar (r) && is_scalar (s)) + error ("sinetone: rate and sec must be scalars"); + endif + + n = length (f); + + retval = zeros (r * s, n); + for k = 1:n + retval (:, k) = a(k) * sin (2 * pi * (1:r*s) / r * f(k))'; + endfor + +endfunction +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/sinewave.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,37 @@ +## Copyright (C) 1995, 1996, 1997 Andreas Weingessel +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: sinewave (m, n [, d]) +## +## Computes an (m x 1)-vector X with i-th element X(i) given by sin (2 * +## pi * (i+d-1) / n). +## +## The default value for d is 0. + +## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at> +## Description: Compute a sine wave + +function x = sinewave (m, n, d) + + if (nargin == 2) + d = 0; + elseif (nargin != 3) + usage ("sinewave (m, n [, d])"); + endif + + x = sin ( ((1 : m) + d - 1) * 2 * pi / n); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/spectral_adf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,61 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: retval = spectral_adf (c, [win, [b]]) +## +## Returns the spectral density estimator. +## c ....... vector of autocovariances (starting at lag 0) +## win ..... window name, eg. "triangle" or "rectangle" +## spectral_adf searches for a function called win_lw () +## b ....... bandwidth +## +## If win is omitted, the triangle window is used as default. +## If b is omitted, 1 / sqrt( length (c)) is used as default. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Spectral density estimation + +function retval = spectral_adf (c, win, b) + + cr = length (c); + + if (columns (c) > 1) + c=c'; + endif + + if (nargin < 3) + b = 1 / ceil (sqrt (cr)); + endif + + if (nargin == 1) + w = triangle_lw (cr, b); + else + win = [win, "_lw"]; + w = feval (win, cr, b); + endif + + c = c .* w; + + retval = 2 * real (fft (c)) - c(1); + retval = [zeros (cr, 1) retval]; + retval(:, 1) = (0 : xr-1)' / xr; + +endfunction + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/spectral_xdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: retval = spectral_xdf (X, [win, [b]]) +## +## Returns the spectral density estimator. +## X ....... data vector +## win ..... window name, eg. "triangle" or "rectangle" +## spectral_adf searches for a function called win_sw () +## b ....... bandwidth +## +## If win is omitted, the triangle window is used as default. +## If b is omitted, 1 / sqrt (length (X)) is used as default. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Spectral density estimation + +function retval = spectral_xdf (X, win, b) + + xr = length (X); + + if (columns (X) > 1) + X = X'; + endif + + if (nargin < 3) + b = 1 / ceil (sqrt (xr)); + endif + + if (nargin == 1) + w = triangle_sw (xr, b); + else + win = [win, "_sw"]; + w = feval (win, xr, b); + endif + + X = X - sum (X) / xr; + + retval = (abs (fft (X)) / xr).^2; + retval = real (ifft (fft(retval) .* fft(w))); + + retval = [zeros (xr, 1) retval]; + retval(:, 1) = (0 : xr-1)' / xr; + +endfunction + + + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/spencer.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,59 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: spencer (X) +## +## returns Spencer's 15 point moving average of every single column of X + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Apply Spencer's 15-point MA filter + +function retval = spencer (X) + + if (nargin != 1) + usage ("spencer (X)"); + endif + + [xr xc] = size(X); + + n = xr; + c = xc; + + if (is_vector(X)) + n = length(X); + c = 1; + X = reshape(X, n, 1); + endif + + W = [ -3 -6 -5 3 21 46 67 74 67 46 21 3 -5 -6 -3 ] / 320; + + retval = fftfilt (W, X); + retval = [zeros(7,c); retval(15:n,:); zeros(7,c);]; + + retval = reshape(retval, xr, xc); + +endfunction + + + + + + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/stft.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,114 @@ +## Copyright (C) 1995, 1996, 1997 Andreas Weingessel +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: [Y, c] = stft (X [, win_size [, inc [, num_coef [, w_type]]]]) +## +## Computes the short-term Fourier transform of the vector X with +## "num_coef" coefficients by applying a window of "win_size" data +## points and an increment of "inc" points. +## +## Before computing the Fourier transform, one of the following windows +## is applied: "hanning" (w_type = 1), "hamming" (w_type = 2), +## "rectangle" (w_type = 3). The window names can be passed as strings +## or by the w_type number. +## +## If not all arguments are specified, the following defaults are used: +## win_size = 80, inc = 24, num_coef = 64, w_type = 1. +## +## Y = stft (X [, ...]) returns the absolute values of the Fourier +## coefficients according to the num_coef positive frequencies. +## [Y, c] = stft (X [, ...]) returns the entire STFT-matrix Y and a +## vector c = [win_size, inc, w_type] which is needed by the synthesis +## function. + +## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at> +## Description: Short-term Fourier transform + +function [Y, c] = stft(X, win, inc, coef, w_type) + + ## default values of unspecified arguments + if (nargin < 5) + w_type = 1; + if (nargin < 4) + coef = 64; + if (nargin < 3) + inc = 24; + if (nargin < 2) + win = 80; + endif + endif + endif + elseif (nargin == 5) + if (isstr (w_type)) + if (strcmp (w_type, "hanning")) + w_type = 1; + elseif (strcmp (w_type, "hamming")) + w_type = 2; + elseif (strcmp (w_type, "rectangle")) + w_type = 3; + else + error (["stft: unknown window type `", w_type, "'"]) + endif + endif + else + usage ("[Y [, c]] = ", ... + "stft(X [, win_size [, inc [, num_coef [, w_type]]]])"); + endif + + ## check whether X is a vector + [nr, nc] = size (X); + if (nc != 1) + if (nr == 1) + X = X'; + nr = nc; + else + error ("stft: X must be a vector"); + endif + endif + + num_coef = 2 * coef; + if (win > num_coef) + win = num_coef; + printf ("stft: window size adjusted to %f\n", win); + endif + num_win = fix ((nr - win) / inc); + + ## compute the window coefficients + if (w_type == 3) # rectangular window + WIN_COEF = ones (win, 1); + elseif (w_type == 2) # Hamming window + WIN_COEF = hamming (win); + else # Hanning window + WIN_COEF = hanning (win); + endif + + ## create a matrix Z whose columns contain the windowed time-slices + Z = zeros (num_coef, num_win + 1); + start = 1; + for i = 0:num_win + Z(1:win, i+1) = X(start:start+win-1) .* WIN_COEF; + start = start + inc; + endfor + + Y = fft (Z); + + if (nargout == 1) + Y = abs (Y(1:coef, :)); + else + c = [win, inc, w_type]; + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/synthesis.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 1995, 1996, 1997 Andreas Weingessel +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: X = synthesis (Y, c) +## +## recovers a signal X from its short-time Fourier transform Y. c = +## [win_size, increment, window_type]. +## +## Y and c can be derived by [Y, c] = stft (X [, ...]). + +## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at> +## Description: Recover a signal from its short-term Fourier transform + +function X = synthesis (Y, c) + + if (nargin != 2) + usage ("X = synthesis (Y, c)"); + endif + + [nr nc] = size (c); + if (nr * nc != 3) + error ("synthesis: c must contain exactly 3 elements"); + endif + + ## not necessary, enables better reading + win = c(1); + inc = c(2); + w_type = c(3); + + if (w_type == 1) + H = hanning (win); + elseif (w_type == 2) + H = hamming (win); + elseif (w_type == 3) + H = ones (win, 1); + else + error ("synthesis: window_type must be 1, 2, or 3"); + endif + + Z = real (ifft (Y)); + st = fix ((win-inc) / 2); + Z = Z(st:st+inc-1, :); + H = H(st:st+inc-1); + + nc = columns(Z); + for i = 1:nc + Z(:, i) = Z(:, i) ./ H; + endfor + + X = reshape(Z, inc * nc, 1); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/triangle_lw.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,31 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: retval = triangle_lw (n, b) +## +## Triangular lag window. Subfunction used for spectral density +## estimation. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Triangular lag window + +function retval = triangle_lw (n, b) + + retval = 1 - (0 : n-1)' * b; + retval = max ([retval'; zeros (1, n)])'; + +endfunction +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/triangle_sw.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: retval = triangle_sw (n, b) +## +## Triangular spectral window. Subfunction used for spectral density +## estimation. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Triangular spectral window + +function retval = triangle_sw (n, b) + + retval = zeros(n,1); + retval(1) = 1 / b; + + l = (2:n)' - 1; + l = 2 * pi * l / n; + + retval(2:n) = b * (sin (l / (2*b)) ./ sin (l / 2) ).^2; + +endfunction + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/signal/yulewalker.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,52 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: [a, v] = yulewalker (c) +## +## fits an AR (p)-model with Yule-Walker estimates. +## c = [gamma_0, ..., gamma_p] autocovariances +## a .... AR coefficients +## v .... variance of white noise + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Fit AR model by Yule-Walker method + +function [a, v] = yulewalker (c) + + p = length (c) - 1; + + if (columns (c) > 1) + c = c'; + endif + + cp = c(2 : p+1); + CP = zeros(p, p); + + for i = 1:p + for j = 1:p + CP (i, j) = c (abs (i-j) + 1); + endfor + endfor + + a = inv (CP) * cp; + v = c(1) - a' * cp; + +endfunction + + + + +
--- a/scripts/statistics/Makefile.in Mon Oct 19 22:37:57 1998 +0000 +++ b/scripts/statistics/Makefile.in Fri Oct 23 05:44:01 1998 +0000 @@ -20,56 +20,56 @@ INSTALL_PROGRAM = @INSTALL_PROGRAM@ INSTALL_DATA = @INSTALL_DATA@ -SOURCES = *.m +SOURCES = # *.m DISTFILES = Makefile.in $(SOURCES) -FCN_FILES = $(wildcard $(srcdir)/*.m) -FCN_FILES_NO_DIR = $(notdir $(FCN_FILES)) +SUBDIRS = base distributions models tests + +DISTSUBDIRS = $(SUBDIRS) -BINDISTFILES = $(FCN_FILES) +BINDISTSUBDIRS = $(SUBDIRS) -all: +FCN_FILES = # $(wildcard $(srcdir)/*.m) +FCN_FILES_NO_DIR = # $(notdir $(FCN_FILES)) + +all: $(SUBDIRS) .PHONY: all -install install-strip: - $(top_srcdir)/mkinstalldirs $(fcnfiledir)/$(script_sub_dir) - for f in $(FCN_FILES_NO_DIR); do \ - rm -f $(fcnfiledir)/$(script_sub_dir)/$$f; \ - $(INSTALL_DATA) $(srcdir)/$$f $(fcnfiledir)/$(script_sub_dir)/$$f; \ - done -.PHONY: install install-strip +$(SUBDIRS): + echo making all in $@ + cd $@; $(MAKE) all +.PHONY: $(SUBDIRS) -uninstall: - for f in $(FCN_FILES_NO_DIR); \ - do rm -f $(fcnfiledir)/$(script_sub_dir)/$$f; \ - done -.PHONY: uninstall +install install-strip uninstall clean mostlyclean distclean maintainer-clean:: + @$(subdir-for-command) +.PHONY: install install-strip uninstall +.PHONY: clean mostlyclean distclean maintainer-clean + +install install-strip:: -clean: -.PHONY: clean +uninstall:: -tags: $(SOURCES) +tags TAGS:: $(SOURCES) + $(SUBDIR_FOR_COMMAND) + +tags:: ctags $(SOURCES) -TAGS: $(SOURCES) +TAGS:: $(SOURCES) etags $(SOURCES) -mostlyclean: clean -.PHONY: mostlyclean - -distclean: clean +distclean:: rm -f Makefile -.PHONY: distclean -maintainer-clean: distclean - rm -f tags TAGS -.PHONY: maintainer-clean +maintainer-clean:: + rm -f tags TAGS Makefile dist: - ln $(DISTFILES) ../../`cat ../../.fname`/scripts/statistics + ln $(DISTFILES) ../../`cat ../../.fname`/scripts + for dir in $(DISTSUBDIRS); do mkdir ../../`cat ../../.fname`/scripts/statistics/$$dir; cd $$dir; $(MAKE) $@; cd ..; done .PHONY: dist bin-dist: - ln $(BINDISTFILES) ../../`cat ../../.fname`/scripts/statistics + for dir in $(BINDISTSUBDIRS); do mkdir ../../`cat ../../.fname`/scripts/statistics/$$dir; cd $$dir; $(MAKE) $@; cd ..; done .PHONY: bin-dist
--- a/scripts/statistics/corrcoef.m Mon Oct 19 22:37:57 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,48 +0,0 @@ -## Copyright (C) 1996, 1997 John W. Eaton -## -## 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 2, 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, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: corrcoef (x, y) -## -## The (i,j)-th entry of corrcoef (x, y) is the correlation between the -## i-th variable in x and the j-th variable in y. -## For matrices, each row is an observation and each column a variable; -## vectors are always observations and may be row or column vectors. -## corrcoef (x) is corrcoef (x, x). - -## Author: Kurt Hornik <hornik@ci.tuwien.ac.at> -## Created: March 1993 -## Adapted-By: jwe - -function retval = corrcoef (x, y) - - if (nargin < 1 || nargin > 2) - usage ("corrcoef (x, y)"); - endif - - if (nargin == 2) - c = cov (x, y); - s = std (x)' * std (y); - retval = c ./ s; - elseif (nargin == 1) - c = cov (x); - s = reshape (sqrt (diag (c)), 1, columns (c)); - retval = c ./ (s' * s); - endif - -endfunction
--- a/scripts/statistics/cov.m Mon Oct 19 22:37:57 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,51 +0,0 @@ -## Copyright (C) 1996, 1997 John W. Eaton -## -## 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 2, 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, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: cov (X [, Y]) -## -## If each row of X and Y is an observation and each column is a -## variable, the (i,j)-th entry of cov(X, Y) is the covariance -## between the i-th variable in X and the j-th variable in Y. -## cov(X) is cov(X, X). - -## Author: Kurt Hornik <hornik@ci.tuwien.ac.at> -## Created: March 1993 -## Adapted-By: jwe - -function retval = cov (X, Y) - - if (nargin < 1 || nargin > 2) - usage ("cov (X [, Y])"); - endif - - [Tx, kx] = size (X); - if (nargin == 2) - [Ty, ky] = size (Y); - if (Tx != Ty) - error ("cov: X and Y must have the same number of rows."); - endif - X = X - ones (Tx, 1) * sum (X) / Tx; - Y = Y - ones (Tx, 1) * sum (Y) / Tx; - retval = conj (X' * Y / (Tx - 1)); - elseif (nargin == 1) - X = X - ones (Tx, 1) * sum (X) / Tx; - retval = conj (X' * X / (Tx - 1)); - endif - -endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/beta_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,60 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: beta_cdf (x, a, b) +## +## For each element of x, returns the CDF at x of the beta distribution +## with parameters a and b, i.e., PROB( beta(a,b) <= x). + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the Beta distribution + +function cdf = beta_cdf (x, a, b) + + if (nargin != 3) + usage ("beta_cdf (a, b, x)"); + endif + + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("beta_cdf: x, a and b must be of common size or scalar"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, s, 1); + a = reshape (a, s, 1); + b = reshape (b, s, 1); + cdf = zeros (s, 1); + + k = find (!(a > 0) | !(b > 0) | isnan (x)); + if any (k) + cdf (k) = NaN * ones (length (k), 1); + endif + + k = find ((x >= 1) & (a > 0) & (b > 0)); + if any (k) + cdf (k) = ones (length (k), 1); + endif + + k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0)); + if any (k) + cdf (k) = betai (a(k), b(k), x(k)); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/beta_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,92 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: beta_inv (x, a, b) +## +## For each component of x, compute the quantile (the inverse of the +## CDF) at x of the Beta distribution with parameters a and b. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Beta distribution + +function inv = beta_inv (x, a, b) + + if (nargin != 3) + usage ("beta_inv (x, a, b)"); + endif + + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("beta_inv: x, a and b must be of common size or scalars"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, s, 1); + a = reshape (a, s, 1); + b = reshape (b, s, 1); + inv = zeros (s, 1); + + k = find ((x < 0) | (x > 1) | !(a > 0) | !(b > 0) | isnan (x)); + if any (k) + inv (k) = NaN * ones (length (k), 1); + endif + + k = find ((x == 1) & (a > 0) & (b > 0)); + if any (k) + inv (k) = ones (length (k), 1); + endif + + k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0)); + if any (k) + a = a (k); + b = b (k); + x = x (k); + y = a ./ b; + l = find (y < eps); + if any (l) + y(l) = sqrt (eps) * ones (length (l), 1); + endif + l = find (y > 1 - eps); + if any (l) + y(l) = 1 - sqrt (eps) * ones (length (l), 1); + endif + + y_old = y; + for i = 1 : 100 + h = (beta_cdf (y_old, a, b) - x) ./ beta_pdf (y_old, a, b); + y_new = y_old - h; + ind = find (y_new <= eps); + if (any (ind)) + y_new (ind) = y_old (ind) / 10; + endif + ind = find (y_new >= 1 - eps); + if any (ind) + y_new (ind) = 1 - (1 - y_old (ind)) / 10; + endif + h = y_old - y_new; + if (max (abs (h)) < sqrt (eps)) + break; + endif + y_old = y_new; + endfor + + inv (k) = y_new; + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/beta_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,56 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: beta_pdf (x, a, b) +## +## For each element of x, returns the PDF at x of the beta distribution +## with parameters a and b. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the Beta distribution + +function pdf = beta_pdf (x, a, b) + + if (nargin != 3) + usage ("beta_pdf (a, b, x)"); + endif + + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("beta_pdf: x, a and b must be of common size or scalar"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, s, 1); + a = reshape (a, s, 1); + b = reshape (b, s, 1); + pdf = zeros (s, 1); + + k = find (!(a > 0) | !(b > 0) | isnan (x)); + if any (k) + pdf (k) = NaN * ones (length (k), 1); + endif + + k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0)); + if any (k) + pdf(k) = exp ((a(k) - 1) .* log (x(k)) ... + + (b(k) - 1) .* log (1 - x(k))) ./ beta (a(k), b(k)); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/beta_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,72 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: beta_rnd (a, b [, r, c]) +## +## beta_rnd (a, b) returns a matrix of random samples from the Beta +## distribution with parameters a and b. The size of the matrix is the +## common size of a and b. +## +## beta_rnd (a, b, r, c) returns an r by c matrix of random samples from +## the Beta distribution with parameters a and b. Both a and b must be +## scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the Beta distribution + +function rnd = beta_rnd (a, b, r, c) + + if (nargin == 4) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("beta_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("beta_rnd: c must be a positive integer"); + endif + [retval, a, b] = common_size (a, b, zeros (r, c)); + if (retval > 0) + error (strcat("beta_rnd: ", + "a and b must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin == 2) + [retval, a, b] = common_size (a, b); + if (retval > 0) + error ("beta_rnd: a and b must be of common size or scalar"); + endif + else + usage ("beta_rnd (a, b [, r, c])"); + endif + + [r, c] = size (a); + s = r * c; + a = reshape (a, 1, s); + b = reshape (b, 1, s); + rnd = zeros (1, s); + + k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)); + if (any (k)) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf)); + if (any (k)) + rnd(k) = beta_inv (rand (1, length (k)), a(k), b(k)); + endif + + rnd = reshape (rnd, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/binomial_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: binomial_cdf (x, n, p) +## +## For each element of x, compute the CDF at x of the binomial +## distribution with parameters n and p. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the binomial distribution + +function cdf = binomial_cdf (x, n, p) + + if (nargin != 3) + usage ("binomial_cdf (x, n, p)"); + endif + + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error (["binomial_cdf: ", ... + "x, n and p must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + n = reshape (n, 1, s); + p = reshape (p, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | !(n >= 0) | (n != round (n)) ... + | !(p >= 0) | !(p <= 1)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x >= n) & (n >= 0) & (n == round (n)) ... + & (p >= 0) & (p <= 1)); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + k = find ((x >= 0) & (x < n) & (n == round (n)) ... + & (p >= 0) & (p <= 1)); + if any (k) + tmp = floor (x(k)); + cdf(k) = 1 - betai (tmp + 1, n(k) - tmp, p(k)); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/binomial_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,67 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: binomial_inv (x, n, p) +## +## For each element of x, compute the quantile at x of the binomial +## distribution with parameters n and p. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the binomial distribution + +function inv = binomial_inv (x, n, p) + + if (nargin != 3) + usage ("binomial_inv (x, n, p)"); + endif + + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error (["binomial_inv: ", ... + "x, n and p must be of common size or scalars"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + n = reshape (n, 1, s); + p = reshape (p, 1, s); + inv = zeros (1, s); + + k = find (!(x >= 0) | !(x <= 1) | !(n >= 0) | (n != round (n)) ... + | !(p >= 0) | !(p <= 1)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find ((x >= 0) & (x <= 1) & (n >= 0) & (n == round (n)) ... + & (p >= 0) & (p <= 1)); + if any (k) + cdf = binomial_pdf (0, n(k), p(k)); + while (any (inv(k) < n(k))) + m = find (cdf < x(k)); + if any (m) + inv(k(m)) = inv(k(m)) + 1; + cdf(m) = cdf(m) + binomial_pdf (inv(k(m)), n(k(m)), p(k(m))); + else + break; + endif + endwhile + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/binomial_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,59 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: binomial_pdf (x, n, p) +## +## For each element of x, compute the probability density function (PDF) +## at x of the binomial distribution with parameters n and p. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the binomial distribution + +function pdf = binomial_pdf (x, n, p) + + if (nargin != 3) + usage ("binomial_pdf (x, n, p)"); + endif + + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error (["binomial_pdf: ", ... + "x, n and p must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + n = reshape (n, 1, s); + p = reshape (p, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | !(n >= 0) | (n != round (n)) ... + | !(p >= 0) | !(p <= 1)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x >= 0) & (x <= n) & (x == round (x)) ... + & (n == round (n)) & (p >= 0) & (p <= 1)); + if any (k) + pdf(k) = bincoeff (n(k), x(k)) .* (p(k) .^ x(k)) ... + .* ((1 - p(k)) .^ (n(k) - x(k))); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/binomial_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,78 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: binomial_rnd (n, p [, r, c]) +## +## binomial_rnd (n, p) returns a matrix of random samples from the +## binomial distribution with parameters n and p. The size of the +## matrix is the common size of n and p. +## +## binomial_rnd (n, p, r, c) returns an r by c matrix of random samples +## from the binomial distribution with parameters n and p. Both n and p +## must be scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the binomial distribution + +function rnd = binomial_rnd (n, p, r, c) + + if (nargin == 4) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("binomial_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("binomial_rnd: c must be a positive integer"); + endif + [retval, n, p] = common_size (n, p, zeros (r, c)); + if (retval > 0) + error (strcat("binomial_rnd: ", + "n and p must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin == 2) + [retval, n, p] = common_size (n, p); + if (retval > 0) + error ("binomial_rnd: n and p must be of common size or scalar"); + endif + else + usage ("binomial_rnd (n, p [, r, c])"); + endif + + [r, c] = size (n); + s = r * c; + n = reshape (n, 1, s); + p = reshape (p, 1, s); + rnd = zeros (1, s); + + k = find (!(n > 0) | !(n < Inf) | !(n == round (n)) | + !(p <= 0) | !(p >= 1)); + if any (k) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find ((n > 0) & (n < Inf) & (n == round (n)) & (p >= 0) & (p <= 1)); + if any (k) + N = max (n(k)); + L = length (k); + tmp = rand (N, L); + ind = (1 : N)' * ones (1, L); + rnd(k) = sum ((tmp < ones (N, 1) * p(k)) & + (ind <= ones (N, 1) * n(k))); + endif + + rnd = reshape (rnd, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/cauchy_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,58 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: cauchy_cdf (x [, lambda, sigma]) +## +## For each element of x, compute the cumulative distribution function +## (CDF) at x of the Cauchy distribution with location parameter lambda +## and scale parameter sigma. Default values are lambda = 0, sigma = 1. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the Cauchy distribution + +function cdf = cauchy_cdf (x, location, scale) + + if !(nargin == 1 || nargin == 3) + usage ("cauchy_cdf (x [, lambda, sigma])"); + endif + + if (nargin == 1) + location = 0; + scale = 1; + endif + + [retval, x, location, scale] = common_size (x, location, scale); + if (retval > 0) + error (["cauchy_cdf: ", ... + "x, lambda and sigma must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + location = reshape (location, 1, s); + scale = reshape (scale, 1, s); + cdf = NaN * ones (1, s); + + k = find ((x > -Inf) & (x < Inf) & (location > -Inf) & + (location < Inf) & (scale > 0) & (scale < Inf)); + if any (k) + cdf(k) = 0.5 + atan ((x(k) - location(k)) ./ scale(k)) / pi; + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/cauchy_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,71 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: cauchy_inv (x [, lambda, sigma]) +## +## For each element of x, compute the quantile (the inverse of the CDF) +## at x of the Cauchy distribution with location parameter lambda and +## scale parameter sigma. Default values are lambda = 0, sigma = 1. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Cauchy distribution + +function inv = cauchy_inv (x, location, scale) + + if !(nargin == 1 || nargin == 3) + usage ("cauchy_inv (x [, lambda, sigma])"); + endif + + if (nargin == 1) + location = 0; + scale = 1; + endif + + [retval, x, location, scale] = common_size (x, location, scale); + if (retval > 0) + error (["cauchy_inv: ", ... + "x, lambda and sigma must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + location = reshape (location, 1, s); + scale = reshape (scale, 1, s); + + inv = NaN * ones (1, s); + + ok = ((location > -Inf) & (location < Inf) & + (scale > 0) & (scale < Inf)); + + k = find ((x == 0) & ok); + if any (k) + inv(k) = -Inf * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < 1) & ok); + if any (k) + inv(k) = location(k) - scale(k) .* cot (pi * x(k)); + endif + + k = find ((x == 1) & ok); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/cauchy_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,60 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: cauchy_pdf (x [, lambda, sigma]) +## +## For each element of x, compute the probability density function (PDF) +## at x of the Cauchy distribution with location parameter lambda and +## scale parameter sigma > 0. Default values are lambda = 0, sigma = 1. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the Cauchy distribution + +function pdf = cauchy_pdf (x, location, scale) + + if !(nargin == 1 || nargin == 3) + usage ("cauchy_pdf (x [, lambda, sigma])"); + endif + + if (nargin == 1) + location = 0; + scale = 1; + endif + + [retval, x, location, scale] = common_size (x, location, scale); + if (retval > 0) + error (["cauchy_pdf: ", ... + "x, lambda and sigma must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + location = reshape (location, 1, s); + scale = reshape (scale, 1, s); + + pdf = NaN * ones (1, s); + + k = find ((x > -Inf) & (x < Inf) & (location > -Inf) & + (location < Inf) & (scale > 0) & (scale < Inf)); + if any (k) + pdf(k) = (1 ./ (1 + ((x(k) - location(k)) ./ scale(k)) .^ 2)) ... + / pi ./ scale(k); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/cauchy_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,69 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: cauchy_rnd (lambda, sigma [, r, c]) +## +## cauchy_rnd (lambda, sigma) returns a matrix of random samples from +## the Cauchy distribution with parameters lambda and sigma. The size +## of the matrix is the common size of the parameters. +## +## cauchy_rnd (lambda, sigma, r, c) returns an r by c matrix of random +## samples from the Cauchy distribution with parameters lambda and sigma +## which must both be scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the Cauchy distribution + +function rnd = cauchy_rnd (l, scale, r, c) + + if (nargin == 4) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("cauchy_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("cauchy_rnd: c must be a positive integer"); + endif + [retval, l, scale] = common_size (l, scale, zeros (r, c)); + if (retval > 0) + error (strcat("cauchy_rnd: ", + "lambda and sigma must be scalar or of size", + sprintf ("%d by %d.", r, c))); + endif + elseif (nargin == 2) + [retval, l, scale] = common_size (l, scale); + if (retval > 0) + error (["cauchy_rnd: ", ... + "lambda and sigma must be of common size or scalar"]); + endif + [r, c] = size (l); + else + usage ("cauchy_rnd (lambda, sigma [, r, c])"); + endif + + s = r * c; + l = reshape (l, 1, s); + scale = reshape (scale, 1, s); + + rnd = NaN * ones (1, s); + + k = find ((l > -Inf) & (l < Inf) & (scale > 0) & (scale < Inf)); + if (any (k)) + rnd(k) = l(k) - cot (pi * rand (1, length (k))) .* scale(k); + endif + + rnd = reshape (rnd, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/chisquare_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,50 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: chisquare_cdf (x, n) +## +## For each element of x, compute the cumulative distribution function +## (CDF) at x of the chisquare distribution with n degrees of freedom. + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: CDF of the chi-square distribution + +function cdf = chisquare_cdf (x, n) + + if (nargin != 2) + usage ("chisquare_cdf (x, n)"); + endif + + [retval, x, n] = common_size (x, n); + if (retval > 0) + error (["chisquare_cdf: ", ... + "x and n must be of common size or scalar"]); + endif + + cdf = gamma_cdf (x, n / 2, 1 / 2); + + ## should we really only allow for positive integer n? + k = find (n != round (n)); + if any (k) + fprintf (stderr, ... + "WARNING: n should be positive integer\n"); + [r, c] = size (x); + cdf = reshape (cdf, 1, r * c); + cdf(k) = NaN * ones (1, length (k)); + cdf = reshape (cdf, r, c); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/chisquare_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,49 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: chisquare_inv (x, n) +## +## For each element of x, compute the quantile (the inverse of the CDF) +## at x of the chisquare distribution with n degrees of freedom. + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: Quantile function of the chi-square distribution + +function inv = chisquare_inv (x, n) + + if (nargin != 2) + usage ("chisquare_inv (x, n)"); + endif + + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ("chisquare_inv: x and n must be of common size or scalar"); + endif + + inv = gamma_inv (x, n / 2, 1 / 2); + + ## Allow only for (positive) integer n. + k = find (n != round (n)); + if any (k) + fprintf (stderr, ... + "WARNING: n should be positive integer\n"); + [r, c] = size (x); + inv = reshape (inv, 1, r * c); + inv(k) = NaN * ones (1, length (k)); + inv = reshape (inv, r, c); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/chisquare_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,50 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: chisquare_pdf (x, n) +## +## For each element of x, compute the probability density function (PDF) +## at x of the chisquare distribution with k degrees of freedom. + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: PDF of the chi-sqaure distribution + +function pdf = chisquare_pdf (x, n) + + if (nargin != 2) + usage ("chisquare_pdf (x, n)"); + endif + + [retval, x, n] = common_size (x, n); + if (retval > 0) + error (["chisquare_pdf: ", ... + "x and n must be of common size or scalar"]); + endif + + pdf = gamma_pdf (x, n / 2, 1 / 2); + + ## should we really only allow for positive integer n? + k = find (n != round (n)); + if any (k) + fprintf (stderr, ... + "WARNING: n should be positive integer\n"); + [r, c] = size (x); + pdf = reshape (pdf, 1, r * c); + pdf(k) = NaN * ones (1, length (k)); + pdf = reshape (pdf, r, c); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/chisquare_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: chisquare_rnd (n [, r, c]) +## +## chisquare_rnd (n) returns a matrix of random samples from the +## chisquare distribution with n degrees of freedom. The size of the +## matrix is the size of n. +## +## chisquare_rnd (n, r, c) returns an r by c matrix of random samples +## from the chisquare distribution with n degrees of freedom. n must be +## a scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the chi-square distribution + +function rnd = chisquare_rnd (n, r, c) + + if (nargin == 3) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("chisquare_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("chisquare_rnd: c must be a positive integer"); + endif + [retval, n] = common_size (n, zeros (r, c)); + if (retval > 0) + error (strcat("chisquare_rnd: ", + "n must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin != 1) + usage ("chisquare_rnd (n [, r, c])"); + endif + + [r, c] = size (n); + s = r * c; + n = reshape (n, 1, s); + rnd = zeros (1, s); + + k = find (!(n > 0) | !(n < Inf) | !(n == round (n))); + if (any (k)) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find ((n > 0) & (n < Inf) & (n == round (n))); + if (any (k)) + rnd(k) = chisquare_inv (rand (1, length (k)), n(k)); + endif + + rnd = reshape (rnd, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/discrete_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,61 @@ +## Copyright (C) 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: discrete_cdf (X, V, P) +## +## For each element of X, compute the cumulative distribution function +## (CDF) at X of a univariate discrete distribution which assumes the +## values in V with probabilities P. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of a discrete distribution + +function cdf = discrete_cdf (X, V, P) + + if (nargin != 3) + usage ("discrete_cdf (X, V, P)"); + endif + + [r, c] = size (X); + + if (! is_vector (V)) + error ("discrete_cdf: V must be a vector"); + elseif (! is_vector (P) || (length (P) != length (V))) + error ("discrete_cdf: P must be a vector with length (V) elements"); + elseif (! (all (P >= 0) && any (P))) + error ("discrete_cdf: P must be a nonzero, nonnegative vector"); + endif + + n = r * c; + m = length (V); + X = reshape (X, n, 1); + V = reshape (V, 1, m); + P = reshape (P / sum (P), m, 1); + + cdf = zeros (n, 1); + k = find (isnan (X)); + if any (k) + cdf (k) = NaN * ones (length (k), 1); + endif + k = find (!isnan (X)); + if any (k) + n = length (k); + cdf (k) = ((X(k) * ones (1, m)) >= (ones (n, 1) * V)) * P; + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/discrete_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: discrete_inv (X, V, P) +## +## For each component of X, compute the quantile (the inverse of the +## CDF) at X of the univariate distribution which assumes the values in +## V with probabilities P. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of a discrete distribution + +function inv = discrete_inv (X, V, P) + + if (nargin != 3) + usage ("discrete_inv (X, V, P)"); + endif + + [r, c] = size (X); + + if (! is_vector (V)) + error ("discrete_inv: V must be a vector"); + elseif (! is_vector (P) || (length (P) != length (V))) + error ("discrete_inv: P must be a vector with length (V) elements"); + elseif (! (all (P >= 0) && any (P))) + error ("discrete_inv: P must be a nonzero, nonnegative vector"); + endif + + n = r * c; + X = reshape (X, 1, n); + m = length (V); + [V, ind] = sort (V); + s = reshape (cumsum (P / sum (P)), m, 1); + + inv = NaN * ones (n, 1); + if any (k = find (X == 0)) + inv(k) = -Inf * ones (1, length (k)); + endif + if any (k = find (X == 1)) + inv(k) = V(m) * ones (1, length (k)); + endif + if any (k = find ((X > 0) & (X < 1))) + n = length (k); + ## --FIXME-- + ## This does not work! + inv(k) = V(sum ((ones (m, 1) * X(k)) > (s * ones (1, n))) + 1); + endif + + inv = reshape (inv, r, c); + +endfunction + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/discrete_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,61 @@ +## Copyright (C) 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: discrete_pdf (X, V, P) +## +## For each element of X, compute the probability density function (PDF) +## at X of a univariate discrete distribution which assumes the values +## in V with probabilities P. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of a discrete distribution + +function pdf = discrete_pdf (X, V, P) + + if (nargin != 3) + usage ("discrete_pdf (X, V, P)"); + endif + + [r, c] = size (X); + + if (! is_vector (V)) + error ("discrete_pdf: V must be a vector"); + elseif (! is_vector (P) || (length (P) != length (V))) + error ("discrete_pdf: P must be a vector with length (V) elements"); + elseif (! (all (P >= 0) && any (P))) + error ("discrete_pdf: P must be a nonzero, nonnegative vector"); + endif + + n = r * c; + m = length (V); + X = reshape (X, n, 1); + V = reshape (V, 1, m); + P = reshape (P / sum (P), m, 1); + + pdf = zeros (n, 1); + k = find (isnan (X)); + if any (k) + pdf (k) = NaN * ones (length (k), 1); + endif + k = find (!isnan (X)); + if any (k) + n = length (k); + pdf (k) = ((X(k) * ones (1, m)) == (ones (n, 1) * V)) * P; + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/discrete_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,55 @@ +## Copyright (C) 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: discrete_rnd (N, V, P) +## +## Generate a random sample of size N from the univariate distribution +## which assumes the values in V with probabilities P. +## +## Currently, N must be a scalar. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from a discrete distribution + +function rnd = discrete_rnd (N, V, P) + + if (nargin != 3) + usage ("discrete_rnd (N, V, P)"); + endif + + if !is_scalar (N) + error ("discrete_rnd: N must be a scalar"); + endif + + if (! is_vector (V)) + error ("discrete_rnd: V must be a vector"); + elseif (! is_vector (P) || (length (P) != length (V))) + error ("discrete_rnd: P must be a vector with length (V) elements"); + elseif (! (all (P >= 0) && any (P))) + error ("discrete_rnd: P must be a nonzero, nonnegative vector"); + endif + + u = rand (1, N); + m = length (P); + s = reshape (cumsum (P / sum (P)), m, 1); + + rnd = V (1 + sum ((s * ones (1, N)) <= ((ones (m, 1) * u)))); + + if (prefer_column_vectors) + rnd = rnd'; + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/empirical_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,34 @@ +## Copyright (C) 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: empirical_cdf (X, DATA) +## +## For each element of X, compute the cumulative distribution function +## (CDF) at X of the empirical distribution obtained from the univariate +## sample DATA. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the empirical distribution + +function cdf = empirical_cdf (X, DATA) + + if (! is_vector (DATA)) + error ("empirical_cdf: DATA must be a vector"); + endif + + cdf = discrete_cdf (X, DATA, ones (size (DATA)) / length (DATA)); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/empirical_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,34 @@ +## Copyright (C) 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: empirical_inv (X, DATA) +## +## For each element of X, compute the quantile (the inverse of the CDF) +## at X of the empirical distribution obtained from the univariate +## sample DATA. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the empirical distribution + +function inv = empirical_inv (X, DATA) + + if (! is_vector (DATA)) + error ("empirical_inv: DATA must be a vector"); + endif + + inv = discrete_inv (X, DATA, ones (size (DATA)) / length (DATA)); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/empirical_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,34 @@ +## Copyright (C) 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: empirical_pdf (X, DATA) +## +## For each element of X, compute the probability density function (PDF) +## at X of the empirical distribution obtained from the univariate +## sample DATA. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the empirical distribution + +function pdf = empirical_pdf (X, DATA) + + if (! is_vector (DATA)) + error ("empirical_pdf: DATA must be a vector"); + endif + + pdf = discrete_pdf (X, DATA, ones (size (DATA)) / length (DATA)); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/empirical_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,33 @@ +## Copyright (C) 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: empirical_rnd (N, DATA) +## +## Generate a bootstrap sample of size N from the empirical distribution +## obtained from the univariate sample DATA. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Bootstrap samples from the empirical distribution + +function rnd = empirical_rnd (N, DATA) + + if (! is_vector (DATA)) + error ("empirical_rnd: DATA must be a vector"); + endif + + rnd = discrete_rnd (N, DATA, ones (size (DATA)) / length (DATA)); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/exponential_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,62 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: exponential_cdf (x, lambda) +## +## For each element of x, compute the cumulative distribution function +## (CDF) at x of the exponential distribution with parameter lambda. +## +## The arguments can be of common size or scalar. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the exponential distribution + +function cdf = exponential_cdf (x, l) + + if (nargin != 2) + usage ("exponential_cdf (x, lambda)"); + endif + + [retval, x, l] = common_size (x, l); + if (retval > 0) + error (["exponential_cdf: ", + "x and lambda must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + l = reshape (l, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | !(l > 0)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == Inf) & (l > 0)); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + k = find ((x > 0) & (x < Inf) & (l > 0)); + if any (k) + cdf (k) = 1 - exp (- l(k) .* x(k)); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/exponential_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,60 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: exponential_inv (x, lambda) +## +## For each element of x, compute the quantile (the inverse of the CDF) +## at x of the exponential distribution with parameter lambda. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the exponential distribution + +function inv = exponential_inv (x, l) + + if (nargin != 2) + usage ("exponential_inv (x, lambda)"); + endif + + [retval, x, l] = common_size (x, l); + if (retval > 0) + error (["exponential_inv: ", ... + "x and lambda must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + l = reshape (l, 1, s); + inv = zeros (1, s); + + k = find (!(l > 0) | (x < 0) | (x > 1) | isnan (x)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == 1) & (l > 0)); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < 1) & (l > 0)); + if any (k) + inv(k) = - log (1 - x(k)) ./ l(k); + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/exponential_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,55 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: exponential_pdf (x, lambda) +## +## For each element of x, compute the probability density function (PDF) +## of the exponential distribution with parameter lambda. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the exponential distribution + +function pdf = exponential_pdf (x, l) + + if (nargin != 2) + usage ("exponential_pdf (x, lambda)"); + endif + + [retval, x, l] = common_size (x, l); + if (retval > 0) + error (["exponential_pdf: ", ... + "x and lambda must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + l = reshape (l, 1, s); + pdf = zeros (1, s); + + k = find (!(l > 0) | isnan (x)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < Inf) & (l > 0)); + if any (k) + pdf(k) = l(k) .* exp (- l(k) .* x(k)); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/exponential_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: exponential_rnd (lambda [, r, c]) +## +## exponential_rnd (lambda) returns a matrix of random samples from the +## exponential distribution with parameter lambda. The size of the +## matrix is the size of lambda. +## +## exponential_rnd (lambda, r, c) returns an r by c matrix of random +## samples from the exponential distribution with parameter lambda, +## which must be a scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the exponential distribution + +function rnd = exponential_rnd (l, r, c) + + if (nargin == 3) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("exponential_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("exponential_rnd: c must be a positive integer"); + endif + [retval, l] = common_size (l, zeros (r, c)); + if (retval > 0) + error (strcat("exponential_rnd: ", + "lambda must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin != 1) + usage ("exponential_rnd (lambda [, r, c])"); + endif + + [r, c] = size (l); + s = r * c; + l = reshape (l, 1, s); + rnd = zeros (1, s); + + k = find (!(l > 0) | !(l < Inf)); + if (any (k)) + rnd(k) = NaN * ones (1, length (k)); + endif + k = find ((l > 0) & (l < Inf)); + if (any (k)) + rnd(k) = - log (1 - rand (1, length (k))) ./ l(k); + endif + + rnd = reshape (rnd, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/f_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,69 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: f_cdf (x, m, n) +## +## For each element of x, compute the CDF at x of the F distribution +## with m and n degrees of freedom, i.e., PROB( F(m,n) <= x ). + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the F distribution + +function cdf = f_cdf (x, m, n) + + if (nargin != 3) + usage ("f_cdf (x, m, n)"); + endif + + [retval, x, m, n] = common_size (x, m, n); + if (retval > 0) + error ("f_cdf: x, m and n must be of common size or scalar"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + m = reshape (m, 1, s); + n = reshape (n, 1, s); + cdf = zeros (1, s); + + k = find (!(m > 0) | !(n > 0) | isnan (x)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == Inf) & (m > 0) & (n > 0)); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + k = find ((x > 0) & (x < Inf) & (m > 0) & (n > 0)); + if any (k) + cdf(k) = 1 - betai (n(k) / 2, m(k) / 2, ... + 1 ./ (1 + m(k) .* x(k) ./ n(k))); + endif + + ## should we really only allow for positive integer m, n? + k = find ((m != round (m)) | (n != round (n))); + if any (k) + fprintf (stderr, ... + "WARNING: m and n should be positive integers\n"); + cdf(k) = NaN * ones (1, length (k)); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/f_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,69 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: f_inv (x, m, n) +## +## For each component of x, compute the quantile (the inverse of the +## CDF) at x of the F distribution with parameters m and n. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the F distribution + +function inv = f_inv (x, m, n) + + if (nargin != 3) + usage ("f_inv (x, m, n)"); + endif + + [retval, x, m, n] = common_size (x, m, n); + if (retval > 0) + error ("f_inv: x, m and n must be of common size or scalar"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + m = reshape (m, 1, s); + n = reshape (n, 1, s); + inv = zeros (1, s); + + k = find ((x < 0) | (x > 1) | isnan (x) | !(m > 0) | !(n > 0)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == 1) & (m > 0) & (n > 0)); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < 1) & (m > 0) & (n > 0)); + if any (k) + inv(k) = (1 ./ beta_inv (1 - x(k), n(k) / 2, m(k) / 2) - 1) ... + .* n(k) ./ m(k); + endif + + ## should we really only allow for positive integer m, n? + k = find ((m != round (m)) | (n != round (n))); + if any (k) + fprintf (stderr, ... + "WARNING: m and n should be positive integers\n"); + inv(k) = NaN * ones (1, length (k)); + endif + + inv = reshape (inv, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/f_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: f_pdf (x, m, n) +## +## For each element of x, compute the probability density function (PDF) +## at x of the F distribution with m and n degrees of freedom. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the F distribution + +function pdf = f_pdf (x, m, n) + + if (nargin != 3) + usage ("f_pdf (x, m, n)."); + endif + + [retval, x, m, n] = common_size (x, m, n); + if (retval > 0) + error ("f_pdf: x, m and n must be of common size or scalar"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + m = reshape (m, 1, s); + n = reshape (n, 1, s); + pdf = zeros (1, s); + + k = find (isnan (x) | !(m > 0) | !(n > 0)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < Inf) & (m > 0) & (n > 0)); + if any (k) + tmp = m(k) .* x(k) ./ n(k); + pdf(k) = exp ((m(k) / 2 - 1) .* log (tmp) ... + - ((m(k) + n(k)) / 2) .* log (1 + tmp)) ... + .* (m(k) ./ n(k)) ./ beta (m(k) / 2, n(k) / 2); + endif + + ## should we really only allow for positive integer m, n? + k = find ((m != round (m)) | (n != round (n))); + if any (k) + fprintf (stderr, ... + "WARNING: m and n should be positive integers\n"); + pdf(k) = NaN * ones (1, length (k)); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/f_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,74 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: f_rnd (m, n [, r, c]) +## +## f_rnd (m, n) returns a matrix of random samples from the F +## distribution with m and n degrees of freedom. The size of the matrix +## is the common size of m and n. +## +## f_rnd (m, n, r, c) returns an r by c matrix of random samples from +## the F distribution with m and n degrees of freedom. Both m and n +## must be scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the F distribution + +function rnd = f_rnd (m, n, r, c) + + if (nargin == 4) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("f_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("f_rnd: c must be a positive integer"); + endif + [retval, m, n] = common_size (m, n, zeros (r, c)); + if (retval > 0) + error (strcat("f_rnd: ", + "m and n must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin == 2) + [retval, m, n] = common_size (m, n); + if (retval > 0) + error ("f_rnd: m and n must be of common size or scalar"); + endif + else + usage ("f_rnd (m, n [, r, c])"); + endif + + [r, c] = size (m); + s = r * c; + m = reshape (m, 1, s); + n = reshape (n, 1, s); + rnd = zeros (1, s); + + k = find (!(m > 0) | !(m < Inf) | !(m == round (m)) | + !(n > 0) | !(n < Inf) | !(n == round (n))); + if (any (k)) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find ((m > 0) & (m < Inf) & (m == round (m)) & + (n > 0) & (n < Inf) & (n == round (n))); + if (any (k)) + rnd(k) = f_inv (rand (1, length (k)), m(k), n(k)); + endif + + rnd = reshape (rnd, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/gamma_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,55 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: gamma_cdf (x, a, b) +## +## For each element of x, compute the cumulative distribution function +## (CDF) at x of the Gamma distribution with parameters a and b. + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: CDF of the Gamma distribution + +function cdf = gamma_cdf (x, a, b) + + if (nargin != 3) + usage ("gamma_cdf (x, a, b)"); + endif + + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("gamma_cdf: x, a and b must be of common size or scalars"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, s, 1); + a = reshape (a, s, 1); + b = reshape (b, s, 1); + cdf = zeros (s, 1); + + k = find (!(a > 0) | !(b > 0) | isnan (x)); + if any (k) + cdf (k) = NaN * ones (length (k), 1); + endif + + k = find ((x > 0) & (a > 0) & (b > 0)); + if any (k) + cdf (k) = gammai (a(k), b(k) .* x(k)); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/gamma_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,84 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: gamma_inv (x, a, b) +## +## For each component of x, compute the quantile (the inverse of the +## CDF) at x of the Gamma distribution with parameters a and b. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Gamma distribution + +function inv = gamma_inv (x, a, b) + + if (nargin != 3) + usage ("gamma_inv (x, a, b)"); + endif + + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("gamma_inv: x, a and b must be of common size or scalars"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, s, 1); + a = reshape (a, s, 1); + b = reshape (b, s, 1); + inv = zeros (s, 1); + + k = find ((x < 0) | (x > 1) | isnan (x) | !(a > 0) | !(b > 0)); + if any (k) + inv (k) = NaN * ones (length (k), 1); + endif + + k = find ((x == 1) & (a > 0) & (b > 0)); + if any (k) + inv (k) = Inf * ones (length (k), 1); + endif + + k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0)); + if any (k) + a = a (k); + b = b (k); + x = x (k); + y = a ./ b; + l = find (x < eps); + if any (l) + y(l) = sqrt (eps) * ones (length (l), 1); + endif + + y_old = y; + for i = 1 : 100 + h = (gamma_cdf (y_old, a, b) - x) ./ gamma_pdf (y_old, a, b); + y_new = y_old - h; + ind = find (y_new <= eps); + if (any (ind)) + y_new (ind) = y_old (ind) / 10; + h = y_old - y_new; + endif + if (max (abs (h)) < sqrt (eps)) + break; + endif + y_old = y_new; + endfor + + inv (k) = y_new; + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/gamma_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,56 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: gamma_pdf (x, a, b) +## +## For each element of x, return the probability density function (PDF) +## at x of the Gamma distribution with parameters a and b. + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: PDF of the Gamma distribution + +function pdf = gamma_pdf (x, a, b) + + if (nargin != 3) + usage ("gamma_pdf (x, a, b)"); + endif + + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("gamma_pdf: x, a and b must be of common size or scalars"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, s, 1); + a = reshape (a, s, 1); + b = reshape (b, s, 1); + pdf = zeros (s, 1); + + k = find (!(a > 0) | !(b > 0) | isnan (x)); + if any (k) + pdf (k) = NaN * ones (length (k), 1); + endif + + k = find ((x > 0) & (a > 0) & (b > 0)); + if any (k) + pdf (k) = (b(k) .^ a(k)) .* (x(k) .^ (a(k) - 1)) ... + .* exp(-b(k) .* x(k)) ./ gamma (a(k)); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/gamma_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,72 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: gamma_rnd (a, b [, r, c]) +## +## gamma_rnd (a, b) returns a matrix of random samples from the Gamma +## distribution with parameters a and b. The size of the matrix is the +## common size of a and b. +## +## gamma_rnd (a, b, r, c) returns an r by c matrix of random samples +## from the Gamma distribution with parameters a and b. Both a and b +## must be scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the Gamma distribution + +function rnd = gamma_rnd (a, b, r, c) + + if (nargin == 4) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("gamma_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("gamma_rnd: c must be a positive integer"); + endif + [retval, a, b] = common_size (a, b, zeros (r, c)); + if (retval > 0) + error (strcat("gamma_rnd: ", + "a and b must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin == 2) + [retval, a, b] = common_size (a, b); + if (retval > 0) + error (["gamma_rnd: ", ... + "a and b must be of common size or scalar"]); + endif + else + usage ("gamma_rnd (a, b [, r, c])"); + endif + + [r, c] = size (a); + s = r * c; + a = reshape (a, 1, s); + b = reshape (b, 1, s); + rnd = zeros (1, s); + + k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)); + if (any (k)) + rnd(k) = NaN * ones (1, length (k)); + endif + k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf)); + if (any (k)) + rnd(k) = gamma_inv (rand (1, length (k)), a(k), b(k)); + endif + + rnd = reshape (rnd, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/geometric_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,61 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: geometric_cdf (x, p) +## +## For each element of x, compute the CDF at x of the geometric +## distribution with parameter p. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the geometric distribution + +function cdf = geometric_cdf (x, p) + + if (nargin != 2) + usage ("geometric_cdf (x, p)"); + endif + + [retval, x, p] = common_size (x, p); + if (retval > 0) + error (["geometric_cdf: ", ... + "x and p must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + p = reshape (p, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | !(p >= 0) | !(p <= 1)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == Inf) & (p >= 0) & (p <= 1)); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + k = find ((x >= 0) & (x < Inf) & (x == round (x)) ... + & (p > 0) & (p <= 1)); + if any (k) + cdf(k) = 1 - ((1 - p(k)) .^ (x(k) + 1)); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/geometric_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,61 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: geometric_inv (x, p) +## +## For each element of x, compute the quantile at x of the geometric +## distribution with parameter p. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the geometric distribution + +function inv = geometric_inv (x, p) + + if (nargin != 2) + usage ("geometric_inv (x, p)"); + endif + + [retval, x, p] = common_size (x, p); + if (retval > 0) + error (["geometric_inv: ", ... + "x and p must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + p = reshape (p, 1, s); + inv = zeros (1, s); + + k = find (!(x >= 0) | !(x <= 1) | !(p >= 0) | !(p <= 1)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == 1) & (p >= 0) & (p <= 1)); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < 1) & (p > 0) & (p <= 1)); + if any (k) + inv(k) = max (ceil (log (1 - x(k)) ./ log (1 - p(k))) - 1, ... + zeros (1, length (k))); + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/geometric_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,62 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: geometric_pdf (x, p) +## +## For each element of x, compute the probability density function (PDF) +## at x of the geometric distribution with parameter p. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the geometric distribution + +function pdf = geometric_pdf (x, p) + + if (nargin != 2) + usage ("geometric_pdf (x, p)"); + endif + + [retval, x, p] = common_size (x, p); + if (retval > 0) + error (["geometric_pdf: ", ... + "x and p must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + p = reshape (p, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | !(p >= 0) | !(p <= 1)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + ## Just for the fun of it ... + k = find ((x == Inf) & (p == 0)); + if any (k) + pdf(k) = ones (1, length (k)); + endif + + k = find ((x >= 0) & (x < Inf) & (x == round (x)) ... + & (p > 0) & (p <= 1)); + if any (k) + pdf(k) = p(k) .* ((1 - p(k)) .^ x(k)); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/geometric_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,71 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: geometric_rnd (p [, r, c]) +## +## geometric_rnd (p) returns a matrix of random samples from the +## geometric distribution with parameter p. The size of the matrix is +## the size of p. +## +## geometric_rnd (p, r, c) returns an r by c matrix of random samples +## from the geometric distribution with parameter p, which must be a +## scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the geometric distribution + +function rnd = geometric_rnd (p, r, c) + + if (nargin == 3) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("geometric_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("geometric_rnd: c must be a positive integer"); + endif + [retval, p] = common_size (p, zeros (r, c)); + if (retval > 0) + error (strcat("geometric_rnd: ", + "p must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin != 1) + usage ("geometric_rnd (p [, r, c])"); + endif + + [r, c] = size (p); + s = r * c; + p = reshape (p, 1, s); + rnd = zeros (1, s); + + k = find (!(p >= 0) | !(p <= 1)); + if (any (k)) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find (p == 0); + if (any (k)) + rnd(k) = Inf * ones (1, length (k)); + endif + + k = find ((p > 0) & (p < 1)); + if (any (k)) + rnd(k) = floor (log (rand (1, length (k))) ./ log (1 - p(k))); + endif + + rnd = reshape (rnd, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/hypergeometric_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,45 @@ +## Copyright (C) 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## Compute the cumulative distribution function (CDF) at x of the +## hypergeometric distribution with parameters m, t, and n. This is the +## probability of obtaining not more than x marked items when randomly +## drawing a sample of size n without replacement from a population of +## total size t containing m marked items. +## +## The parameters m, t, and n must positive integers with m and n not +## greater than t. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the hypergeometric distribution + +function cdf = hypergeometric_cdf (x, m, t, n) + + if (nargin != 4) + usage ("hypergeometrix_cdf (x, m, t, n)"); + endif + + if ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) | + (t != round (t)) | (n != round (n)) | (m > t) | (n > t)) + cdf = NaN * ones (size (x)) + else + cdf = discrete_cdf (x, 0 : n, hypergeometric_pdf (0 : n, m, t, n)); + endif + +endfunction + + + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/hypergeometric_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,39 @@ +## Copyright (C) 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## For each element of x, compute the quantile at x of the hypergeometric +## distribution with parameters m, t, and n. +## +## The parameters m, t, and n must positive integers with m and n not +## greater than t. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the hypergeometric distribution + +function inv = hypergeometric_inv (x, m, t, n) + + if (nargin != 4) + usage ("hypergeometric_inv (x, m, t, n)"); + endif + + if ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) | + (t != round (t)) | (n != round (n)) | (m > t) | (n > t)) + inv = NaN * ones (size (x)) + else + inv = discrete_inv (x, 0 : n, hypergeometric_pdf (0 : n, m, t, n)); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/hypergeometric_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,64 @@ +## Copyright (C) 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: hypergeometric_pdf (x, m, t, n) +## +## Compute the probability density function (PDF) at x of the +## hypergeometric distribution with parameters m, t, and n. This is the +## probability of obtaining x marked items when randomly drawing a +## sample of size n without replacement from a population of total size +## t containing m marked items. +## +## The arguments must be of common size or scalar. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the hypergeometric distribution + +function pdf = hypergeometric_pdf (x, m, t, n) + + if (nargin != 4) + usage ("hypergeometric_pdf (x, m, t, n)"); + endif + + [retval, x, m, t, n] = common_size (x, m, t, n); + if (retval > 0) + error (["hypergeometric_pdf: ", ... + "x, m, t, and n must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + m = reshape (m, 1, s); + t = reshape (t, 1, s); + n = reshape (n, 1, s); + pdf = zeros * ones (1, s); + ## everything in i1 gives NaN + i1 = ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) | + (t != round (t)) | (n != round (n)) | (m > t) | (n > t)); + ## everything in i2 gives 0 unless in i1 + i2 = ((x != round (x)) | (x < 0) | (x > m) | (n < x) | (n-x > t-m)); + k = find (i1); + if any (k) + pdf (k) = NaN * ones (size (k)); + endif + k = find (!i1 & !i2); + if any (k) + pdf (k) = (bincoeff (m(k), x(k)) .* bincoeff (t(k)-m(k), n(k)-x(k)) + ./ bincoeff (t(k), n(k))); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/hypergeometric_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,40 @@ +## Copyright (C) 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## Generate a random sample of size N from the hypergeometric +## distribution with parameters m, t, and n. +## +## The parameters m, t, and n must positive integers with m and n not +## greater than t. + +function rnd = hypergeometric_rnd (N, m, t, n) + + if (nargin != 4) + usage ("hypergeometric_rnd (N, m, t, n)"); + endif + + if ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) | + (t != round (t)) | (n != round (n)) | (m > t) | (n > t)) + if (prefer_column_vectors) + rnd = NaN * ones (N, 1) + else + rnd = NaN * ones (1, N) + endif + else + rnd = discrete_rnd (N, 0 : n, hypergeometric_pdf (0 : n, m, t, n)); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/kolmogorov_smirnov_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,64 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: kolmogorov_smirnov_cdf (x [, tol]) +## +## Returns the CDF at x of the Kolmogorov-Smirnov distribution, +## i.e. Q(x) = sum_{k=-\infty}^\infty (-1)^k exp(-2 k^2 x^2), x > 0. +## +## The optional tol specifies the precision up to which the series +## should be evaluated; the default is tol = eps. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the Kolmogorov-Smirnov distribution + +function cdf = kolmogorov_smirnov_cdf (x, tol) + + if (nargin < 1 || nargin > 2) + usage ("kolmogorov_smirnov_cdf (x [, tol])"); + endif + + if (nargin == 1) + tol = eps; + else + if (!is_scalar (tol) || !(tol > 0)) + error (["kolmogorov_smirnov_cdf: ", ... + "tol has to be a positive scalar."]); + endif + endif + + [nr, nc] = size(x); + if (min (nr, nc) == 0) + error ("kolmogorov_smirnov_cdf: x must not be empty."); + endif + + n = nr * nc; + x = reshape (x, 1, n); + cdf = zeros (1, n); + ind = find (x > 0); + if (length (ind) > 0) + y = x(ind); + K = ceil( sqrt( - log (tol) / 2 ) / min (y) ); + k = (1:K)'; + A = exp( - 2 * k.^2 * y.^2 ); + odd = find (rem (k, 2) == 1); + A(odd, :) = -A(odd, :); + cdf(ind) = 1 + 2 * sum (A); + endif + + cdf = reshape (cdf, nr, nc); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/laplace_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,53 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: laplace_cdf (x) +## +## For each element of x, compute the cumulative distribution function +## (CDF) at x of the Laplace distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the Laplace distribution + +function cdf = laplace_cdf (x) + + if (nargin != 1) + usage ("laplace_cdf (x)"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find (x == Inf); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + k = find ((x > -Inf) & (x < Inf)); + if any (k) + cdf(k) = (1 + sign (x(k)) .* (1 - exp (- abs (x(k))))) / 2; + endif + + cdf = reshape (cdf, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/laplace_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,54 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: laplace_inv (x) +## +## For each element of x, compute the quantile (the inverse of the CDF) +## at x of the Laplace distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Laplace distribution + +function inv = laplace_inv (x) + + if (nargin != 1) + usage ("laplace_inv (x)"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + inv = (-Inf) * ones (1, s); + + k = find (isnan (x) | (x < 0) | (x > 1)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find (x == 1); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < 1)); + if any (k) + inv(k) = (x(k) < 1/2) .* log (2 * x(k)) ... + - (x(k) > 1/2) .* log (2 * (1 - x(k))); + endif + + inv = reshape (inv, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/laplace_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,48 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: laplace_pdf (x) +## +## For each element of x, compute the probability density function (PDF) +## at x of the Laplace distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the Laplace distribution + +function pdf = laplace_pdf (x) + + if (nargin != 1) + usage ("laplace_pdf (x)"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + pdf = zeros (1, s); + + k = find (isnan (x)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x > -Inf) & (x < Inf)); + if any (k) + pdf(k) = exp (- abs (x(k))) / 2; + endif + + pdf = reshape (pdf, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/laplace_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,42 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: laplace_rnd (r, c) +## +## Return an r by c matrix of random numbers from the Laplace +## distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the Laplace distribution + +function rnd = laplace_rnd (r, c) + + if (nargin != 2) + usage ("laplace_rnd (r, c)"); + endif + + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("laplace_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("laplace_rnd: c must be a positive integer"); + endif + + tmp = rand (r, c); + rnd = ((tmp < 1/2) .* log (2 * tmp) + - (tmp > 1/2) .* log (2 * (1 - tmp))); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/logistic_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,33 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: logistic_cdf (x) +## +## For each component of x, compute the CDF at x of the logistic +## distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the logistic distribution + +function cdf = logistic_cdf (x) + + if (nargin != 1) + usage ("logistic_cdf (x)"); + endif + + cdf = 1 ./ (1 + exp (- x)); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/logistic_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,58 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: logistic_inv (x) +## +## For each component of x, compute the quantile (the inverse of the +## CDF) at x of the logistic distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the logistic distribution + +function inv = logistic_inv (x) + + if (nargin != 1) + usage ("logistic_inv (x)"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + inv = zeros (1, s); + + k = find ((x < 0) | (x > 1) | isnan (x)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find (x == 0); + if any (k) + inv(k) = (-Inf) * ones (1, length (k)); + endif + + k = find (x == 1); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < 1)); + if any (k) + inv (k) = - log (1 ./ x(k) - 1); + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/logistic_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,34 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: logistic_pdf (x) +## +## For each component of x, compute the PDF at x of the logistic +## distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the logistic distribution + +function pdf = logistic_pdf (x) + + if (nargin != 1) + usage ("logistic_pdf (x)"); + endif + + cdf = logistic_cdf (x); + pdf = cdf .* (1 - cdf); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/logistic_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,40 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: logistic_rnd (r, c) +## +## Return an r by c matrix of random numbers from the logistic +## distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the logistic distribution + +function rnd = logistic_rnd (r, c) + + if (nargin != 2) + usage ("logistic_rnd (r, c)"); + endif + + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("logistic_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("logistic_rnd: c must be a positive integer"); + endif + + rnd = - log (1 ./ rand (r, c) - 1); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/lognormal_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,77 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: lognormal_cdf (x [, a, v]) +## +## For each element of x, compute the cumulative distribution function +## (CDF) at x of the lognormal distribution with parameters a and v. If +## a random variable follows this distribution, its logarithm is +## normally distributed with mean log (a) and variance v. +## +## Default values are a = 1, v = 1. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the log normal distribution + +function cdf = lognormal_cdf (x, a, v) + + if !((nargin == 1) || (nargin == 3)) + usage ("lognormal_cdf (x [, a, v])"); + endif + + if (nargin == 1) + a = 1; + v = 1; + endif + + ## The following "straightforward" implementation unfortunately does + ## not work (because exp (Inf) -> NaN etc): + ## cdf = normal_cdf (log (x), log (a), v); + ## Hence ... + + [retval, x, a, v] = common_size (x, a, v); + if (retval > 0) + error (["lognormal_cdf: ", ... + "x, a and v must be of common size or scalars"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + a = reshape (a, 1, s); + v = reshape (v, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | !(a > 0) | !(a < Inf) ... + | !(v > 0) | !(v < Inf)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + k = find ((x > 0) & (x < Inf) & (a > 0) & (a < Inf) ... + & (v > 0) & (v < Inf)); + if any (k) + cdf(k) = stdnormal_cdf ((log (x(k)) - log (a(k))) ./ sqrt (v(k))); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/lognormal_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,77 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: lognormal_inv (x [, a, v]) +## +## For each element of x, compute the quantile (the inverse of the CDF) +## at x of the lognormal distribution with parameters a and v. If a +## random variable follows this distribution, its logarithm is normally +## distributed with mean log (a) and variance v. +## +## Default values are a = 1, v = 1. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the log normal distribution + +function inv = lognormal_inv (x, a, v) + + if !((nargin == 1) || (nargin == 3)) + usage ("lognormal_inv (x [, a, v])"); + endif + + if (nargin == 1) + a = 1; + v = 1; + endif + + ## The following "straightforward" implementation unfortunately does + ## not work (because exp (Inf) -> NaN): + ## inv = exp (normal_inv (x, log (a), v)); + ## Hence ... + + [retval, x, a, v] = common_size (x, a, v); + if (retval > 0) + error (["lognormal_inv: ", ... + "x, a and v must be of common size or scalars"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + a = reshape (a, 1, s); + v = reshape (v, 1, s); + inv = zeros (1, s); + + k = find (!(x >= 0) | !(x <= 1) | !(a > 0) | !(a < Inf) ... + | !(v > 0) | !(v < Inf)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == 1) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < 1) & (a > 0) & (a < Inf) ... + & (v > 0) & (v < Inf)); + if any (k) + inv(k) = a(k) .* exp (sqrt (v(k)) .* stdnormal_inv (x(k))); + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/lognormal_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,72 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: lognormal_pdf (x [, a, v]) +## +## For each element of x, compute the probability density function (PDF) +## at x of the lognormal distribution with parameters a and v. If a +## random variable follows this distribution, its logarithm is normally +## distributed with mean log (a) and variance v. +## +## Default values are a = 1, v = 1. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the log normal distribution + +function pdf = lognormal_pdf (x, a, v) + + if !((nargin == 1) || (nargin == 3)) + usage ("lognormal_pdf (x [, a, v])"); + endif + + if (nargin == 1) + a = 1; + v = 1; + endif + + ## The following "straightforward" implementation unfortunately does + ## not work for the special cases (Inf, ...) + ## pdf = (x > 0) ./ x .* normal_pdf (log (x), log (a), v); + ## Hence ... + + [retval, x, a, v] = common_size (x, a, v); + if (retval > 0) + error (["lognormal_pdf: ", ... + "x, a and v must be of common size or scalars"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + a = reshape (a, 1, s); + v = reshape (v, 1, s); + pdf = zeros (1, s); + + k = find (isnan (x) | !(a > 0) | !(a < Inf) ... + | !(v > 0) | !(v < Inf)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < Inf) & (a > 0) & (a < Inf) ... + & (v > 0) & (v < Inf)); + if any (k) + pdf(k) = normal_pdf (log (x(k)), log (a(k)), v(k)) ./ x(k); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/lognormal_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,73 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: lognormal_rnd (a, v [, r, c]) +## +## lognormal_rnd (a, v) returns a matrix of random samples from the +## lognormal distribution with parameters a and v. The size of the +## matrix is the common size of a and v. +## +## lognormal_rnd (a, v, r, c) returns an r by c matrix of random samples +## from the lognormal distribution with parameters a and v. Both a and v +## must be scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the log normal distribution + +function rnd = lognormal_rnd (a, v, r, c) + + if (nargin == 4) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("lognormal_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("lognormal_rnd: c must be a positive integer"); + endif + [retval, a, v] = common_size (a, v, zeros (r, c)); + if (retval > 0) + error (strcat("lognormal_rnd: ", + "a and v must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin == 2) + [retval, a, v] = common_size (a, v); + if (retval > 0) + error (strcat("lognormal_rnd: ", + "a and v must be of common size or scalar")); + endif + else + usage ("lognormal_rnd (a, v [, r, c])"); + endif + + [r, c] = size (a); + s = r * c; + a = reshape (a, 1, s); + v = reshape (v, 1, s); + rnd = zeros (1, s); + + k = find (!(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); + if any (k) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find ((a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + if any (k) + rnd(k) = a(k) .* exp (sqrt (v(k)) .* randn (1, length (k))); + endif + + rnd = reshape (rnd, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/normal_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,63 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: normal_cdf (x [, m, v]) +## +## For each element of x, compute the cumulative distribution function +## (CDF) at x of the normal distribution with mean m and variance v. +## +## Default values are m = 0, v = 1. + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: CDF of the normal distribution + +function cdf = normal_cdf (x, m, v) + + if !((nargin == 1) || (nargin == 3)) + usage ("normal_cdf (x [, m, v])"); + endif + + if (nargin == 1) + m = 0; + v = 1; + endif + + [retval, x, m, v] = common_size (x, m, v); + if (retval > 0) + error (["normal_cdf: ", ... + "x, m and v must be of common size or scalars"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + m = reshape (m, 1, s); + v = reshape (v, 1, s); + cdf = zeros (1, s); + + k = find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find (!isinf (m) & !isnan (m) & (v > 0) & (v < Inf)); + if any (k) + cdf(k) = stdnormal_cdf ((x(k) - m(k)) ./ sqrt (v(k))); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/normal_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,63 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: normal_inv (x [, m, v]) +## +## For each element of x, compute the quantile (the inverse of the CDF) +## at x of the normal distribution with mean m and variance v. +## +## Default values are m = 0, v = 1. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the normal distribution + +function inv = normal_inv (x, m, v) + + if !((nargin == 1) || (nargin == 3)) + usage ("normal_inv (x [, m, v])"); + endif + + if (nargin == 1) + m = 0; + v = 1; + endif + + [retval, x, m, v] = common_size (x, m, v); + if (retval > 0) + error (["normal_inv: ", ... + "x, m and v must be of common size or scalars"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + m = reshape (m, 1, s); + v = reshape (v, 1, s); + inv = zeros (1, s); + + k = find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find (!isinf (m) & !isnan (m) & (v > 0) & (v < Inf)); + if any (k) + inv(k) = m(k) + sqrt (v(k)) .* stdnormal_inv (x(k)); + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/normal_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,64 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: normal_pdf (x [, m, v]) +## +## For each element of x, compute the probability density function (PDF) +## at x of the normal distribution with mean m and variance v. +## +## Default values are m = 0, v = 1. + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: PDF of the normal distribution + +function pdf = normal_pdf (x, m, v) + + if !((nargin == 1) || (nargin == 3)) + usage ("normal_pdf (x [, m, v])"); + endif + + if (nargin == 1) + m = 0; + v = 1; + endif + + [retval, x, m, v] = common_size (x, m, v); + if (retval > 0) + error (["normal_pdf: ", ... + "x, m and v must be of common size or scalars"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + m = reshape (m, 1, s); + v = reshape (v, 1, s); + pdf = zeros (1, s); + + k = find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + k = find (!isinf (m) & !isnan (m) & (v > 0) & (v < Inf)); + if any (k) + pdf(k) = stdnormal_pdf ((x(k) - m(k)) ./ sqrt (v(k))) ... + ./ sqrt (v(k)); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/normal_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,72 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: normal_rnd (m, v [, r, c]) +## +## normal_rnd (m, v) returns a matrix of random samples from the normal +## distribution with parameters m and v. The size of the matrix is the +## common size of m and v. +## +## normal_rnd (m, v, r, c) returns an r by c matrix of random samples +## from the normal distribution with parameters m and v. Both m and v +## must be scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the normal distribution + +function rnd = normal_rnd (m, v, r, c) + + if (nargin == 4) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("normal_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("normal_rnd: c must be a positive integer"); + endif + [retval, m, v] = common_size (m, v, zeros (r, c)); + if (retval > 0) + error (strcat("normal_rnd: ", + "m and v must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin == 2) + [retval, m, v] = common_size (m, v); + if (retval > 0) + error ("normal_rnd: m and v must be of common size or scalar"); + endif + else + usage ("normal_rnd (m, v [, r, c])"); + endif + + [r, c] = size (m); + s = r * c; + m = reshape (m, 1, s); + v = reshape (v, 1, s); + rnd = zeros (1, s); + + k = find (isnan (m) | isinf (m) | !(v > 0) | !(v < Inf)); + if any (k) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find ((m > -Inf) & (m < Inf) & (v > 0) & (v < Inf)); + if any (k) + rnd(k) = m(k) + sqrt (v(k)) .* randn (1, length (k)); + endif + + rnd = reshape (rnd, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/pascal_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,82 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: pascal_cdf (x, n, p) +## +## For each element of x, compute the CDF at x of the Pascal (negative +## binomial) distribution with parameters n and p. +## +## The number of failures in a Bernoulli experiment with success +## probability p before the n-th success follows this distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the Pascal (negative binomial) distribution + +function cdf = pascal_cdf (x, n, p) + + if (nargin != 3) + usage ("pascal_cdf (x, n, p)"); + endif + + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error (["pascal_cdf: ", ... + "x, n and p must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + n = reshape (n, 1, s); + p = reshape (p, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n)) ... + | (p < 0) | (p > 1)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n)) ... + & (p >= 0) & (p <= 1)); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0) ... + & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); + if any (k) + ## Does anyone know a better way to do the summation? + m = zeros (1, length (k)); + x = floor (x(k)); + n = n(k); + p = p(k); + y = cdf(k); + while (1) + l = find (m <= x); + if any (l) + y(l) = y(l) + pascal_pdf (m(l), n(l), p(l)); + m(l) = m(l) + 1; + else + break; + endif + endwhile + cdf(k) = y; + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/pascal_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,81 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: pascal_inv (x, n, p) +## +## For each element of x, compute the quantile at x of the Pascal +## (negative binomial) distribution with parameters n and p. +## +## The number of failures in a Bernoulli experiment with success +## probability p before the n-th success follows this distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Pascal distribution + +function inv = pascal_inv (x, n, p) + + if (nargin != 3) + usage ("pascal_inv (x, n, p)"); + endif + + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error (["pascal_inv: ", ... + "x, n and p must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + n = reshape (n, 1, s); + p = reshape (p, 1, s); + inv = zeros (1, s); + + k = find (isnan (x) | (x < 0) | (x > 1) | (n < 1) | (n == Inf) ... + | (n != round (n)) | (p < 0) | (p > 1)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == 1) & (n > 0) & (n < Inf) & (n == round (n)) ... + & (p >= 0) & (p <= 1)); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + k = find ((x >= 0) & (x < 1) & (n > 0) & (n < Inf) ... + & (n == round (n)) & (p > 0) & (p <= 1)); + if any (k) + x = x(k); + n = n(k); + p = p(k); + m = zeros (1, length (k)); + s = p .^ n; + while (1) + l = find (s < x); + if any (l) + m(l) = m(l) + 1; + s(l) = s(l) + pascal_pdf (m(l), n(l), p(l)); + else + break; + endif + endwhile + inv(k) = m; + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/pascal_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,70 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: pascal_pdf (x, n, p) +## +## For each element of x, compute the probability density function (PDF) +## at x of the Pascal (negative binomial) distribution with parameters n +## and p. +## +## The number of failures in a Bernoulli experiment with success +## probability p before the n-th success follows this distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the Pascal (negative binomial) distribution + +function pdf = pascal_pdf (x, n, p) + + if (nargin != 3) + usage ("pascal_pdf (x, n, p)"); + endif + + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error (["pascal_pdf: ", ... + "x, n and p must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + n = reshape (n, 1, s); + p = reshape (p, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n)) ... + | (p < 0) | (p > 1)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + ## Just for the fun of it ... + k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n)) ... + & (p == 0)); + if any (k) + pdf(k) = ones (1, length (k)); + endif + + k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0) ... + & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); + if any (k) + pdf(k) = bincoeff (-n(k), x(k)) .* (p(k) .^ n(k)) ... + .* ((p(k) - 1) .^ x(k)); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/pascal_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,78 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: pascal_rnd (n, p [, r, c]) +## +## pascal_rnd (n, p) returns a matrix of random samples from the Pascal +## (negative binomial) distribution with parameters n and p. The size of +## the matrix is the common size of n and p. +## +## pascal_rnd (n, p, r, c) returns an r by c matrix of random samples +## from the Pascal distribution with parameters n and p. Both n and p +## must be scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the Pascal distribution + +function rnd = pascal_rnd (n, p, r, c) + + if (nargin == 4) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("pascal_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("pascal_rnd: c must be a positive integer"); + endif + [retval, n, p] = common_size (n, p, zeros (r, c)); + if (retval > 0) + error (strcat("pascal_rnd: ", + "n and p must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin == 2) + [retval, n, p] = common_size (n, p); + if (retval > 0) + error ("pascal_rnd: n and p must be of common size or scalar"); + endif + else + usage ("pascal_rnd (n, p [, r, c])"); + endif + + [r, c] = size (n); + s = r * c; + n = reshape (n, 1, s); + p = reshape (p, 1, s); + rnd = zeros (1, s); + + k = find (!(n > 0) | !(n < Inf) | !(n == round (n)) ... + | !(p <= 0) | !(p >= 1)); + if (any (k)) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find ((n > 0) & (n < Inf) & (n == round (n)) ... + & (p >= 0) & (p <= 1)); + if (any (k)) + N = max (n(k)); + L = length (k); + tmp = floor (log (rand (N, L)) ./ (ones (N, 1) * log (1 - p(k)))); + ind = (1 : N)' * ones (1, L); + rnd(k) = sum (tmp .* (ind <= ones (N, 1) * n(k))); + endif + + rnd = reshape (rnd, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/poisson_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,60 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: poisson_cdf (x, lambda) +## +## For each element of x, compute the cumulative distribution function +## (CDF) at x of the Poisson distribution with parameter lambda. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the Poisson distribution + +function cdf = poisson_cdf (x, l) + + if (nargin != 2) + usage ("poisson_cdf (x, lambda)"); + endif + + [retval, x, l] = common_size (x, l); + if (retval > 0) + error (["poisson_cdf: ", + "x and lambda must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + l = reshape (l, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | !(l > 0)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == Inf) & (l > 0)); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + k = find ((x >= 0) & (x < Inf) & (l > 0)); + if any (k) + cdf(k) = 1 - gammai (floor (x(k)) + 1, l(k)); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/poisson_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,69 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: poisson_inv (x, lambda) +## +## For each component of x, compute the quantile (the inverse of the +## CDF) at x of the Poisson distribution with parameter lambda. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Poisson distribution + +function inv = poisson_inv (x, l) + + if (nargin != 2) + usage ("poisson_inv (x, lambda)"); + endif + + [retval, x, l] = common_size (x, l); + if (retval > 0) + error (["poisson_inv: ", ... + "x and lambda must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + l = reshape (l, 1, s); + inv = zeros (1, s); + + k = find ((x < 0) | (x > 1) | isnan (x) | !(l > 0)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == 1) & (l > 0)); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < 1) & (l > 0)); + if any (k) + cdf = exp (-l(k)); + while (1) + m = find (cdf < x(k)); + if any (m) + inv(k(m)) = inv(k(m)) + 1; + cdf(m) = cdf(m) + poisson_pdf (inv(k(m)), l(k(m))); + else + break; + endif + endwhile + endif + + inv = reshape (inv, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/poisson_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,55 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: poisson_pdf (x, lambda) +## +## For each element of x, compute the probability density function (PDF) +## at x of the poisson distribution with parameter lambda. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the Poisson distribution + +function pdf = poisson_pdf (x, l) + + if (nargin != 2) + usage ("poisson_pdf (x, lambda)"); + endif + + [retval, x, l] = common_size (x, l); + if (retval > 0) + error (["poisson_pdf: ", ... + "x and lambda must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + l = reshape (l, 1, s); + pdf = zeros (1, s); + + k = find (!(l > 0) | isnan (x)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (l > 0)); + if any (k) + pdf(k) = exp (x(k) .* log (l(k)) - l(k) - lgamma (x(k) + 1)); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/poisson_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,80 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: poisson_rnd (lambda [, r, c]) +## +## poisson_rnd (lambda) returns a matrix of random samples from the +## Poisson distribution with parameter lambda. The size of the matrix +## is the size of lambda. +## +## poisson_rnd (lambda, r, c) returns an r by c matrix of random samples +## from the Poisson distribution with parameter lambda, which must be a +## scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the Poisson distribution + +function rnd = poisson_rnd (l, r, c) + + if (nargin == 3) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("poisson_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("poisson_rnd: c must be a positive integer"); + endif + [retval, l] = common_size (l, zeros (r, c)); + if (retval > 0) + error (strcat("poisson_rnd: ", + "lambda must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin != 1) + usage ("poisson_rnd (lambda [, r, c])"); + endif + + [r, c] = size (l); + s = r * c; + l = reshape (l, 1, s); + rnd = zeros (1, s); + + k = find (!(l > 0) | !(l < Inf)); + if (any (k)) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find ((l > 0) & (l < Inf)); + if (any (k)) + l = l(k); + len = length (k); + num = zeros (1, len); + sum = - log (1 - rand (1, len)) ./ l; + while (1) + ind = find (sum < 1); + if (any (ind)) + sum(ind) = (sum(ind) + - log (1 - rand (1, length (ind))) ./ l(ind)); + num(ind) = num(ind) + 1; + else + break; + endif + endwhile + rnd(k) = num; + endif + + rnd = reshape (rnd, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/stdnormal_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,42 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: stdnormal_cdf (x) +## +## For each component of x, compute the CDF of the standard normal +## distribution at x. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the standard normal distribution + +function cdf = stdnormal_cdf (x) + + if (nargin != 1) + usage ("stdnormal_cdf (x)"); + endif + + [r_x, c_x] = size (x); + if (r_x * c_x == 0) + error ("stdnormal_cdf: x must not be empty."); + endif + + cdf = ( ones (r_x, c_x) + erf (x / sqrt(2)) ) / 2; + +endfunction + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/stdnormal_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,33 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: stdnormal_inv (x) +## +## For each component of x, compute compute the quantile (the inverse of +## the CDF) at x of the standard normal distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the standard normal distribution + +function inv = stdnormal_inv (x) + + if (nargin != 1) + usage ("stdnormal_inv (x)"); + endif + + inv = sqrt (2) * erfinv (2 * x - 1); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/stdnormal_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,48 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: stdnormal_pdf (x) +## +## For each element of x, compute the probability density function (PDF) +## of the standard normal distribution at x. + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: PDF of the standard normal distribution + +function pdf = stdnormal_pdf (x) + + if (nargin != 1) + usage ("stdnormal_pdf (x)"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + pdf = zeros (1, s); + + k = find (isnan (x)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + k = find (!isinf (x)); + if any (k) + pdf (k) = (2 * pi)^(- 1/2) * exp( - x(k) .^ 2 / 2); + endif + + pdf = reshape (pdf, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/stdnormal_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,40 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: stdnormal_rnd (r, c) +## +## Return an r by c matrix of random numbers from the standard normal +## distribution. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the standard normal distribution + +function rnd = stdnormal_rnd (r, c) + + if (nargin != 2) + usage ("stdnormal_rnd (r, c)"); + endif + + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("stdnormal_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("stdnormal_rnd: c must be a positive integer"); + endif + + rnd = randn (r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/t_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,71 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: t_cdf (x, n) +## +## For each element of x, compute the CDF at x of the Student (t) +## distribution with n degrees of freedom, i.e., PROB( t(n) <= x ). + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the t distribution + +function cdf = t_cdf (x, n) + + if (nargin != 2) + usage ("t_cdf (x, n)"); + endif + + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ("t_cdf: x and n must be of common size or scalar"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + n = reshape (n, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | !(n > 0)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == Inf) & (n > 0)); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + k = find ((x > -Inf) & (x < Inf) & (n > 0)); + if any (k) + cdf(k) = betai (n(k) / 2, 1 / 2, 1 ./ (1 + x(k) .^ 2 ./ n(k))) / 2; + ind = find (x(k) > 0); + if any (ind) + cdf(k(ind)) = 1 - cdf(k(ind)); + endif + endif + + ## should we really only allow for positive integer n? + k = find (n != round (n)); + if any (k) + fprintf (stderr, ... + "WARNING: n should be positive integer\n"); + cdf(k) = NaN * ones (1, length (k)); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/t_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,83 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: t_inv (x, n) +## +## For each component of x, compute the quantile (the inverse of the +## CDF) at x of the t (Student) distribution with parameter n. + +## For very large n, the "correct" formula does not really work well, +## and the quantiles of the standard normal distribution are used +## directly. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the t distribution + +function inv = t_inv (x, n) + + if (nargin != 2) + usage ("t_inv (x, n)"); + endif + + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ("t_inv: x and n must be of common size or scalar"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + n = reshape (n, 1, s); + inv = zeros (1, s); + + k = find ((x < 0) | (x > 1) | isnan (x) | !(n > 0)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find ((x == 0) & (n > 0)); + if any (k) + inv(k) = (-Inf) * ones (1, length (k)); + endif + + k = find ((x == 1) & (n > 0)); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < 1) & (n > 0) & (n < 10000)); + if any (k) + inv(k) = sign (x(k) - 1/2) .* sqrt (n(k) .* (1 ... + ./ beta_inv (2 * min (x(k), 1 - x(k)), n(k) / 2, 1 / 2) - 1)); + endif + + ## For large n, use the quantiles of the standard normal + k = find ((x > 0) & (x < 1) & (n >= 10000)); + if any (k) + inv(k) = stdnormal_inv (x(k)); + endif + + ## should we really only allow for positive integer n? + k = find (n != round (n)); + if any (k) + fprintf (stderr, ... + "WARNING: n should be positive integer\n"); + inv(k) = NaN * ones (1, length (k)); + endif + + inv = reshape (inv, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/t_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,63 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: t_pdf (x, n) +## +## For each element of x, compute the probability density function (PDF) +## at x of the t (Student) distribution with n degrees of freedom. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the t distribution + +function pdf = t_pdf (x, n) + + if (nargin != 2) + usage ("t_pdf (x, n)"); + endif + + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ("t_pdf: x and n must be of common size or scalar"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + n = reshape (n, 1, s); + pdf = zeros (1, s); + + k = find (isnan (x) | !(n > 0) | !(n < Inf)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + k = find (!isinf (x) & !isnan (x) & (n > 0) & (n < Inf)); + if any (k) + pdf(k) = exp (- (n(k) + 1) .* log (1 + x(k) .^ 2 ./ n(k)) / 2) ... + ./ (sqrt (n(k)) .* beta (n(k) / 2, 1 / 2)); + endif + + ## should we really only allow for positive integer n? + k = find (n != round (n)); + if any (k) + fprintf (stderr, ... + "WARNING: n should be positive integer\n"); + pdf(k) = NaN * ones (1, length (k)); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/t_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: t_rnd (n [, r, c]) +## +## t_rnd (n) returns a matrix of random samples from the t (Student) +## distribution with n degrees of freedom. The size of the matrix is +## the size of n. +## +## t_rnd (n, r, c) returns an r by c matrix of random samples from the t +## distribution with n degrees of freedom. n must be a scalar or of +## size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the t distribution + +function rnd = t_rnd (n, r, c) + + if (nargin == 3) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("t_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("t_rnd: c must be a positive integer"); + endif + [retval, n] = common_size (n, zeros (r, c)); + if (retval > 0) + error (strcat("t_rnd: ", + "n must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin != 1) + usage ("t_rnd (n [, r, c])"); + endif + + [r, c] = size (n); + s = r * c; + n = reshape (n, 1, s); + rnd = zeros (1, s); + + k = find (!(n > 0) | !(n < Inf) | !(n == round (n))); + if any (k) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find ((n > 0) & (n < Inf) & (n == round (n))); + if any (k) + rnd(k) = t_inv (rand (1, length (k)), n(k)); + endif + + rnd = reshape (rnd, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/uniform_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,68 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: uniform_cdf (x [, a, b]) +## +## Returns the CDF at x of the uniform distribution on [a, b], i.e., +## PROB( uniform(a,b) <= x ). +## +## Default values are a = 0, b = 1. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the uniform distribution + +function cdf = uniform_cdf (x, a, b) + + if !(nargin == 1 || nargin == 3) + usage ("uniform_cdf (x [, a, b])"); + endif + + if (nargin == 1) + a = 0; + b = 1; + endif + + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error (["uniform_cdf: ", ... + "x, a and b must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + a = reshape (a, 1, s); + b = reshape (b, 1, s); + cdf = zeros (1, s); + + k = find (isnan (x) | !(a < b)); + if any (k) + cdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x >= b) & (a < b)); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + k = find ((x > a) & (x < b)); + if any (k) + cdf(k) = (x(k) < b(k)) .* (x(k) - a(k)) ./ (b(k) - a(k)); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/uniform_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,62 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: uniform_inv (x [, a, b]) +## +## For each element of x, compute the quantile (the inverse of the CDF) +## at x of the uniform distribution on [a, b]. +## +## Default values are a = 0, b = 1. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the uniform distribution + +function inv = uniform_inv (x, a, b) + + if !(nargin == 1 || nargin == 3) + usage ("uniform_inv (x [, a, b])"); + endif + + if (nargin == 1) + a = 0; + b = 1; + endif + + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("uniform_inv: x, a and b must be of common size or scalars"); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + a = reshape (a, 1, s); + b = reshape (b, 1, s); + inv = zeros (1, s); + + k = find ((x < 0) | (x > 1) | isnan (x) | !(a < b)); + if any (k) + inv(k) = NaN * ones (1, length (k)); + endif + + k = find ((x >= 0) & (x <= 1) & (a < b)); + if any (k) + inv(k) = a(k) + x(k) .* (b(k) - a(k)); + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/uniform_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,63 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: uniform_pdf (x [, a, b]) +## +## For each element of x, compute the PDF at x of the uniform +## distribution on [a, b]. +## +## Default values are a = 0, b = 1. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the uniform distribution + +function pdf = uniform_pdf (x, a, b) + + if !(nargin == 1 || nargin == 3) + usage ("uniform_pdf (x [, a, b])"); + endif + + if (nargin == 1) + a = 0; + b = 1; + endif + + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error (["uniform_pdf: ", ... + "x, a and b must be of common size or scalars"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + a = reshape (a, 1, s); + b = reshape (b, 1, s); + pdf = zeros (1, s); + + k = find (isnan (x) | !(a < b)); + if any (k) + pdf(k) = NaN * ones (1, length (k)); + endif + + k = find ((x > a) & (x < b)); + if any (k) + pdf(k) = 1 ./ (b(k) - a(k)); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/uniform_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,72 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: uniform_rnd (a, b [, r, c]) +## +## uniform_rnd (a, b) returns a matrix of random samples from the +## uniform distribution on [a, b]. The size of the matrix is the common +## size of a and b. +## +## uniform_rnd (a, b, r, c) returns an r by c matrix of random samples +## from the uniform distribution on [a, b]. Both a and b must be scalar +## or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the uniform distribution + +function rnd = uniform_rnd (a, b, r, c) + + if (nargin == 4) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("uniform_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("uniform_rnd: c must be a positive integer"); + endif + [retval, a, b] = common_size (a, b, zeros (r, c)); + if (retval > 0) + error (strcat("uniform_rnd: ", + "a and b must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin == 2) + [retval, a, b] = common_size (a, b); + if (retval > 0) + error ("uniform_rnd: a and b must be of common size or scalar"); + endif + else + usage ("uniform_rnd (a, b [, r, c])"); + endif + + [r, c] = size (a); + s = r * c; + a = reshape (a, 1, s); + b = reshape (b, 1, s); + rnd = zeros (1, s); + + k = find (!(-Inf < a) | !(a < b) | !(b < Inf)); + if any (k) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find ((-Inf < a) & (a < b) & (b < Inf)); + if any (k) + rnd(k) = a(k) + (b(k) - a(k)) .* rand (1, length (k)); + endif + + rnd = reshape (rnd, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/weibull_cdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: weibull_cdf (x, alpha, sigma) +## +## Compute the cumulative distribution function (CDF) at x of the +## Weibull distribution with shape parameter alpha and scale parameter +## sigma, which is 1 - exp(-(x/sigma)^alpha), x >= 0. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the Weibull distribution + +function cdf = weibull_cdf (x, shape, scale) + + if (nargin != 3) + usage ("weibull_cdf (x, alpha, sigma)"); + endif + + [retval, x, shape, scale] = common_size (x, shape, scale); + if (retval > 0) + error (["weibull_cdf: ", ... + "x, alpha and sigma must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + shape = reshape (shape, 1, s); + scale = reshape (scale, 1, s); + + cdf = NaN * ones (1, s); + + ok = ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); + + k = find ((x <= 0) & ok); + if any (k) + cdf(k) = zeros (1, length (k)); + endif + + k = find ((x > 0) & (x < Inf) & ok); + if any (k) + cdf(k) = 1 - exp (- (x(k) ./ scale(k)) .^ shape(k)); + endif + + k = find ((x == Inf) & ok); + if any (k) + cdf(k) = ones (1, length (k)); + endif + + cdf = reshape (cdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/weibull_inv.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,63 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: weibull_inv (x, lambda, alpha) +## +## Compute the quantile (the inverse of the CDF) at x of the Weibull +## distribution with shape parameter alpha and scale parameter sigma. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Weibull distribution + +function inv = weibull_inv (x, shape, scale) + + if (nargin != 3) + usage ("weibull_inv (x, alpha, sigma)"); + endif + + [retval, x, shape, scale] = common_size (x, shape, scale); + if (retval > 0) + error (["weibull_inv: ", ... + "x, alpha and sigma must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + shape = reshape (shape, 1, s); + scale = reshape (scale, 1, s); + + inv = NaN * ones (1, s); + ok = ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); + + k = find ((x == 0) & ok); + if any (k) + inv(k) = -Inf * ones (1, length (k)); + endif + + k = find ((x > 0) & (x < 1) & ok); + if any (k) + inv(k) = scale(k) .* (- log (1 - x(k))) .^ (1 ./ shape(k)); + endif + + k = find ((x == 1) & ok); + if any (k) + inv(k) = Inf * ones (1, length (k)); + endif + + inv = reshape (inv, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/weibull_pdf.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,63 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: weibull_pdf (x, alpha, sigma) +## +## Compute the probability density function (PDF) at x of the Weibull +## distribution with shape parameter alpha and scale parameter sigma +## which is given by +## alpha * sigma^(-alpha) * x^(alpha-1) * exp(-(x/sigma)^alpha) +## for x > 0. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the Weibull distribution + +function pdf = weibull_pdf (x, shape, scale) + + if (nargin != 3) + usage ("weibull_pdf (x, alpha, sigma)"); + endif + + [retval, x, shape, scale] = common_size (x, shape, scale); + if (retval > 0) + error (["weibull_pdf: ", ... + "x, alpha and sigma must be of common size or scalar"]); + endif + + [r, c] = size (x); + s = r * c; + x = reshape (x, 1, s); + shape = reshape (shape, 1, s); + scale = reshape (scale, 1, s); + + pdf = NaN * ones (1, s); + ok = ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); + + k = find ((x > -Inf) & (x <= 0) & ok); + if any (k) + pdf(k) = zeros (1, length (k)); + endif + + k = find ((x > 0) & (x < Inf) & ok); + if any (k) + pdf(k) = (shape(k) .* (scale(k) .^ shape(k)) + .* (x(k) .^ (shape(k) - 1)) + .* exp(- (x(k) ./ scale(k)) .^ shape(k))); + endif + + pdf = reshape (pdf, r, c); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/weibull_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,69 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: weibull_rnd (alpha, sigma [, r, c]) +## +## weibull_rnd (alpha, sigma) returns a matrix of random samples from +## the Weibull distribution with parameters alpha and sigma. The size of +## the matrix is the common size of alpha and sigma. +## +## weibull_rnd (alpha, sigma, r, c) returns an r by c matrix of random +## samples from the Weibull distribution with parameters alpha and sigma +## which must be scalar or of size r by c. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the Weibull distribution + +function rnd = weibull_rnd (shape, scale, r, c) + + if (nargin == 4) + if ( !(is_scalar (r) && (r > 0) && (r == round (r))) ) + error ("weibull_rnd: r must be a positive integer"); + endif + if ( !(is_scalar (c) && (c > 0) && (c == round (c))) ) + error ("weibull_rnd: c must be a positive integer"); + endif + [retval, shape, scale] = common_size (shape, scale, zeros (r, c)); + if (retval > 0) + error (strcat("weibull_rnd: ", + "alpha and sigma must be scalar or of size ", + sprintf ("%d by %d", r, c))); + endif + elseif (nargin == 2) + [retval, shape, scale] = common_size (shape, scale); + if (retval > 0) + error (strcat("weibull_rnd: ", + "alpha and sigma must be of common size or scalar")); + endif + else + usage ("weibull_rnd (alpha, sigma [, r, c])"); + endif + + [r, c] = size (shape); + s = r * c; + shape = reshape (shape, 1, s); + scale = reshape (scale, 1, s); + + rnd = NaN * ones (1, s); + k = find ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); + if any (k) + rnd(k) = (scale(k) + .* (- log (1 - rand (1, length (k)))) .^ (1 ./ shape(k))); + endif + + rnd = reshape (rnd, r, c); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/wiener_rnd.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,46 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: wiener_rnd (t [, d [,n]]) +## +## Returns a simulated realization of the d-dimensional Wiener Process +## on the interval [0,t]. If d is omitted, d=1 is used. The first +## column of the return matrix contains time, the remaining columns +## contain the Wiener process. +## +## The optional parameter n gives the number of summands used for +## simulating the process over an interval of length 1. If n is +## omitted, n=1000 is used. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Simulate a Wiener process + +function retval = wiener_rnd (t, d, n) + + if (nargin == 1) + d = 1; + n = 1000; + elseif (nargin == 2) + n = 1000; + elseif (nargin > 3) + usage ("wiener_rnd (t [, d [,n]])"); + endif + + retval = randn (n * t, d); + retval = cumsum (retval) / sqrt (n); + + retval = [(1: n*t)' / n retval]; +endfunction
--- a/scripts/statistics/gls.m Mon Oct 19 22:37:57 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,69 +0,0 @@ -## Copyright (C) 1996, 1997 John W. Eaton -## -## 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 2, 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, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: [BETA, v [,R]] = gls (Y, X, O) -## -## Generalized Least Squares (GLS) estimation for the multivariate model -## -## Y = X*B + E, mean(E) = 0, cov(vec(E)) = (s^2)*O -## -## with Y ... T x p As usual, each row of Y and X is an observation -## X ... T x k and each column a variable. -## B ... k x p -## E ... T x p -## O ... Tp x Tp. -## -## BETA is the GLS estimator for B. -## v is the GLS estimator for s^2. -## R = Y - X*BETA is the matrix of GLS residuals. - -## Author: Teresa Twaroch <twaroch@ci.tuwien.ac.at> -## Created: May 1993 -## Adapted-By: jwe - -function [BETA, v, R] = gls (Y, X, O) - - if (nargin != 3) - usage ("[BETA, v [, R]] = gls (Y, X, O)"); - endif - - [rx, cx] = size (X); - [ry, cy] = size (Y); - if (rx != ry) - error ("gls: incorrect matrix dimensions"); - endif - - O = O^(-1/2); - Z = kron (eye (cy), X); - Z = O * Z; - Y1 = O * reshape (Y, ry*cy, 1); - U = Z' * Z; - r = rank (U); - - if (r == cx*cy) - B = inv (U) * Z' * Y1; - else - B = pinv (Z) * Y1; - endif - - BETA = reshape (B, cx, cy); - R = Y - X * BETA; - v = (reshape (R, ry*cy, 1))' * (O^2) * reshape (R, ry*cy, 1) / (rx*cy - r); - -endfunction
--- a/scripts/statistics/kurtosis.m Mon Oct 19 22:37:57 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,59 +0,0 @@ -## Copyright (C) 1996, 1997 John W. Eaton -## -## 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 2, 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, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: kurtosis (x) -## -## If x is a vector of length N, return the kurtosis -## -## kurtosis(x) = N^(-1) std(x)^(-4) SUM_i (x(i)-mean(x))^4 - 3 -## -## of x. -## -## If x is a matrix, return a row vector containing the kurtosis for each -## column. - -## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> -## Created: 29 July 1994 -## Adapted-By: jwe - -function retval = kurtosis (x) - - if (nargin != 1) - usage ("kurtosis (x)"); - endif - - if (is_vector (x)) - x = x - mean (x); - if (! any (x)) - retval = 0; - else - retval = sum (x .^ 4) / (length (x) * std (x) ^ 4) - 3; - endif - elseif (is_matrix (x)) - [nr, nc] = size (x); - x = x - ones (nr, 1) * mean (x); - retval = zeros (1, nc); - s = std (x); - ind = find (s > 0); - retval (ind) = sum (x (:, ind) .^ 4) ./ (nr * s (ind) .^ 4) - 3; - else - error ("kurtosis: x has to be a matrix or a vector."); - endif - -endfunction
--- a/scripts/statistics/mahalanobis.m Mon Oct 19 22:37:57 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,55 +0,0 @@ -## Copyright (C) 1996, 1997 John W. Eaton -## -## 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 2, 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, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: mahalanobis (X, Y) -## -## Returns Mahalanobis' D-square distance between the multivariate -## samples X and Y, which must have the same number of components -## (columns), but may have a different number of observations (rows). - -## Author: Friedrich Leisch <leisch@ci.tuwien.ac.at> -## Created: July 1993 -## Adapted-By: jwe - -function retval = mahalanobis (X, Y) - - if (nargin != 2) - usage ("mahalanobis (X, Y)"); - endif - - [xr, xc] = size (X); - [yr, yc] = size (Y); - - if (xc != yc) - error ("mahalanobis: X and Y must have the same number of columns."); - endif - - Xm = sum (X) / xr; - Ym = sum (Y) / yr; - - X = X - ones (xr, 1) * Xm; - Y = Y - ones (yr, 1) * Ym; - - W = (X' * X + Y' * Y) / (xr + yr - 2); - - Winv = inv (W); - - retval = (Xm - Ym) * Winv * (Xm - Ym)'; - -endfunction
--- a/scripts/statistics/mean.m Mon Oct 19 22:37:57 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,46 +0,0 @@ -## Copyright (C) 1996, 1997 John W. Eaton -## -## 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 2, 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, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: mean (a) -## -## For vector arguments, return the mean the values. -## -## For matrix arguments, return a row vector containing the mean for -## each column. -## -## See also: median, std - -## Author: jwe - -function retval = mean (a) - - if (nargin != 1) - usage ("mean (a)"); - endif - - [nr, nc] = size (a); - if (nr == 1 || nc == 1) - retval = sum (a) / length (a); - elseif (nr > 0 && nc > 0) - retval = sum (a) / nr; - else - error ("mean: invalid matrix argument"); - endif - -endfunction
--- a/scripts/statistics/median.m Mon Oct 19 22:37:57 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,59 +0,0 @@ -## Copyright (C) 1996, 1997 John W. Eaton -## -## 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 2, 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, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: median (a) -## -## For vector arguments, return the median of the values. -## -## For matrix arguments, return a row vector containing the median for -## each column. -## -## See also: std, mean - -## Author: jwe - -function retval = median (a) - - if (nargin != 1) - usage ("median (a)"); - endif - - [nr, nc] = size (a); - s = sort (a); - if (nr == 1 && nc > 0) - if (rem (nc, 2) == 0) - i = nc/2; - retval = (s (i) + s (i+1)) / 2; - else - i = ceil (nc/2); - retval = s (i); - endif - elseif (nr > 0 && nc > 0) - if (rem (nr, 2) == 0) - i = nr/2; - retval = (s (i,:) + s (i+1,:)) / 2; - else - i = ceil (nr/2); - retval = s (i,:); - endif - else - error ("median: invalid matrix argument"); - endif - -endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/models/logistic_regression.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,170 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## Performs ordinal logistic regression. +## +## Suppose Y takes values in k ordered categories, and let gamma_i (x) +## be the cumulative probability that Y falls in one of the first i +## categories given the covariate x. Then +## [theta, beta] = +## logistic_regression (y, x) +## fits the model +## logit (gamma_i (x)) = theta_i - beta' * x, i = 1, ..., k-1. +## The number of ordinal categories, k, is taken to be the number of +## distinct values of round (y) . If k equals 2, y is binary and the +## model is ordinary logistic regression. X is assumed to have full +## column rank. +## +## theta = logistic_regression (y) +## fits the model with baseline logit odds only. +## +## The full form is +## [theta, beta, dev, dl, d2l, gamma] = +## logistic_regression (y, x, print, theta, beta) +## in which all output arguments and all input arguments except y are +## optional. +## +## print = 1 requests summary information about the fitted model to be +## displayed; print = 2 requests information about convergence at each +## iteration. Other values request no information to be displayed. The +## input arguments `theta' and `beta' give initial estimates for theta +## and beta. +## +## `dev' holds minus twice the log-likelihood. +## +## `dl' and `d2l' are the vector of first and the matrix of second +## derivatives of the log-likelihood with respect to theta and beta. +## +## `p' holds estimates for the conditional distribution of Y given x. + +## Original for MATLAB written by Gordon K Smyth <gks@maths.uq.oz.au>, +## U of Queensland, Australia, on Nov 19, 1990. Last revision Aug 3, +## 1992. + +## Author: Gordon K Smyth <gks@maths.uq.oz.au>, +## Adapted-By: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Ordinal logistic regression + +## Uses the auxiliary functions logistic_regression_derivatives and +## logistic_regression_likelihood. + +function [theta, beta, dev, dl, d2l, p] ... + = logistic_regression (y, x, print, theta, beta) + + ## check input + y = round (vec (y)); + [my ny] = size (y); + if (nargin < 2) + x = zeros (my, 0); + endif; + [mx nx] = size (x); + if (mx != my) + error ("x and y must have the same number of observations"); + endif + + ## initial calculations + x = -x; + tol = 1e-6; incr = 10; decr = 2; + ymin = min (y); ymax = max (y); yrange = ymax - ymin; + z = (y * ones (1, yrange)) == ((y * 0 + 1) * (ymin : (ymax - 1))); + z1 = (y * ones (1, yrange)) == ((y * 0 + 1) * ((ymin + 1) : ymax)); + z = z(:, any (z)); + z1 = z1 (:, any(z1)); + [mz nz] = size (z); + + ## starting values + if (nargin < 3) + print = 0; + endif; + if (nargin < 4) + beta = zeros (nx, 1); + endif; + if (nargin < 5) + g = cumsum (sum (z))' ./ my; + theta = log (g ./ (1 - g)); + endif; + tb = [theta; beta]; + + ## likelihood and derivatives at starting values + [g, g1, p, dev] = logistic_regression_likelihood (y, x, tb, z, z1); + [dl, d2l] = logistic_regression_derivatives (x, z, z1, g, g1, p); + epsilon = std (vec (d2l)) / 1000; + + ## maximize likelihood using Levenberg modified Newton's method + iter = 0; + while (abs (dl' * (d2l \ dl) / length (dl)) > tol) + iter = iter + 1; + tbold = tb; + devold = dev; + tb = tbold - d2l \ dl; + [g, g1, p, dev] = logistic_regression_likelihood (y, x, tb, z, z1); + if ((dev - devold) / (dl' * (tb - tbold)) < 0) + epsilon = epsilon / decr; + else + while ((dev - devold) / (dl' * (tb - tbold)) > 0) + epsilon = epsilon * incr; + if (epsilon > 1e+15) + error ("epsilon too large"); + endif + tb = tbold - (d2l - epsilon * eye (size (d2l))) \ dl; + [g, g1, p, dev] = logistic_regression_likelihood (y, x, tb, z, z1); + disp ("epsilon"); disp (epsilon); + endwhile + endif + [dl, d2l] = logistic_regression_derivatives (x, z, z1, g, g1, p); + if (print == 2) + disp ("Iteration"); disp (iter); + disp ("Deviance"); disp (dev); + disp ("First derivative"); disp (dl'); + disp ("Eigenvalues of second derivative"); disp (eig (d2l)'); + endif + endwhile + + ## tidy up output + + theta = tb (1 : nz, 1); + beta = tb ((nz + 1) : (nz + nx), 1); + + if (print >= 1) + printf ("\n"); + printf ("Logistic Regression Results:\n"); + printf ("\n"); + printf ("Number of Iterations: %d\n", iter); + printf ("Deviance: %f\n", dev); + printf ("Parameter Estimates:\n"); + printf (" Theta S.E.\n"); + se = sqrt (diag (inv (-d2l))); + for i = 1 : nz + printf (" %8.4f %8.4f\n", tb (i), se (i)); + endfor + if (nx > 0) + printf (" Beta S.E.\n"); + for i = (nz + 1) : (nz + nx) + printf (" %8.4f %8.4f\n", tb (i), se (i)); + endfor + endif + endif + + if (nargout == 6) + if (nx > 0) + e = ((x * beta) * ones (1, nz)) + ((y * 0 + 1) * theta'); + else + e = (y * 0 + 1) * theta'; + endif + gamma = diff ([(y * 0) exp (e) ./ (1 + exp (e)) (y * 0 + 1)]')'; + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/models/logistic_regression_derivatives.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,37 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## Called by logistic_regression. Calculates derivates of the +## log-likelihood for ordinal logistic regression model. + +## Author: Gordon K. Smyth <gks@maths.uq.oz.au> +## Adapted-By: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Derivates of log-likelihood in logistic regression + +function [dl, d2l] ... + = logistic_regression_derivatives (x, z, z1, g, g1, p) + + ## first derivative + v = g .* (1 - g) ./ p; v1 = g1 .* (1 - g1) ./ p; + dlogp = [dmult (v, z) - dmult (v1, z1) dmult (v - v1, x)]; + dl = sum (dlogp)'; + + ## second derivative + w = v .* (1 - 2 * g); w1 = v1 .* (1 - 2 * g1); + d2l = [z x]' * dmult (w, [z x]) - [z1 x]' * dmult (w1, [z1 x]) ... + - dlogp' * dlogp; + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/models/logistic_regression_likelihood.m Fri Oct 23 05:44:01 1998 +0000 @@ -0,0 +1,34 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## Calculates likelihood for the ordinal logistic regression model. +## Called by logistic_regression. + +## Author: Gordon K. Smyth <gks@maths.uq.oz.au> +## Adapted-By: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Likelihood in logistic regression + +function [g, g1, p, dev] ... + = logistic_regression_likelihood (y, x, beta, z, z1) + + e = exp ([z x] * beta); e1 = exp ([z1 x] * beta); + g = e ./ (1 + e); g1 = e1 ./ (1 + e1); + g = max (y == max (y), g); g1 = min (y > min(y), g1); + + p = g - g1; + dev = -2 * sum (log (p)); + +endfunction
--- a/scripts/statistics/ols.m Mon Oct 19 22:37:57 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,70 +0,0 @@ -## Copyright (C) 1996, 1997 John W. Eaton -## -## 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 2, 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, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: [BETA, SIGMA [, R]] = ols (Y, X) -## -## Ordinary Least Squares (OLS) estimation for the multivariate model -## -## Y = X*B + E, mean(E) = 0, cov(vec(E)) = kron(S,I) -## -## with Y ... T x p As usual, each row of Y and X is an observation -## X ... T x k and each column a variable. -## B ... k x p -## E ... T x p. -## -## BETA is the OLS estimator for B, i.e. -## -## BETA = pinv(X)*Y, -## -## where pinv(X) denotes the pseudoinverse of X. -## SIGMA is the OLS estimator for the matrix S, i.e. -## -## SIGMA = (Y - X*BETA)'*(Y - X*BETA) / (T - rank(X)). -## -## R = Y - X*BETA is the matrix of OLS residuals. - -## Author: Teresa Twaroch <twaroch@ci.tuwien.ac.at> -## Created: May 1993 -## Adapted-By: jwe - -function [BETA, SIGMA, R] = ols (Y, X) - - if (nargin != 2) - error("usage : [BETA, SIGMA [, R]] = ols (Y, X)"); - endif - - [nr, nc] = size (X); - [ry, cy] = size (Y); - if (nr != ry) - error ("ols: incorrect matrix dimensions"); - endif - - Z = X' * X; - r = rank (Z); - - if (r == nc) - BETA = inv (Z) * X' * Y; - else - BETA = pinv (X) * Y; - endif - - R = Y - X * BETA; - SIGMA = R' * R / (nr - r); - -endfunction
--- a/scripts/statistics/skewness.m Mon Oct 19 22:37:57 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,59 +0,0 @@ -## Copyright (C) 1996, 1997 John W. Eaton -## -## 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 2, 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, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: skewness (x) -## -## If x is a vector of length N, return the skewness -## -## skewness (x) = N^(-1) std(x)^(-3) SUM_i (x(i)-mean(x))^3 -## -## of x. -## -## If x is a matrix, return a row vector containing the skewness for each -## column. - -## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> -## Created: 29 July 1994 -## Adapted-By: jwe - -function retval = skewness (x) - - if (nargin != 1) - usage ("skewness (x)"); - endif - - if (is_vector (x)) - x = x - mean (x); - if (! any (x)) - retval = 0; - else - retval = sum (x .^ 3) / (length (x) * std (x) ^ 3); - endif - elseif (is_matrix (x)) - [nr, nc] = size (x); - x = x - ones (nr, 1) * mean (x); - retval = zeros (1, nc); - s = std (x); - ind = find (s > 0); - retval (ind) = sum (x (:, ind) .^ 3) ./ (nr * s (ind) .^ 3); - else - error ("skewness: x has to be a matrix or a vector."); - endif - -endfunction
--- a/scripts/statistics/std.m Mon Oct 19 22:37:57 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,49 +0,0 @@ -## Copyright (C) 1996, 1997 John W. Eaton -## -## 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 2, 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, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: std (a) -## -## For vector arguments, std returns the standard deviation of the -## values. For matrix arguments, std returns a row vector containing -## the standard deviation for each column. -## -## See also: mean, median - -## Author: jwe - -function retval = std (a) - - if (nargin != 1) - usage ("std (a)"); - endif - - nr = rows (a); - nc = columns (a); - if (nc == 1 && nr == 1) - retval = 0; - elseif (nc == 1 || nr == 1) - n = length (a); - retval = sqrt (sumsq (a - mean (a)) / (n - 1)); - elseif (nr > 1 && nc > 0) - retval = sqrt (sumsq (a - ones (nr, 1) * mean (a)) / (nr - 1)); - else - error ("std: invalid matrix argument"); - endif - -endfunction