dochia-cycle-model / src /MicroModel.php

Dochia v0.1.0: model card, params and reference implementations

e215d3c 9 days ago

14.9 kB

<?php

declare(strict_types=1);

namespace Dochia;

/**
 * Dochia v1 — reference implementation (pure PHP, no deps).
 *
 * Pure PHP reference implementation; mirrors python/dochia.py. The golden
 * fixtures in tests/fixtures/golden.json bind the two implementations.
 */
final class MicroModel
{
    public function __construct(public readonly ModelParams $params)
    {
    }

    /**
     * Personalise a user's cycle-length distribution from observed lengths.
     *
     * @param  list<int|float>  $lengths
     * @return array<string, mixed>
     */
    public function personalState(array $lengths): array
    {
        $lengths = $this->filterLengths($lengths);
        $n = count($lengths);

        if ($n === 0) {
            $withinVar = $this->params->withinSdPrior ** 2;
            $betweenVar = $this->params->popBetweenSd ** 2;
            $posteriorVar = ($withinVar * $betweenVar) / max(1e-9, $withinVar);
            return [
                'mean' => $this->params->popMean,
                'within_var' => $withinVar,
                'cycle_sd' => sqrt($withinVar + $posteriorVar),
                'n' => 0,
                'tail_prob' => $this->params->tailProbBase,
                'trend' => 0.0,
                'recent_sd' => 0.0,
                'recent_range' => 0,
                'gap' => false,
                'abstain' => false,
                'confidence' => 'new',
            ];
        }

        $mean = array_sum($lengths) / $n;
        $sampleVar = 0.0;
        if ($n > 1) {
            $sumSq = 0.0;
            foreach ($lengths as $y) {
                $sumSq += ($y - $mean) ** 2;
            }
            $sampleVar = $sumSq / ($n - 1);
        }
        $median = self::median($lengths);
        $absDev = [];
        foreach ($lengths as $y) {
            $absDev[] = abs($y - $median);
        }
        $robustSd = 1.4826 * self::median($absDev);
        $robustVar = $robustSd * $robustSd;
        $sampleVar = $this->params->varianceBlend * $sampleVar + (1.0 - $this->params->varianceBlend) * $robustVar;

        $priorVar = $this->params->withinSdPrior ** 2;
        $priorWeight = $this->params->priorWeight;
        // Trust the user's data more when observed variance clearly exceeds the
        // prior; this widens intervals for genuinely variable users instead of
        // shrinking them back to the population prior.
        if ($sampleVar > $priorVar) {
            $priorWeight = max(0.5, $this->params->priorWeight * $priorVar / $sampleVar);
        }
        $withinVar = ($n * $sampleVar + $priorWeight * $priorVar) / ($n + $priorWeight);
        $withinVar = min(max($withinVar, $this->params->minWithinVar), $this->params->maxWithinVar);

        $betweenVar = $this->params->popBetweenSd ** 2;
        $weight = ($n * $betweenVar) / max(1e-9, $n * $betweenVar + $withinVar);
        $shrunkMean = $weight * $mean + (1.0 - $weight) * $this->params->popMean;
        $shrunkMean = min(max($shrunkMean, $this->params->minMean), $this->params->maxMean);

        $recent = array_slice($lengths, -6);
        $trend = 0.0;
        if ($n >= 4 && count($recent) >= 3) {
            $weights = [];
            $w = 1.0;
            foreach (array_reverse($recent) as $ignored) {
                $weights[] = $w;
                $w *= (1.0 - $this->params->trendEwma);
            }
            $weights = array_reverse($weights);
            $totalW = array_sum($weights);
            $acc = 0.0;
            foreach ($recent as $i => $y) {
                $acc += ($y - $shrunkMean) * $weights[$i];
            }
            $trend = $acc / $totalW;
            $shift = abs($trend) < $this->params->trendThreshold
                ? 0.0
                : max(-$this->params->trendMaxShift, min($this->params->trendMaxShift, $trend));
            $shrunkMean = min(max($shrunkMean + $shift, $this->params->minMean), $this->params->maxMean);
        }

        $recentMean = array_sum($recent) / count($recent);
        $recentSumSq = 0.0;
        foreach ($recent as $y) {
            $recentSumSq += ($y - $recentMean) ** 2;
        }
        $recentSd = sqrt($recentSumSq / max(1, count($recent) - 1));
        $recentRange = (int) (max($recent) - min($recent));
        $gap = max($recent) > $this->params->abstainGapDays;

        $pi = $this->params->tailProbBase
            + $this->params->tailProbSdGain * max(0.0, $recentSd - 2.0)
            + $this->params->tailProbRangeGain * max(0.0, $recentRange - 7);
        if ($gap) {
            $pi += $this->params->tailProbGapBonus;
        }
        $pi = min(max($pi, $this->params->tailProbBase), $this->params->tailProbMax);

        $abstain = $gap || $recentSd > $this->params->abstainRecentSd || $n < $this->params->abstainMinCycles;
        $confidence = $abstain
            ? 'abstain'
            : ($n >= 6 ? 'established' : ($n >= 3 ? 'building' : ($n >= 1 ? 'learning' : 'new')));

        return [
            'mean' => $shrunkMean,
            'within_var' => $withinVar,
            'cycle_sd' => sqrt($withinVar + ($withinVar * $betweenVar) / max(1e-9, $n * $betweenVar + $withinVar)),
            'n' => $n,
            'tail_prob' => $pi,
            'trend' => $trend,
            'recent_sd' => $recentSd,
            'recent_range' => $recentRange,
            'gap' => $gap,
            'abstain' => $abstain,
            'confidence' => $confidence,
        ];
    }

