File Coverage

blib/lib/Math/Vector/Real/MultiNormalMixture.pm

Criterion	Covered	Total	%
statement	15	71	21.1
branch	0	10	0.0
condition	0	8	0.0
subroutine	5	10	50.0
pod	5	5	100.0
total	25	104	24.0

package Math::Vector::Real::MultiNormalMixture;

our $VERSION = '0.02';

use 5.010;
use strict;
use warnings;
use Carp;

use Math::Vector::Real;

my $PI = 3.14159265358979323846264338327950288419716939937510;

sub new {
    my ($class, %opts) = @_;
    my $mu = delete $opts{mu} // croak "required argument 'mu' missing";
    @$mu >= 1 or die "'mu' must containt one point at least";
    my $dim = @{$mu->[0]};
    my $n = @$mu;
    my $alpha = delete $opts{alpha};
    my @alpha = (defined $alpha ? @$alpha : ((1/$n) x $n));
    my $sigma = delete $opts{sigma} // 1.0;
    my @sigma = (ref $sigma ? @$sigma : (($sigma) x $n));
    croak "bad array size" unless (@alpha == $n and @sigma == $n);
    my $c1 = (2 * $PI) ** (-0.5 * $dim);
    my @c = map $alpha[$_] * $c1 * ($sigma[$_] ** (-$dim)), (0..$#alpha);
    my @isigma2 = map 1.0/($_ * $_), @sigma;
    my $c2 = -sqrt(2.0 / $PI);
    my @second = map $c2/($_ * $_ * $_), @sigma;

    my $self = { mu => [map Math::Vector::Real::clone($_), @$mu],
                 alpha => \@alpha,
                 sigma => \@sigma,
                 isigma2 => \@isigma2,
                 c => \@c,
                 second => \@second,
                 dim => $dim,
               };
    bless $self, $class;
}

sub density {
    my ($self, $p) = @_;
    my $c = $self->{c};
    my $mu = $self->{mu};
    my $isigma2 = $self->{sigma};
    my $acu = 0;
    for (0..$#$mu) {
        $acu += $c->[$_] * exp(-$isigma2->[$_] * $mu->[$_]->dist2($p));
    }
    $acu;
}

sub density_portion {
    my $self = shift;
    my $p = shift;
        my $c = $self->{c};
    my $mu = $self->{mu};
    my $isigma2 = $self->{sigma};
    my $acu = 0;
    for (@_) {
        $acu += $c->[$_] * exp(-$isigma2->[$_] * $mu->[$_]->dist2($p));
    }
    $acu;
}

sub density_and_gradient {
    my ($self, $p) = @_;
    my $c = $self->{c};
    my $mu = $self->{mu};
    my $isigma2 = $self->{sigma};
    my $d = 0;
    my $g = $Math::Vector::Real->zero(scalar @{$mu->[0]});
    for (0..$#$mu) {
        my $mu_p = $p - $mu->[$_];
        my $isigma2 = $isigma2->[$_];
        my $Dd = $c->[$_] * exp(-$isigma2 * $mu_p->abs2);
        $d += $Dd;
        $g += -2 * $isigma2 * $mu_p * $Dd;
    }
    return ($d, $g);
}

sub max_density_estimation {
    my $self = shift;
    my $mu = $self->{mu};
    my $max = $self->density($mu->[0]);
    for my $ix (1..$#$mu) {
        my $d = $self->density($mu->[$ix]);
        # print "d: $d\n";
        $max = $d if $d > $max;
    }
    return $max;
}

1;
__END__