sg_random_stable

SG_RANDOM_STABLE  create MM x NN matrix of random numbers from stable dist
   S_PHI = SG_RANDOM_STABLE(A_STAB, B_STAB, G_STAB, D_STAB, MM, NN)

   Input:
       a_stab - alpha, exponent
       b_stab - beta, skewness
       g_stab - gama, scale parameter
       d_stab - delta, location parameter
       mm,nn  - a mm by nn matrix

   Output:
       s_phi - matrix (mm by nn) of random numbers based on stable
       distribution

0001 function s_phi = sg_random_stable (a_stab, b_stab, g_stab, d_stab, mm, nn)
0002 %SG_RANDOM_STABLE  create MM x NN matrix of random numbers from stable dist
0003 %   S_PHI = SG_RANDOM_STABLE(A_STAB, B_STAB, G_STAB, D_STAB, MM, NN)
0004 %
0005 %   Input:
0006 %       a_stab - alpha, exponent
0007 %       b_stab - beta, skewness
0008 %       g_stab - gama, scale parameter
0009 %       d_stab - delta, location parameter
0010 %       mm,nn  - a mm by nn matrix
0011 %
0012 %   Output:
0013 %       s_phi - matrix (mm by nn) of random numbers based on stable
0014 %       distribution
0015 
0016 %   SynGrid
0017 %   Copyright (c) 2018, Electric Power and Energy Systems (EPES) Research Lab
0018 %   by Hamidreza Sadeghian and Zhifang Wang, Virginia Commonwealth University
0019 %
0020 %   This file is part of SynGrid.
0021 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0022 
0023 W = - log( rand(mm,nn) );
0024 V = pi/2 * (2*rand(mm,nn) - 1);
0025 Cnst = b_stab * tan(pi*a_stab/2);
0026 B = atan( Cnst );
0027 S = (1 + Cnst * Cnst).^(1/(2*a_stab));
0028 Rand = S * sin( a_stab*V + B ) ./ ( cos(V) ).^(1/a_stab) .*( cos( (1-a_stab) * V - B ) ./ W ).^((1-a_stab)/a_stab);
0029 s_phi= g_stab * Rand + d_stab;