    /** @param array<string, mixed> $state */
    public function predictive(array $state): CalibratedPredictive
    {
        $sd = ((float) $state['cycle_sd']) * $this->params->bodyScale;
        $base = LognormalMixture::fromMoments(
            mean: (float) $state['mean'],
            sd: $sd,
            pi: (float) $state['tail_prob'],
            tailShift: $this->params->tailShiftDays,
            tailSdMultiplier: $this->params->tailSdMultiplier,
            momentMatch: $this->params->momentMatch,
        );
        return CalibratedPredictive::fromBase($base, $this->params);
    }

    /**
     * @param  list<int|float>  $lengths
     * @return array<string, mixed>
     */
    public function forecast(array $lengths, int $currentDay = 1): array
    {
        $state = $this->personalState($lengths);
        $dist = $this->predictive($state);
        $hazard = $this->onsetProbabilities($dist, $currentDay);
        $ovulation = $this->ovulationDistribution($dist);

        $quantiles = [];
        foreach ([0.025, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.975] as $q) {
            $quantiles[(string) $q] = $dist->quantile($q);
        }

        $explanations = [];
        if ($state['abstain']) {
            $explanations[] = 'Insufficient or highly variable history for a confident forecast.';
        }
        if ($state['gap']) {
            $explanations[] = 'A long gap was detected; the forecast is intentionally wide.';
        }
        if ($state['recent_sd'] > 4.0) {
            $explanations[] = 'Recent cycle lengths are unusually variable; the interval is wide.';
        }
        if ($state['n'] < 3) {
            $explanations[] = 'Only a few cycles recorded; the population prior dominates.';
        }
        if ($state['trend'] != 0.0) {
            $explanations[] = 'A recent trend was applied to the expected length.';
        }

        return [
            'params_version' => $this->params->version,
            'mean_length' => (float) $state['mean'],
            'sd_length' => ((float) $state['cycle_sd']) * $this->params->bodyScale,
            'cycles_used' => (int) $state['n'],
            'confidence' => (string) $state['confidence'],
            'abstain' => (bool) $state['abstain'],
            'abstain_reason' => $this->abstainReason($state),
            'tail_prob' => (float) $state['tail_prob'],
            'trend' => (float) $state['trend'],
            'quantiles' => $quantiles,
            'hazard' => $hazard,
            'ovulation' => [
                'mean_day' => $ovulation['mean_day'],
                'sd_day' => $ovulation['sd_day'],
                'fertile_probabilities' => $ovulation['fertile_probabilities'],
                'peak_days' => $ovulation['peak_days'],
                'window_start' => $ovulation['window_start'],
                'window_end' => $ovulation['window_end'],
            ],
            'prompt' => $this->prompt($state, $hazard),
            'explanations' => $explanations,
        ];
    }

