123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687 |
- function alpha = circ_vmrnd(theta, kappa, n)
- % 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
- % By Philipp Berens and Marc J. Velasco, 2009
- % velasco@ccs.fau.edu
- % default parameter
- if nargin < 3
- n = 10;
- end
- if nargin < 2
- kappa = 1;
- end
- if nargin < 1
- theta = 0;
- end
- if numel(n) > 2
- error('n must be a scalar or two-entry vector!')
- elseif numel(n) == 2
- m = n;
- n = n(1) * n(2);
- end
- % if kappa is small, treat as uniform distribution
- if kappa < 1e-6
- alpha = 2*pi*rand(n,1);
- return
- end
- % other cases
- a = 1 + sqrt((1+4*kappa.^2));
- b = (a - sqrt(2*a))/(2*kappa);
- r = (1 + b^2)/(2*b);
- alpha = zeros(n,1);
- for j = 1:n
- while true
- u = rand(3,1);
- z = cos(pi*u(1));
- f = (1+r*z)/(r+z);
- c = kappa*(r-f);
- if u(2) < c * (2-c) || ~(log(c)-log(u(2)) + 1 -c < 0)
- break
- end
-
- end
- alpha(j) = theta + sign(u(3) - 0.5) * acos(f);
- alpha(j) = angle(exp(i*alpha(j)));
- end
- if exist('m','var')
- alpha = reshape(alpha,m(1),m(2));
- end
|