From bug-request at octave dot org Wed Nov 23 10:41:37 2005 Subject: lognormal distribution defs From: Ben Barrowes To: bug at octave dot org Date: Wed, 23 Nov 2005 11:38:39 -0500 The logn***.m functions in octave are different from their counterparts in ML. Both ML: http://www.mathworks.com/access/helpdesk/help/toolbox/stats/lognpdf.html and mathworld: http://mathworld.wolfram.com/LogNormalDistribution.html define this dist not with log(mean), but with the mean itself. For example on lognpdf.m ML: >> x=1:.25:3 x = Columns 1 through 4 1 1.25 1.5 1.75 Columns 5 through 8 2 2.25 2.5 2.75 Column 9 3 >> y=lognpdf(x,1.2,.2) y = Columns 1 through 4 3.03794142491165e-08 1.05355986355732e-05 0.000497439354500118 0.00676992245065161 Columns 5 through 8 0.0401997276397916 0.133637810505383 0.29172951838526 0.465436670152427 Column 9 0.58472930838224 Octave: octave:42> x=1:.25:3 x = Columns 1 through 4: 1.00000000000000 1.25000000000000 1.50000000000000 1.75000000000000 Columns 5 through 8: 2.00000000000000 2.25000000000000 2.50000000000000 2.75000000000000 Column 9: 3.00000000000000 octave:43> y=lognpdf(x,1.2,.2) y = Columns 1 through 3: 1.31651093144260e+00 1.56287236440818e+00 7.13638473120990e-01 Columns 4 through 6: 1.92337599676002e-01 3.82185979303903e-02 6.34688320301560e-03 Columns 7 through 9: 9.49404390608672e-04 1.34035664837834e-04 1.84034626064715e-05 octave:44> y=lognpdf(x,exp(1.2),.2) y = Columns 1 through 3: 3.03794142491166e-08 1.05355986355732e-05 4.97439354500118e-04 Columns 4 through 6: 6.76992245065160e-03 4.01997276397916e-02 1.33637810505383e-01 Columns 7 through 9: 2.91729518385260e-01 4.65436670152427e-01 5.84729308382240e-01 octave:45> Here is a patch for lognpdf.m to bring it into conformity with ML: --- lognpdf.m.old 2005-11-23 10:12:21.000000000 -0500 +++ lognpdf.m 2005-11-23 10:42:37.000000000 -0500 at @ -18,56 +18,56 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## at deftypefn {Function File} {} lognpdf (@var{x}, @var{a}, @var{v}) +## at deftypefn {Function File} {} lognpdf (@var{x}, @var{mu}, @var{sigma}) ## For each element of at var{x}, compute the probability density function ## (PDF) at at var{x} of the lognormal distribution with parameters -## at var{a} and @var{v}. If a random variable follows this distribution, -## its logarithm is normally distributed with mean at code{log (@var{a})} -## and variance at var{v} dot +## at var{mu} and @var{sigma}. If a random variable follows this distribution, +## its logarithm is normally distributed with mean at var{mu} +## and standard deviation at var{sigma} dot ## -## Default values are at var{a} = 1, @var{v} = 1. +## Default values are at var{mu} = 1, @var{sigma} = 1. ## at end deftypefn ## Author: KH ## Description: PDF of the log normal distribution -function pdf = lognpdf (x, a, v) +function pdf = lognpdf (x, mu, sigma) if (! ((nargin == 1) || (nargin == 3))) - usage ("lognpdf (x, a, v)"); + usage ("lognpdf (x, mu, sigma)"); endif if (nargin == 1) - a = 1; - v = 1; + mu = 1; + sigma = 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); + ## pdf = (x > 0) ./ x .* normal_pdf (log (x), mu, sigma); ## Hence ... - if (!isscalar (a) || !isscalar (v)) - [retval, x, a, v] = common_size (x, a, v); + if (!isscalar (mu) || !isscalar (sigma)) + [retval, x, mu, sigma] = common_size (x, mu, sigma); if (retval > 0) - error ("lognpdf: x, a and v must be of common size or scalars"); + error ("lognpdf: x, mu and sigma must be of common size or scalars"); endif endif pdf = zeros (size (x)); - k = find (isnan (x) | !(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); + k = find (isnan (x) | !(mu > 0) | !(mu < Inf) | !(sigma > 0) | !(sigma < Inf)); if (any (k)) pdf(k) = NaN; endif - k = find ((x > 0) & (x < Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + k = find ((x > 0) & (x < Inf) & (mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf)); if (any (k)) - if (isscalar (a) && isscalar (v)) - pdf(k) = normpdf (log (x(k)), log (a), v) ./ x(k); + if (isscalar (mu) && isscalar (sigma)) + pdf(k) = normpdf (log (x(k)), mu, sigma) ./ x(k); else - pdf(k) = normpdf (log (x(k)), log (a(k)), v(k)) ./ x(k); - endif + pdf(k) = normpdf (log (x(k)), mu(k), sigma(k)) ./ x(k); + endif endif endfunction logncdf.m has the same problem. It also has the problem that it expects the variance rather than the standard deviation, so here is a patch to update it: --- logncdf.m.old 2005-11-23 10:26:43.000000000 -0500 +++ logncdf.m 2005-11-23 10:36:22.000000000 -0500 at @ -18,60 +18,60 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## at deftypefn {Function File} {} logncdf (@var{x}, @var{a}, @var{v}) +## at deftypefn {Function File} {} logncdf (@var{x}, @var{mu}, @var{sigma}) ## For each element of at var{x}, compute the cumulative distribution ## function (CDF) at at var{x} of the lognormal distribution with -## parameters at var{a} and @var{v}. If a random variable follows this +## parameters at var{mu} and @var{sigma}. If a random variable follows this ## distribution, its logarithm is normally distributed with mean -## at code{log (@var{a})} and variance @var{v}. +## at var{mu} and standard deviation @var{sigma}. ## -## Default values are at var{a} = 1, @var{v} = 1. +## Default values are at var{mu} = 1, @var{sigma} = 1. ## at end deftypefn ## Author: KH ## Description: CDF of the log normal distribution -function cdf = logncdf (x, a, v) +function cdf = logncdf (x, mu, sigma) if (! ((nargin == 1) || (nargin == 3))) - usage ("logncdf (x, a, v)"); + usage ("logncdf (x, mu, sigma)"); endif if (nargin == 1) - a = 1; - v = 1; + mu = 1; + sigma = 1; endif ## The following "straightforward" implementation unfortunately does ## not work (because exp (Inf) -> NaN etc): - ## cdf = normal_cdf (log (x), log (a), v); + ## cdf = normal_cdf (log (x), log (mu), sigma); ## Hence ... - if (!isscalar (a) || !isscalar (v)) - [retval, x, a, v] = common_size (x, a, v); + if (!isscalar (mu) || !isscalar (sigma)) + [retval, x, mu, sigma] = common_size (x, mu, sigma); if (retval > 0) - error ("logncdf: x, a and v must be of common size or scalars"); + error ("logncdf: x, mu and sigma must be of common size or scalars"); endif endif cdf = zeros (size (x)); - k = find (isnan (x) | !(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); + k = find (isnan (x) | !(mu > 0) | !(mu < Inf) | !(sigma > 0) | !(sigma < Inf)); if (any (k)) cdf(k) = NaN; endif - k = find ((x == Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + k = find ((x == Inf) & (mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf)); if (any (k)) cdf(k) = 1; endif - k = find ((x > 0) & (x < Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + k = find ((x > 0) & (x < Inf) & (mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf)); if (any (k)) - if (isscalar (a) && isscalar (v)) - cdf(k) = stdnormal_cdf ((log (x(k)) - log (a)) / sqrt (v)); + if (isscalar (mu) && isscalar (sigma)) + cdf(k) = stdnormal_cdf ((log (x(k)) - mu) / sigma); else - cdf(k) = stdnormal_cdf ((log (x(k)) - log (a(k))) ./ sqrt (v(k))); + cdf(k) = stdnormal_cdf ((log (x(k)) - mu(k)) ./ sigma(k)); endif endif logninv.m has the same issues: ML: >> y=logninv(.2:.1:.8,1.2,.2) y = Columns 1 through 4 2.80576366507724 2.98954110028598 3.15607945768458 3.32011692273655 Columns 5 through 7 3.49268024726115 3.68724697566029 3.92876154105379 Octave: octave:73> y=logninv(.2:.1:.8,1.2,.2) y = Columns 1 through 4: 0.823606096892220 0.949141410012455 1.071459018008457 1.200000000000000 Columns 5 through 7: 1.343961808895461 1.517160651520939 1.748408620860954 octave:74> y=logninv(.2:.1:.8,exp(1.2),.2^2) y = Columns 1 through 4: 2.80576366507724 2.98954110028598 3.15607945768458 3.32011692273655 Columns 5 through 7: 3.49268024726115 3.68724697566029 3.92876154105379 octave:75> --- logninv.m.old 2005-11-23 11:20:29.000000000 -0500 +++ logninv.m 2005-11-23 11:27:34.000000000 -0500 at @ -18,61 +18,61 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## at deftypefn {Function File} {} logninv (@var{x}, @var{a}, @var{v}) +## at deftypefn {Function File} {} logninv (@var{x}, @var{mu}, @var{sigma}) ## For each element of at var{x}, compute the quantile (the inverse of the -## CDF) at at var{x} of the lognormal distribution with parameters @var{a} -## and at var{v} dot If a random variable follows this distribution, its -## logarithm is normally distributed with mean at code{log (@var{a})} and -## variance at var{v} dot +## CDF) at at var{x} of the lognormal distribution with parameters @var{mu} +## and at var{sigma} dot If a random variable follows this distribution, its +## logarithm is normally distributed with mean at var{mu} and +## standard deviation at var{sigma} dot ## -## Default values are at var{a} = 1, @var{v} = 1. +## Default values are at var{mu} = 1, @var{sigma} = 1. ## at end deftypefn ## Author: KH ## Description: Quantile function of the log normal distribution -function inv = logninv (x, a, v) +function inv = logninv (x, mu, sigma) if (! ((nargin == 1) || (nargin == 3))) - usage ("logninv (x, a, v)"); + usage ("logninv (x, mu, sigma)"); endif if (nargin == 1) - a = 1; - v = 1; + mu = 1; + sigma = 1; endif ## The following "straightforward" implementation unfortunately does ## not work (because exp (Inf) -> NaN): - ## inv = exp (normal_inv (x, log (a), v)); + ## inv = exp (normal_inv (x, mu, sigma)); ## Hence ... - if (!isscalar (a) || !isscalar (v)) - [retval, x, a, v] = common_size (x, a, v); + if (!isscalar (mu) || !isscalar (sigma)) + [retval, x, mu, sigma] = common_size (x, mu, sigma); if (retval > 0) - error ("logninv: x, a and v must be of common size or scalars"); + error ("logninv: x, mu and sigma must be of common size or scalars"); endif endif inv = zeros (size (x)); - k = find (!(x >= 0) | !(x <= 1) | !(a > 0) | !(a < Inf) - | !(v > 0) | !(v < Inf)); + k = find (!(x >= 0) | !(x <= 1) | !(mu > 0) | !(mu < Inf) + | !(sigma > 0) | !(sigma < Inf)); if (any (k)) inv(k) = NaN; endif - k = find ((x == 1) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + k = find ((x == 1) & (mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf)); if (any (k)) inv(k) = Inf; endif - k = find ((x > 0) & (x < 1) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + k = find ((x > 0) & (x < 1) & (mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf)); if (any (k)) - if (isscalar (a) && isscalar (v)) - inv(k) = a .* exp (sqrt (v) .* stdnormal_inv (x(k))); + if (isscalar (mu) && isscalar (sigma)) + inv(k) = exp(mu) .* exp (sigma .* stdnormal_inv (x(k))); else - inv(k) = a(k) .* exp (sqrt (v(k)) .* stdnormal_inv (x(k))); + inv(k) = exp(mu(k)) .* exp (sigma(k) .* stdnormal_inv (x(k))); endif endif and finally lognrnd.m --- lognrnd.m.old 2005-11-23 11:36:23.000000000 -0500 +++ lognrnd.m 2005-11-23 11:36:37.000000000 -0500 at @ -1,3 +1,4 @@ + ## Copyright (C) 1995, 1996, 1997 Kurt Hornik ## ## This file is part of Octave. at @ -18,27 +19,27 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## at deftypefn {Function File} {} lognrnd (@var{a}, @var{v}, @var{r}, @var{c}) -## at deftypefnx {Function File} {} lognrnd (@var{a}, @var{v}, @var{sz}) +## at deftypefn {Function File} {} lognrnd (@var{mu}, @var{sigma}, @var{r}, @var{c}) +## at deftypefnx {Function File} {} lognrnd (@var{a}, @var{sigma}, @var{sz}) ## Return an at var{r} by @var{c} matrix of random samples from the -## lognormal distribution with parameters at var{a} and @var{v}. Both -## at var{a} and @var{v} must be scalar or of size @var{r} by @var{c}. +## lognormal distribution with parameters at var{mu} and @var{sigma}. Both +## at var{mu} and @var{sigma} must be scalar or of size @var{r} by @var{c}. ## Or if at var{sz} is a vector, create a matrix of size @var{sz}. ## ## If at var{r} and @var{c} are omitted, the size of the result matrix is -## the common size of at var{a} and @var{v}. +## the common size of at var{mu} and @var{sigma}. ## at end deftypefn ## Author: KH ## Description: Random deviates from the log normal distribution -function rnd = lognrnd (a, v, r, c) +function rnd = lognrnd (mu, sigma, r, c) if (nargin > 1) - if (!isscalar(a) || !isscalar(v)) - [retval, a, v] = common_size (a, v); + if (!isscalar(mu) || !isscalar(sigma)) + [retval, mu, sigma] = common_size (mu, sigma); if (retval > 0) - error ("lognrnd: a and v must be of common size or scalar"); + error ("lognrnd: mu and sigma must be of common size or scalar"); endif endif endif at @ -52,9 +53,9 @@ endif sz = [r, c]; - if (any (size (a) != 1) && - ((length (size (a)) != length (sz)) || any (size (a) != sz))) - error ("lognrnd: a and b must be scalar or of size [r, c]"); + if (any (size (mu) != 1) && + ((length (size (mu)) != length (sz)) || any (size (mu) != sz))) + error ("lognrnd: mu and sigma must be scalar or of size [r, c]"); endif elseif (nargin == 3) at @ -66,35 +67,36 @@ error ("lognrnd: r must be a postive integer or vector"); endif - if (any (size (a) != 1) && - ((length (size (a)) != length (sz)) || any (size (a) != sz))) - error ("lognrnd: a and b must be scalar or of size sz"); + if (any (size (mu) != 1) && + ((length (size (mu)) != length (sz)) || any (size (mu) != sz))) + error ("lognrnd: mu and sigma must be scalar or of size sz"); endif elseif (nargin == 2) - sz = size(a); + sz = size(mu); else - usage ("lognrnd (a, v, r, c)"); + usage ("lognrnd (mu, sigma, r, c)"); endif - if (isscalar (a) && isscalar (v)) - if (!(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)) + if (isscalar (mu) && isscalar (sigma)) + if (!(mu > 0) | !(mu < Inf) | !(sigma > 0) | !(sigma < Inf)) rnd = NaN * ones (sz); - elseif find ((a > 0) & (a < Inf) & (v > 0) & (v < Inf)); - rnd = a * exp (sqrt (v) .* randn (sz)); + elseif find ((mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf)); + rnd = exp (mu) * exp (sigma .* randn (sz)); else rnd = zeros (sz); endif else rnd = zeros (sz); - k = find (!(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); + k = find (!(mu > 0) | !(mu < Inf) | !(sigma > 0) | !(sigma < Inf)); if (any (k)) rnd(k) = NaN * ones (1, length (k)); endif - k = find ((a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + k = find ((mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf)); if (any (k)) - rnd(k) = a(k) .* exp (sqrt (v(k)) .* randn (1, length (k))); + rnd(k) = exp (mu(k)) .* exp (sigma(k) .* randn (1, length (k))); endif endif endfunction + bb ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------