    /** @return array<int, float> day => hazard */
    public function hazardCurve(CalibratedPredictive $dist): array
    {
        $curve = [];
        $prevSurvival = 1.0;
        for ($day = 1; $day <= $this->params->hazardMaxDay; $day++) {
            if ($prevSurvival <= 1e-12) {
                $curve[$day] = 1.0;
                continue;
            }
            $f = $dist->cdf((float) $day) - $dist->cdf((float) ($day - 1));
            $h = max(0.0, $f) / $prevSurvival;
            $curve[$day] = min(max($h, 0.0), 1.0);
            $prevSurvival -= max(0.0, $f);
            if ($prevSurvival < 0.0) {
                $prevSurvival = 0.0;
            }
        }
        return $curve;
    }

    /** @return array<string, mixed> */
    public function onsetProbabilities(CalibratedPredictive $dist, int $currentDay): array
    {
        $currentDay = max(1, $currentDay);
        $survivalNow = max(1e-12, 1.0 - $dist->cdf((float) ($currentDay - 1)));
        $within = function (int $k) use ($dist, $currentDay, $survivalNow): float {
            $survivalThen = max(0.0, 1.0 - $dist->cdf((float) ($currentDay - 1 + $k)));
            return min(max(1.0 - $survivalThen / $survivalNow, 0.0), 1.0);
        };
        $curve = $this->hazardCurve($dist);
        $pToday = $curve[$currentDay] ?? 0.0;
        $window = [];
        for ($d = $currentDay; $d <= min($this->params->hazardMaxDay, $currentDay + 13); $d++) {
            $window[$d] = $curve[$d] ?? 0.0;
        }
        return [
            'current_day' => $currentDay,
            'p_today' => $pToday,
            'p_within_3' => $within(3),
            'p_within_7' => $within(7),
            'p_within_14' => $within(14),
            'late_probability' => 1.0 - $dist->cdf((float) ($currentDay - 1)),
            'curve' => $window,
        ];
    }

    /** @return array<string, mixed> */
    public function ovulationDistribution(CalibratedPredictive $dist): array
    {
        $sdO = sqrt(max($dist->sd() ** 2 - $this->params->lutealSd ** 2, $this->params->ovulationSdFloor ** 2));
        $meanO = $dist->mean() - $this->params->lutealMean;

        $step = $this->params->ovulationGridStep;
        $maxY = 365.0;
        $edges = [];
        $y = 0.0;
        while ($y <= $maxY) {
            $edges[] = $y;
            $y += $step;
        }
        $cdfEdges = [];
        foreach ($edges as $edge) {
            $cdfEdges[] = $dist->cdf($edge);
        }

        $cdfO = function (float $d) use ($edges, $cdfEdges): float {
            $total = 0.0;
            $count = count($edges);
            for ($i = 0; $i < $count - 1; $i++) {
                $mass = $cdfEdges[$i + 1] - $cdfEdges[$i];
                if ($mass <= 0.0) {
                    continue;
                }
                $yMid = 0.5 * ($edges[$i] + $edges[$i + 1]);
                $z = ($d - $yMid + $this->params->lutealMean) / $this->params->lutealSd;
                $total += $mass * Erf::normCdf($z);
            }
            return $total;
        };

        $probs = [];
        $prev = 0.0;
        for ($day = 1; $day <= $this->params->hazardMaxDay; $day++) {
            $cur = $cdfO((float) $day);
            $probs[$day] = min(max($cur - $prev, 0.0), 1.0);
            $prev = $cur;
        }
        $probs[$this->params->hazardMaxDay] += max(0.0, 1.0 - $prev);

        $fertile = [];
        for ($day = 1; $day <= $this->params->hazardMaxDay; $day++) {
            $window = 0.0;
            $last = min($this->params->hazardMaxDay, $day + $this->params->fertileWindowDays - 1);
            for ($o = $day; $o <= $last; $o++) {
                $window += $probs[$o];
            }
            $fertile[$day] = min(max($window, 0.0), 1.0);
        }

        $maxProb = max($fertile);
        $peakDays = [];
        foreach ($fertile as $day => $p) {
            if ($p >= max(0.25, 0.5 * $maxProb)) {
                $peakDays[] = $day;
            }
        }
        $windowDays = [];
        foreach ($fertile as $day => $p) {
            if ($p >= 0.05) {
                $windowDays[] = $day;
            }
        }

        return [
            'mean_day' => $meanO,
            'sd_day' => $sdO,
            'probabilities' => $probs,
            'fertile_probabilities' => $fertile,
            'peak_days' => $peakDays,
            'window_start' => $windowDays === [] ? null : min($windowDays),
            'window_end' => $windowDays === [] ? null : max($windowDays),
        ];
    }

