# For Emacs: -*- mode:cperl; mode:folding; -*-
#
# (c) Jiri Vaclavik
#
package Math::Random::SkewNormal;

# {{{ use

use 5.008;

use strict;
use warnings;

use Perl6::Export::Attrs;
use Readonly;

# }}}
# {{{ variables

my Readonly $pi = 3.14159265358979323846;

our $VERSION = '0.01';

# }}}

# {{{ generate_sn             exported

sub generate_sn : Export {
    my $skewness = shift // 0;
    my $a = generate_n();
    my $b = generate_n();
    my $sconst = $skewness / sqrt(1 + $skewness ** 2);
    my $c = $sconst * $a + sqrt(1 - $sconst ** 2) * $b;
    return $a > 0 ? $c
         : -$c
         ;
}

# }}}
# {{{ generate_sn_multi       exported

sub generate_sn_multi : Export {
    my $delta = shift // return;
    my $omega = shift // return;

    my $n = @$delta;
    return if($n != @$omega);

    # build correlation matrix
    my $cor = [];
    for (0 .. $n-1){
        $cor->[$_][$_] = 1;
    }
    for (1 .. $n){
        $cor->[0][$_] = $delta->[$_-1];
        $cor->[$_][0] = $delta->[$_-1];
    }
    for my $i (1 .. $n){
        for my $j (1 .. $n){
            $cor->[$i][$j] = $omega->[$i-1][$j-1];
        }
    }

    my $sample = generate_n_multi($cor);

    @$sample = map {-$_} @$sample if shift @$sample > 0;

    return $sample;
}

# }}}

# {{{ generate_r              internal

sub generate_r {
    return rand;
}

# }}}
# {{{ generate_n              internal

sub generate_n {
    my $u = generate_r();
    my $v = generate_r();
    return sqrt(-2 * log $u) * cos(2 * $pi * $v);
}

# }}}
# {{{ generate_n_multi        internal

sub generate_n_multi {
    my $cov = shift // return;
    my $n = @$cov || return;
    my $independent = [];
    my $result = [];
    my $sqrt = _choleski($cov) // return;
    for (0 .. $n-1){
        $independent->[$_] = generate_n;
    }
    for my $i (0 .. $n-1){
        $result->[$i] = 0;
        for my $j (0 .. $n-1){
            $result->[$i] += $independent->[$j] * $sqrt->[$i][$j];
        }
    }

    return $result;
}

# }}}
# {{{ _choleski               internal

# choleski's decomposition
sub _choleski {
    my $A = shift // return;
    my $L = [];

    for my $c (0 .. @$A-1) {
        my $sum = 0;
        for my $j (0 .. $c - 1){
            my $i = $c - 1 - $j;
            $sum += $L->[$c][$i] ** 2;
        }
        $L->[$c][$c] = sqrt($A->[$c][$c] - $sum);
        for my $r ($c+1 .. @$A-1) {
            my $sum = 0;
            for my $j (0 .. $c-1){
                my $i = $c - 1 - $j;
                $sum += $L->[$r][$i] * $L->[$c][$i];
            }
            $L->[$r][$c] = ($A->[$r][$c] - $sum) / $L->[$c][$c];
        }
        for my $r (0 .. $c-1) {
            $L->[$r][$c] = 0;
        }
    }

    return $L;
}

# }}}

1;

__END__

# {{{ POD

=head1 NAME

Math::Random::SkewNormal - Handy, easy-to-use Skew-Normal random number generator

=head1 SYNOPSIS

  use Math::Random::SkewNormal;
  
  $skewness = 1;
  $delta    = [0, 0];
  $omega    = [[1, 0], [0, 1]];

  # following functions generate realizations:
  generate_sn($skewness)            # SkewNormal(0,1,$skewness)
  generate_sn_multi($delta, $omega) # SkewNormal_n($delta, $omega)

=head1 DESCRIPTION

This module transforms uniformly spaced random variable realizations into
realizations that follow the Skew-Normal (I<SN>) Probability Density Function (I<PDF>).

We accept following definition of the Skew-Normal Distribution:

=over

=item 1-dimensional SN is determined by PDF: f(x,a) = 2 phi(x) Phi(a x),
    where phi is the PDF and Phi is the CDF of Normal distribution

=item Let X = (X_0, ..., X_k) follows distribution N_k+1 (0, Omega*), where Omega* is
    the matrix of the form:

       [   1   delta ]
       [ delta Omega ]

Omega is a matrix and delta is a vector. Then (X_1, ..., X_n)|(X_0 > 0) follows
n-dimensional Skew Normal distribution.


=back

=head1 ALGORITHM

We use following algorithms:

=over

=item * Box-Muller transformation; for generation Normal realizations

=item * Choleski's decomposition; for inverting matrices

=item * A lemma mentioned e. g. in article from A. Azalini:
    The skew-normal distribution and related multivariate families;
    for generating Skew-Normal realizations

=back

=head1 FUNCTIONS

=head2 generate_sn

Returns realization of SN(0,1, a), where B<a> is the only parameter.
The meaning of B<a> is "skewness" (see definition above).

=head2 generate_sn_multi

Returns realization of centralized n-dimensional Skew-Normal distribution.

1st parameter is correlation matrix of the form C<[[row_1], [row_2], ..., [row_n]]>,
for example C<[[1, 0], [0, 1]]> for 2-dimensional unit matrix.

2nd parameter is delta vector of the form C<[delta_1, ..., delta_n]>, e. g. C<[0, 1]>.

=head1 SEE ALSO

L<Math::Random>, L<Math::Random::Cauchy>,

The examples in F<examples/> subdirectory of this distribution
contains example of generated data with histogram.

The scripts in F<bin/> directory of this distribution contains
tools for generating realizations on command line.

=head1 AUTHOR

Jiri Vaclavik, E<lt>my name dot my last name at gmail dot comE<gt>

=head1 COPYRIGHT AND LICENSE

Copyright (C) 2011 by Jiri Vaclavik

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself, either Perl version 5.8 or,
at your option, any later version of Perl 5 you may have available.

=cut

# }}}