<?php
require_once dirname(__DIR__) . '/base/regression.php';
require_once dirname(dirname(__DIR__)) . '/zendasmath/basic/describe.php';
require_once dirname(dirname(__DIR__)) . '/zendasmath/distribution/stdnormal.php';
require_once dirname(dirname(__DIR__)) . '/zendasmath/distribution/studentT.php';

class Correlation extends RegressionBase
{
    public $a;

    public $b;

    public $n;

    public $cl;

    public $dataframe;

    /**
     * 构造方法。
     *
     * The construct function.
     *
     * @access public
     * @return void
     */
    public function __construct($dataframe, $points, $cl)
    {
        $this->a  = $points[0];
        $this->b  = $points[1];
        $this->n  = count($points[0]);
        $this->cl = $cl;

        parent::__construct($dataframe);
    }

    /**
     * Pearson correlation method.
     * https://support.minitab.com/zh-cn/minitab/21/help-and-how-to/statistics/basic-statistics/how-to/correlation/methods-and-formulas/methods-and-formulas/#pearson-s-correlation-coefficient
     *
     * @param array $cols
     * @access public
     * @return array
     */
    public function PearsonCorrelation($round = true)
    {
        $r  = $this->r($this->a, $this->b, false);
        $ci = $this->PearsonCI($r, $round);

        $result = new stdclass();
        $result->r   = $this->dataframe->roundDigit($r, $round);
        $result->p   = $this->p($r, $round);
        $result->p_l = $ci[0];
        $result->p_u = $ci[1];

        return $result;
    }

    public function PearsonCI($r, $round = true)
    {
        $standardNormal = new standardNormal();

        $α = $this->α();
        $z = $standardNormal->inverse(1 - $α / 2);

        $part1 = 0.5 * (log(1 + $r) - log(1 - $r));
        if($part1 == INF or $this->n <= 3) return array('-', '-');
        $part2 = $z / sqrt($this->n - 3);

        $z_l = $part1 - $part2;
        $z_u = $part1 + $part2;

        $expZ_L = exp(2 * $z_l);
        $expZ_U = exp(2 * $z_u);

        $p_l = ($expZ_L - 1) / ($expZ_L + 1);
        $p_u = ($expZ_U - 1) / ($expZ_U + 1);

        $p_l = $this->dataframe->roundDigit($p_l, $round);
        $p_u = $this->dataframe->roundDigit($p_u, $round);

        return array($p_l, $p_u);
    }

    public function SpearmanCorrelation($round = true)
    {
        $rankA = Describe::fractionalRanking($this->a);
        $rankB = Describe::fractionalRanking($this->b);
        $r     = $this->r($rankA, $rankB, false);
        $ci    = $this->SpearmanCI($r, $round);

        $result = new stdclass();
        $result->r   = $this->dataframe->roundDigit($r, $round);
        $result->p   = $this->p($r, $round);
        $result->p_l = $ci[0];
        $result->p_u = $ci[1];

        return $result;
    }

    public function SpearmanCI($r, $round = true)
    {
        $standardNormal = new standardNormal();

        $α = $this->α();
        $z = $standardNormal->inverse(1 - $α / 2);

        $se    = sqrt((1 + $r **2 / 2) / ($this->n - 3));
        $part1 = 0.5 * (log(1 + $r) - log(1 - $r));
        $part2 = $z * $se;

        $z_l = $part1 - $part2;
        $z_u = $part1 + $part2;

        $expZ_L = exp(2 * $z_l);
        $expZ_U = exp(2 * $z_u);

        $p_l = ($expZ_L - 1) / ($expZ_L + 1);
        $p_u = ($expZ_U - 1) / ($expZ_U + 1);

        $p_l = $this->dataframe->roundDigit($p_l, $round);
        $p_u = $this->dataframe->roundDigit($p_u, $round);

        return array($p_l, $p_u);
    }

    public function α()
    {
        return (100 - $this->cl) / 100;
    }

    public function r($a, $b, $round = true)
    {
        $σa = Describe::standard($a, false);
        $σb = Describe::standard($b, false);

        $meanA   = Describe::mean($a, false);
        $meanB   = Describe::mean($b, false);

        $sum = 0;
        foreach($a as $index => $aValue)
        {
            $bValue = $b[$index];
            $sum += ($aValue - $meanA) * ($bValue - $meanB);
        }
        $r = $sum / (($this->n - 1) * $σa * $σb);
        $r = $this->dataframe->roundDigit($r, $round);

        return $r;
    }

    public function p($r, $round)
    {
        $freedom  = $this->n - 2;
        $studentT = new StudentT($freedom);

        $part = sqrt(1 - $r ** 2);
        if($part == 0) return '-';

        $t   = $r * sqrt($this->n - 2) / $part;
        $cdf = $studentT->cdf($t);
        $p   = (1 - $cdf) * 2;
        if($p >= 1) $p = 2 - $p;

        $p = $this->dataframe->roundDigit($p, $round);

        return $p;
    }
}