    /** @param array<string, mixed> $state */
    private function abstainReason(array $state): ?string
    {
        if (! $state['abstain']) {
            return null;
        }
        if ($state['gap']) {
            return 'long_gap';
        }
        if ($state['recent_sd'] > $this->params->abstainRecentSd) {
            return 'high_variability';
        }
        return 'insufficient_history';
    }

    /**
     * @param  list<int|float>  $lengths
     * @return list<float>
     */
    private function filterLengths(array $lengths): array
    {
        $out = [];
        foreach ($lengths as $length) {
            $y = (float) $length;
            if ($y >= $this->params->minLength && $y <= $this->params->maxLength) {
                $out[] = $y;
            }
        }
        return $out;
    }

    /** @param list<int|float> $values */
    private static function median(array $values): float
    {
        $n = count($values);
        if ($n === 0) {
            return 0.0;
        }
        sort($values);
        return $n % 2 === 1 ? (float) $values[intdiv($n, 2)] : 0.5 * ((float) $values[$n / 2 - 1] + (float) $values[$n / 2]);
    }

    /** @param array<string, mixed> $state @param array<string, mixed> $hazard */
    private function prompt(array $state, array $hazard): string
    {
        if ($state['abstain']) {
            return 'Not enough signal for a reliable forecast.';
        }
        $p3 = (float) $hazard['p_within_3'];
        if ($p3 < 0.05) {
            return 'Your period is unlikely in the next three days.';
        }
        if ($p3 < 0.25) {
            return 'There is a small chance your period starts in the next three days.';
        }
        if ($p3 < 0.5) {
            return 'Your period is fairly likely within three days.';
        }
        if ($p3 < 0.8) {
            return 'Your period is likely within three days.';
        }
        return 'Your period is very likely within three days.';
    }
}

	<?php

	declare(strict_types=1);

	namespace Dochia;

	/**
	* Dochia v1 — reference implementation (pure PHP, no deps).
	*
	* Pure PHP reference implementation; mirrors python/dochia.py. The golden
	* fixtures in tests/fixtures/golden.json bind the two implementations.
	*/
	final class MicroModel
	{
	public function __construct(public readonly ModelParams $params)
	{
	}

	/**
	* Personalise a user's cycle-length distribution from observed lengths.
	*
	* @param list<int\|float> $lengths
	* @return array<string, mixed>
	*/
	public function personalState(array $lengths): array
	{
	$lengths = $this->filterLengths($lengths);
	$n = count($lengths);

	if ($n === 0) {
	$withinVar = $this->params->withinSdPrior ** 2;
	$betweenVar = $this->params->popBetweenSd ** 2;
	$posteriorVar = ($withinVar * $betweenVar) / max(1e-9, $withinVar);
	return [
	'mean' => $this->params->popMean,
	'within_var' => $withinVar,
	'cycle_sd' => sqrt($withinVar + $posteriorVar),
	'n' => 0,
	'tail_prob' => $this->params->tailProbBase,
	'trend' => 0.0,
	'recent_sd' => 0.0,
	'recent_range' => 0,
	'gap' => false,
	'abstain' => false,
	'confidence' => 'new',
	];
	}

	$mean = array_sum($lengths) / $n;
	$sampleVar = 0.0;
	if ($n > 1) {
	$sumSq = 0.0;
	foreach ($lengths as $y) {
	$sumSq += ($y - $mean) ** 2;
	}
	$sampleVar = $sumSq / ($n - 1);
	}
	$median = self::median($lengths);
	$absDev = [];
	foreach ($lengths as $y) {
	$absDev[] = abs($y - $median);
	}
	$robustSd = 1.4826 * self::median($absDev);
	$robustVar = $robustSd * $robustSd;
	$sampleVar = $this->params->varianceBlend * $sampleVar + (1.0 - $this->params->varianceBlend) * $robustVar;

	$priorVar = $this->params->withinSdPrior ** 2;
	$priorWeight = $this->params->priorWeight;
	// Trust the user's data more when observed variance clearly exceeds the
	// prior; this widens intervals for genuinely variable users instead of
	// shrinking them back to the population prior.
	if ($sampleVar > $priorVar) {
	$priorWeight = max(0.5, $this->params->priorWeight * $priorVar / $sampleVar);
	}
	$withinVar = ($n * $sampleVar + $priorWeight * $priorVar) / ($n + $priorWeight);
	$withinVar = min(max($withinVar, $this->params->minWithinVar), $this->params->maxWithinVar);

	$betweenVar = $this->params->popBetweenSd ** 2;
	$weight = ($n * $betweenVar) / max(1e-9, $n * $betweenVar + $withinVar);
	$shrunkMean = $weight * $mean + (1.0 - $weight) * $this->params->popMean;
	$shrunkMean = min(max($shrunkMean, $this->params->minMean), $this->params->maxMean);

	$recent = array_slice($lengths, -6);
	$trend = 0.0;
	if ($n >= 4 && count($recent) >= 3) {
	$weights = [];
	$w = 1.0;
	foreach (array_reverse($recent) as $ignored) {
	$weights[] = $w;
	$w *= (1.0 - $this->params->trendEwma);
	}
	$weights = array_reverse($weights);
	$totalW = array_sum($weights);
	$acc = 0.0;
	foreach ($recent as $i => $y) {
	$acc += ($y - $shrunkMean) * $weights[$i];
	}
	$trend = $acc / $totalW;
	$shift = abs($trend) < $this->params->trendThreshold
	? 0.0
	: max(-$this->params->trendMaxShift, min($this->params->trendMaxShift, $trend));
	$shrunkMean = min(max($shrunkMean + $shift, $this->params->minMean), $this->params->maxMean);
	}

	$recentMean = array_sum($recent) / count($recent);
	$recentSumSq = 0.0;
	foreach ($recent as $y) {
	$recentSumSq += ($y - $recentMean) ** 2;
	}
	$recentSd = sqrt($recentSumSq / max(1, count($recent) - 1));
	$recentRange = (int) (max($recent) - min($recent));
	$gap = max($recent) > $this->params->abstainGapDays;

	$pi = $this->params->tailProbBase
	+ $this->params->tailProbSdGain * max(0.0, $recentSd - 2.0)
	+ $this->params->tailProbRangeGain * max(0.0, $recentRange - 7);
	if ($gap) {
	$pi += $this->params->tailProbGapBonus;
	}
	$pi = min(max($pi, $this->params->tailProbBase), $this->params->tailProbMax);

	$abstain = $gap \|\| $recentSd > $this->params->abstainRecentSd \|\| $n < $this->params->abstainMinCycles;
	$confidence = $abstain
	? 'abstain'
	: ($n >= 6 ? 'established' : ($n >= 3 ? 'building' : ($n >= 1 ? 'learning' : 'new')));

	return [
	'mean' => $shrunkMean,
	'within_var' => $withinVar,
	'cycle_sd' => sqrt($withinVar + ($withinVar * $betweenVar) / max(1e-9, $n * $betweenVar + $withinVar)),
	'n' => $n,
	'tail_prob' => $pi,
	'trend' => $trend,
	'recent_sd' => $recentSd,
	'recent_range' => $recentRange,
	'gap' => $gap,
	'abstain' => $abstain,
	'confidence' => $confidence,
	];
	}

	/** @param array<string, mixed> $state */
	public function predictive(array $state): CalibratedPredictive
	{
	$sd = ((float) $state['cycle_sd']) * $this->params->bodyScale;
	$base = LognormalMixture::fromMoments(
	mean: (float) $state['mean'],
	sd: $sd,
	pi: (float) $state['tail_prob'],
	tailShift: $this->params->tailShiftDays,
	tailSdMultiplier: $this->params->tailSdMultiplier,
	momentMatch: $this->params->momentMatch,
	);
	return CalibratedPredictive::fromBase($base, $this->params);
	}

	/**
	* @param list<int\|float> $lengths
	* @return array<string, mixed>
	*/
	public function forecast(array $lengths, int $currentDay = 1): array
	{
	$state = $this->personalState($lengths);
	$dist = $this->predictive($state);
	$hazard = $this->onsetProbabilities($dist, $currentDay);
	$ovulation = $this->ovulationDistribution($dist);

	$quantiles = [];
	foreach ([0.025, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.975] as $q) {
	$quantiles[(string) $q] = $dist->quantile($q);
	}

	$explanations = [];
	if ($state['abstain']) {
	$explanations[] = 'Insufficient or highly variable history for a confident forecast.';
	}
	if ($state['gap']) {
	$explanations[] = 'A long gap was detected; the forecast is intentionally wide.';
	}
	if ($state['recent_sd'] > 4.0) {
	$explanations[] = 'Recent cycle lengths are unusually variable; the interval is wide.';
	}
	if ($state['n'] < 3) {
	$explanations[] = 'Only a few cycles recorded; the population prior dominates.';
	}
	if ($state['trend'] != 0.0) {
	$explanations[] = 'A recent trend was applied to the expected length.';
	}

	return [
	'params_version' => $this->params->version,
	'mean_length' => (float) $state['mean'],
	'sd_length' => ((float) $state['cycle_sd']) * $this->params->bodyScale,
	'cycles_used' => (int) $state['n'],
	'confidence' => (string) $state['confidence'],
	'abstain' => (bool) $state['abstain'],
	'abstain_reason' => $this->abstainReason($state),
	'tail_prob' => (float) $state['tail_prob'],
	'trend' => (float) $state['trend'],
	'quantiles' => $quantiles,
	'hazard' => $hazard,
	'ovulation' => [
	'mean_day' => $ovulation['mean_day'],
	'sd_day' => $ovulation['sd_day'],
	'fertile_probabilities' => $ovulation['fertile_probabilities'],
	'peak_days' => $ovulation['peak_days'],
	'window_start' => $ovulation['window_start'],
	'window_end' => $ovulation['window_end'],
	],
	'prompt' => $this->prompt($state, $hazard),
	'explanations' => $explanations,
	];
	}

	/** @return array<int, float> day => hazard */
	public function hazardCurve(CalibratedPredictive $dist): array
	{
	$curve = [];
	$prevSurvival = 1.0;
	for ($day = 1; $day <= $this->params->hazardMaxDay; $day++) {
	if ($prevSurvival <= 1e-12) {
	$curve[$day] = 1.0;
	continue;
	}
	$f = $dist->cdf((float) $day) - $dist->cdf((float) ($day - 1));
	$h = max(0.0, $f) / $prevSurvival;
	$curve[$day] = min(max($h, 0.0), 1.0);
	$prevSurvival -= max(0.0, $f);
	if ($prevSurvival < 0.0) {
	$prevSurvival = 0.0;
	}
	}
	return $curve;
	}

	/** @return array<string, mixed> */
	public function onsetProbabilities(CalibratedPredictive $dist, int $currentDay): array
	{
	$currentDay = max(1, $currentDay);
	$survivalNow = max(1e-12, 1.0 - $dist->cdf((float) ($currentDay - 1)));
	$within = function (int $k) use ($dist, $currentDay, $survivalNow): float {
	$survivalThen = max(0.0, 1.0 - $dist->cdf((float) ($currentDay - 1 + $k)));
	return min(max(1.0 - $survivalThen / $survivalNow, 0.0), 1.0);
	};
	$curve = $this->hazardCurve($dist);
	$pToday = $curve[$currentDay] ?? 0.0;
	$window = [];
	for ($d = $currentDay; $d <= min($this->params->hazardMaxDay, $currentDay + 13); $d++) {
	$window[$d] = $curve[$d] ?? 0.0;
	}
	return [
	'current_day' => $currentDay,
	'p_today' => $pToday,
	'p_within_3' => $within(3),
	'p_within_7' => $within(7),
	'p_within_14' => $within(14),
	'late_probability' => 1.0 - $dist->cdf((float) ($currentDay - 1)),
	'curve' => $window,
	];
	}

	/** @return array<string, mixed> */
	public function ovulationDistribution(CalibratedPredictive $dist): array
	{
	$sdO = sqrt(max($dist->sd() 2 - $this->params->lutealSd 2, $this->params->ovulationSdFloor ** 2));
	$meanO = $dist->mean() - $this->params->lutealMean;

	$step = $this->params->ovulationGridStep;
	$maxY = 365.0;
	$edges = [];
	$y = 0.0;
	while ($y <= $maxY) {
	$edges[] = $y;
	$y += $step;
	}
	$cdfEdges = [];
	foreach ($edges as $edge) {
	$cdfEdges[] = $dist->cdf($edge);
	}

	$cdfO = function (float $d) use ($edges, $cdfEdges): float {
	$total = 0.0;
	$count = count($edges);
	for ($i = 0; $i < $count - 1; $i++) {
	$mass = $cdfEdges[$i + 1] - $cdfEdges[$i];
	if ($mass <= 0.0) {
	continue;
	}
	$yMid = 0.5 * ($edges[$i] + $edges[$i + 1]);
	$z = ($d - $yMid + $this->params->lutealMean) / $this->params->lutealSd;
	$total += $mass * Erf::normCdf($z);
	}
	return $total;
	};

	$probs = [];
	$prev = 0.0;
	for ($day = 1; $day <= $this->params->hazardMaxDay; $day++) {
	$cur = $cdfO((float) $day);
	$probs[$day] = min(max($cur - $prev, 0.0), 1.0);
	$prev = $cur;
	}
	$probs[$this->params->hazardMaxDay] += max(0.0, 1.0 - $prev);

	$fertile = [];
	for ($day = 1; $day <= $this->params->hazardMaxDay; $day++) {
	$window = 0.0;
	$last = min($this->params->hazardMaxDay, $day + $this->params->fertileWindowDays - 1);
	for ($o = $day; $o <= $last; $o++) {
	$window += $probs[$o];
	}
	$fertile[$day] = min(max($window, 0.0), 1.0);
	}

	$maxProb = max($fertile);
	$peakDays = [];
	foreach ($fertile as $day => $p) {
	if ($p >= max(0.25, 0.5 * $maxProb)) {
	$peakDays[] = $day;
	}
	}
	$windowDays = [];
	foreach ($fertile as $day => $p) {
	if ($p >= 0.05) {
	$windowDays[] = $day;
	}
	}

	return [
	'mean_day' => $meanO,
	'sd_day' => $sdO,
	'probabilities' => $probs,
	'fertile_probabilities' => $fertile,
	'peak_days' => $peakDays,
	'window_start' => $windowDays === [] ? null : min($windowDays),
	'window_end' => $windowDays === [] ? null : max($windowDays),
	];
	}

	/** @param array<string, mixed> $state */
	private function abstainReason(array $state): ?string
	{
	if (! $state['abstain']) {
	return null;
	}
	if ($state['gap']) {
	return 'long_gap';
	}
	if ($state['recent_sd'] > $this->params->abstainRecentSd) {
	return 'high_variability';
	}
	return 'insufficient_history';
	}

	/**
	* @param list<int\|float> $lengths
	* @return list<float>
	*/
	private function filterLengths(array $lengths): array
	{
	$out = [];
	foreach ($lengths as $length) {
	$y = (float) $length;
	if ($y >= $this->params->minLength && $y <= $this->params->maxLength) {
	$out[] = $y;
	}
	}
	return $out;
	}

	/** @param list<int\|float> $values */
	private static function median(array $values): float
	{
	$n = count($values);
	if ($n === 0) {
	return 0.0;
	}
	sort($values);
	return $n % 2 === 1 ? (float) $values[intdiv($n, 2)] : 0.5 * ((float) $values[$n / 2 - 1] + (float) $values[$n / 2]);
	}

	/** @param array<string, mixed> $state @param array<string, mixed> $hazard */
	private function prompt(array $state, array $hazard): string
	{
	if ($state['abstain']) {
	return 'Not enough signal for a reliable forecast.';
	}
	$p3 = (float) $hazard['p_within_3'];
	if ($p3 < 0.05) {
	return 'Your period is unlikely in the next three days.';
	}
	if ($p3 < 0.25) {
	return 'There is a small chance your period starts in the next three days.';
	}
	if ($p3 < 0.5) {
	return 'Your period is fairly likely within three days.';
	}
	if ($p3 < 0.8) {
	return 'Your period is likely within three days.';
	}
	return 'Your period is very likely within three days.';
	}
	}