File Coverage

blib/lib/Math/Random/SkewNormal.pm

Criterion	Covered	Total	%
statement	10	12	83.3
branch			n/a
condition			n/a
subroutine	4	4	100.0
pod			n/a
total	14	16	87.5

# 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__