Home > atmlab > circular > circ_vmrnd.m

circ_vmrnd

PURPOSE ^

alpha = circ_vmrnd(theta, kappa, n)

SYNOPSIS ^

function alpha = circ_vmrnd(theta, kappa, n)

DESCRIPTION ^

 alpha = circ_vmrnd(theta, kappa, n)
   Simulates n random angles from a von Mises distribution, with preferred 
   direction thetahat and concentration parameter kappa.

   Input:
     [theta    preferred direction, default is 0]
     [kappa    width, default is 1]
     [n        number of samples, default is 10]

     If n is a vector with two entries (e.g. [2 10]), the function creates
     a matrix output with the respective dimensionality.

   Output:
     alpha     samples from von Mises distribution


   References:
     Statistical analysis of circular data, Fisher, sec. 3.3.6, p. 49

 Circular Statistics Toolbox for Matlab

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

circ_vmrnd.m

SOURCE CODE ^

0001 function alpha = circ_vmrnd(theta, kappa, n)
0002 
0003 % alpha = circ_vmrnd(theta, kappa, n)
0004 %   Simulates n random angles from a von Mises distribution, with preferred
0005 %   direction thetahat and concentration parameter kappa.
0006 %
0007 %   Input:
0008 %     [theta    preferred direction, default is 0]
0009 %     [kappa    width, default is 1]
0010 %     [n        number of samples, default is 10]
0011 %
0012 %     If n is a vector with two entries (e.g. [2 10]), the function creates
0013 %     a matrix output with the respective dimensionality.
0014 %
0015 %   Output:
0016 %     alpha     samples from von Mises distribution
0017 %
0018 %
0019 %   References:
0020 %     Statistical analysis of circular data, Fisher, sec. 3.3.6, p. 49
0021 %
0022 % Circular Statistics Toolbox for Matlab
0023 
0024 % By Philipp Berens and Marc J. Velasco, 2009
0025 % velasco@ccs.fau.edu
0026 
0027 
0028 % default parameter
0029 if nargin < 3
0030     n = 10;
0031 end
0032 
0033 if nargin < 2
0034     kappa = 1;
0035 end
0036 
0037 if nargin < 1
0038     theta = 0;
0039 end
0040 
0041 if numel(n) > 2
0042   error('n must be a scalar or two-entry vector!')
0043 elseif numel(n) == 2
0044   m = n;
0045   n = n(1) * n(2);
0046 end  
0047 
0048 % if kappa is small, treat as uniform distribution
0049 if kappa < 1e-6
0050     alpha = 2*pi*rand(n,1);
0051     return
0052 end
0053 
0054 % other cases
0055 a = 1 + sqrt((1+4*kappa.^2));
0056 b = (a - sqrt(2*a))/(2*kappa);
0057 r = (1 + b^2)/(2*b);
0058 
0059 alpha = zeros(n,1);
0060 for j = 1:n
0061   while true
0062       u = rand(3,1);
0063 
0064       z = cos(pi*u(1));
0065       f = (1+r*z)/(r+z);
0066       c = kappa*(r-f);
0067 
0068       if u(2) < c * (2-c) || ~(log(c)-log(u(2)) + 1 -c < 0)
0069          break
0070       end
0071 
0072       
0073   end
0074 
0075   alpha(j) = theta +  sign(u(3) - 0.5) * acos(f);
0076   alpha(j) = angle(exp(i*alpha(j)));
0077 end
0078 
0079 if exist('m','var')
0080   alpha = reshape(alpha,m(1),m(2));
0081 end
0082 
0083 
0084 
0085 
0086 
0087

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005