*/ protected int $k; /** * The amount of epsilon smoothing added to the variance of each feature. * * @var float */ protected float $smoothing; /** * The maximum number of iterations to run until the algorithm terminates. * * @var int */ protected int $epochs; /** * The minimum shift in the components necessary to continue training. * * @var float */ protected float $minChange; /** * The cluster centroid seeder. * * @var Seeder */ protected Seeder $seeder; /** * The precomputed log prior probabilities of each cluster. * * @var float[] */ protected array $logPriors = [ // ]; /** * The computed means of each feature column for each gaussian. * * @var list> */ protected array $means = [ // ]; /** * The computed variances of each feature column for each gaussian. * * @var list> */ protected array $variances = [ // ]; /** * The loss at each epoch from the last training session. * * @var float[]|null */ protected ?array $losses = null; /** * @param int $k * @param float $smoothing * @param int $epochs * @param float $minChange * @param Seeder|null $seeder * @throws InvalidArgumentException */ public function __construct( int $k, float $smoothing = 1e-9, int $epochs = 100, float $minChange = 1e-3, ?Seeder $seeder = null ) { if ($k < 1) { throw new InvalidArgumentException('K must be greater' . " than 0, $k given."); } if ($smoothing <= 0.0) { throw new InvalidArgumentException('Smoothing must be' . " greater than 0, $smoothing given."); } if ($epochs < 0) { throw new InvalidArgumentException('Number of epochs' . " must be greater than 0, $epochs given."); } if ($minChange < 0.0) { throw new InvalidArgumentException('Minimum change must be' . " greater than 0, $minChange given."); } $this->k = $k; $this->smoothing = $smoothing; $this->epochs = $epochs; $this->minChange = $minChange; $this->seeder = $seeder ?? new PlusPlus(); } /** * Return the estimator type. * * @return EstimatorType */ public function type() : EstimatorType { return EstimatorType::clusterer(); } /** * Return the data types that the estimator is compatible with. * * @return list */ public function compatibility() : array { return [ DataType::continuous(), ]; } /** * Return the settings of the hyper-parameters in an associative array. * * @return mixed[] */ public function params() : array { return [ 'k' => $this->k, 'smoothing' => $this->smoothing, 'epochs' => $this->epochs, 'min change' => $this->minChange, 'seeder' => $this->seeder, ]; } /** * Has the learner been trained? * * @return bool */ public function trained() : bool { return $this->means and $this->variances; } /** * Return the cluster prior probabilities. * * @return float[] */ public function priors() : array { return array_map('exp', $this->logPriors); } /** * Return the mean vectors of each component. * * @return list> */ public function means() : array { return $this->means; } /** * Return the multivariate variance of each component. * * @return list> */ public function variances() : array { return $this->variances; } /** * Return an iterable progress table with the steps from the last training session. * * @return Generator */ public function steps() : Generator { if (!$this->losses) { return; } foreach ($this->losses as $epoch => $loss) { yield [ 'epoch' => $epoch, 'loss' => $loss, ]; } } /** * Return the loss at each epoch of training from the last training session. * * @return float[]|null */ public function losses() : ?array { return $this->losses; } /** * Train the learner with a dataset. * * @param Dataset $dataset */ public function train(Dataset $dataset) : void { SpecificationChain::with([ new DatasetIsNotEmpty($dataset), new SamplesAreCompatibleWithEstimator($dataset, $this), ])->check(); if ($this->logger) { $this->logger->info("Training $this"); } $this->initialize($dataset); $samples = $dataset->samples(); $features = $dataset->features(); $n = $dataset->numSamples(); $minEpsilon = CPU::epsilon(); $prevLoss = INF; $this->losses = []; for ($epoch = 1; $epoch <= $this->epochs; ++$epoch) { $loss = $maxVariance = 0.0; $memberships = []; foreach ($samples as $sample) { $jll = $this->jointLogLikelihood($sample); $total = logsumexp($jll); $loss -= $total; $dist = []; foreach ($jll as $cluster => $likelihood) { $dist[$cluster] = exp($likelihood - $total); } $memberships[] = $dist; } if (is_nan($loss)) { if ($this->logger) { $this->logger->warning('Numerical instability detected'); } break; } for ($cluster = 0; $cluster < $this->k; ++$cluster) { $affinities = array_column($memberships, $cluster); $total = array_sum($affinities); $means = $variances = []; foreach ($features as $column) { $sigma = $ssd = 0.0; foreach ($column as $i => $value) { $sigma += $affinities[$i] * $value; } $mean = $sigma / $total; foreach ($column as $i => $value) { $ssd += $affinities[$i] * ($value - $mean) ** 2; } $variance = $ssd / $total; $means[] = $mean; $variances[] = $variance; } $maxVariance = max($maxVariance, ...$variances); $logPrior = log($total / $n); $this->means[$cluster] = $means; $this->variances[$cluster] = $variances; $this->logPriors[$cluster] = $logPrior; } $epsilon = max($this->smoothing * $maxVariance, $minEpsilon); foreach ($this->variances as &$variances) { foreach ($variances as &$variance) { $variance += $epsilon; } } $loss /= $n; $lossChange = abs($loss - $prevLoss); $this->losses[$epoch] = $loss; if ($this->logger) { $lossDirection = $loss < $prevLoss ? '↓' : '↑'; $message = "Epoch: $epoch, " . "Loss: $loss, " . "Loss Change: {$lossDirection}{$lossChange}"; $this->logger->info($message); } if ($loss <= 0.0) { break; } if ($lossChange < $this->minChange) { break; } $prevLoss = $loss; } if ($this->logger) { $this->logger->info('Training complete'); } } /** * Make predictions from a dataset. * * @param Dataset $dataset * @throws RuntimeException * @return list */ public function predict(Dataset $dataset) : array { if (empty($this->means) or empty($this->variances)) { throw new RuntimeException('Estimator has not been trained.'); } DatasetHasDimensionality::with($dataset, count(current($this->means) ?: []))->check(); return array_map([$this, 'predictSample'], $dataset->samples()); } /** * Predict a single sample and return the result. * * @internal * * @param (int|float)[] $sample * @return int */ public function predictSample(array $sample) : int { return argmax($this->jointLogLikelihood($sample)); } /** * Estimate the joint probabilities for each possible outcome. * * @param Dataset $dataset * @throws RuntimeException * @return list */ public function proba(Dataset $dataset) : array { if (empty($this->means) or empty($this->variances)) { throw new RuntimeException('Estimator has not been trained.'); } DatasetHasDimensionality::with($dataset, count(current($this->means) ?: []))->check(); return array_map([$this, 'probaSample'], $dataset->samples()); } /** * Return the membership of a sample to each of the c centroids. * * @param list $sample * @return float[] */ protected function probaSample(array $sample) : array { $jll = $this->jointLogLikelihood($sample); $total = logsumexp($jll); $dist = []; foreach ($jll as $cluster => $likelihood) { $dist[$cluster] = exp($likelihood - $total); } return $dist; } /** * Calculate the joint log likelihood of a sample being a member * of each of the gaussian components. * * @param list $sample * @return array */ protected function jointLogLikelihood(array $sample) : array { $likelihoods = []; foreach ($this->logPriors as $cluster => $prior) { $means = $this->means[$cluster]; $variances = $this->variances[$cluster]; $likelihood = $prior; foreach ($sample as $column => $feature) { $mean = $means[$column]; $variance = $variances[$column]; $pdf = -0.5 * log(TWO_PI * $variance); $pdf -= 0.5 * (($feature - $mean) ** 2) / $variance; $likelihood += $pdf; } $likelihoods[$cluster] = $likelihood; } return $likelihoods; } /** * Initialize the gaussian components by calculating the means and variances of k initial cluster centroids generated by the seeder. * * @param Dataset $dataset */ protected function initialize(Dataset $dataset) : void { $this->logPriors = $this->means = $this->variances = []; $kernel = new Euclidean(); $n = $dataset->numSamples(); $maxVariance = 0.0; /** @var list> $centroids */ $centroids = $this->seeder->seed($dataset, $this->k); $clusters = array_fill(0, $this->k, []); foreach ($dataset->samples() as $sample) { $bestDistance = INF; $bestCluster = -1; foreach ($centroids as $cluster => $centroid) { $distance = $kernel->compute($sample, $centroid); if ($distance < $bestDistance) { $bestDistance = $distance; $bestCluster = $cluster; } } $clusters[$bestCluster][] = $sample; } foreach ($clusters as $cluster => $samples) { $means = $variances = []; $features = array_transpose($samples); foreach ($features as $values) { [$mean, $variance] = Stats::meanVar($values); $means[] = $mean; $variances[] = $variance; } $maxVariance = max($maxVariance, ...$variances); $logPrior = log(count($samples) / $n); $this->means[$cluster] = $means; $this->variances[$cluster] = $variances; $this->logPriors[$cluster] = $logPrior; } $epsilon = max($this->smoothing * $maxVariance, CPU::epsilon()); foreach ($this->variances as &$variances) { foreach ($variances as &$variance) { $variance += $epsilon; } } } /** * Return an associative array containing the data used to serialize the object. * * @return mixed[] */ public function __serialize() : array { $properties = get_object_vars($this); unset($properties['losses'], $properties['logger']); return $properties; } /** * Return the string representation of the object. * * @internal * * @return string */ public function __toString() : string { return 'Gaussian Mixture (' . Params::stringify($this->params()) . ')'; } }