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

class oneVariance extends hypothesis
{
    public $hypo;

    public $hypoVal;

    public $cl;

    public $α;

    public $alterHypo;

    public $desc;

    public $test;

    public $isContinuous;

    public $dataframe;

    /**
     * constructor
     *
     * @param object $dataframe
     * @param string $cl
     * @param string $hypo
     * @param string $hypoVal
     * @param string $alterHypo
     * @access public
     * @return void
     */
    public function __construct($dataframe, $cl, $hypo, $hypoVal, $alterHypo)
    {
        $this->hypo      = $hypo;
        $this->hypoVal   = $hypoVal;
        $this->cl        = $cl;
        $this->α         = 1 - $cl / 100;
        $this->alterHypo = $alterHypo;

        $this->desc = array();
        $this->test = array();

        parent::__construct($dataframe);
    }

    /**
     * Set Data of sample column.
     *
     * @param array $sampleCol
     * @access public
     * @return void
     */
    public function setData($sampleCol)
    {
        $this->isContinuous = true;
        $this->desc = array();
        $this->test = array();

        $cols = $this->dataframe->cols($sampleCol, 'number');
        foreach($cols as $col)
        {
            $name     = $col->name;
            $variance = Describe::variance($col->trimdata, false);
            $standard = Describe::standard($col->trimdata, false);
            $N        = count($col->trimdata);

            if(!isset($this->desc[$col->name])) $this->desc[$name] = new stdclass();
            if(!isset($this->test[$col->name])) $this->test[$name] = new stdclass();

            $this->desc[$name]->variable = $name;
            $this->desc[$name]->std      = $standard;
            $this->desc[$name]->variance = $variance;
            $this->desc[$name]->N        = $N;
            $this->desc[$name]->data     = $col->trimdata;

            $this->test[$name]->variable = $name;
            $this->test[$name]->freedom  = $N - 1;
            $this->test[$name]->std      = $standard;
            $this->test[$name]->variance = $variance;
            $this->test[$name]->data     = $col->trimdata;
            $this->test[$name]->N        = $N;
        }
    }

    /**
     * SetData of sample summary.
     *
     * @param string $N
     * @param string $type
     * @param string $value
     * @access public
     * @return void
     */
    public function setSummary($N, $type, $value)
    {
        $this->isContinuous = false;
        $this->desc = array();
        $this->test = array();

        $desc    = new stdclass();
        $test    = new stdclass();
        $desc->N = $N;
        if($type == 'variance')
        {
            $desc->variance = $value;
            $desc->std      = sqrt($value);
            $test->variance = $value ** 2;
            $test->std      = sqrt($value);
        }

        if($type == 'standard')
        {
            $desc->std      = $value;
            $desc->variance = $value ** 2;
            $test->std      = $value;
            $test->variance = $value ** 2;
        }

        $this->desc[] = $desc;

        $test->freedom = $N - 1;
        $this->test[] = $test;
    }

    /**
     * Get describe values.
     *
     * @param string $type
     * @param bool $round
     * @access public
     * @return void
     */
    public function getDesc($type, $round = true)
    {
        $result = array();
        foreach($this->desc as $desc)
        {
            $chisqCI = $this->chisqCI($desc->N, $desc->variance, $type, $round);
            $desc->chisqCI  = sprintf('(%s, %s)', $chisqCI[0], $chisqCI[1]);
            if($this->isContinuous)
            {
                $bonettCI = $this->bonettCI($desc->data, $desc->N, $desc->std, $type, $round);
                $desc->bonettCI = sprintf('(%s, %s)', $bonettCI[0], $bonettCI[1]);
            }
            if($round)
            {
                $desc->variance = $this->roundDigit($desc->variance);
                $desc->std      = $this->roundDigit($desc->std);
            }

            $result[] = $desc;
        }

        return $result;
    }

    /**
     * Calculate confidence interval of bonett.
     * see https://support.minitab.com/zh-cn/minitab/21/help-and-how-to/statistics/basic-statistics/how-to/1-variance/methods-and-formulas/methods-and-formulas/#confidence-intervals-and-bounds-for-the-bonett-method
     *
     * @param array $data
     * @param string $N
     * @param string $std
     * @param string $type
     * @param bool $round
     * @access public
     * @return void
     */
    public function bonettCI($data, $N, $std, $type, $round = true)
    {
        $se   = $this->se($N, $data, $std);
        $stdN = new standardNormal();

        $α = $this->alterHypo == 'ne' ? 1 - $this->α / 2 : 1 - $this->α;
        $z = $stdN->inverse($α);
        $c = $N / ($N - $z);

        $part1 = $c * $std ** 2;
        $part2 = $z * $se;
        $eLow = exp(-1 * $c * $part2);
        $eUp  = exp($c * $part2);

        //              -Cα · Zα · se                        Cα · Zα · se
        // Cα · S² · e                  ≤  σ² ≤   Cα/2 · S² · e
        $lowLimit = $part1 * $eLow;
        $upLimit  = $part1 * $eUp;

        return $this->formatCI($lowLimit, $upLimit, $type, $round);
    }

    /**
     * Calculate se of bonett.
     * see https://support.minitab.com/zh-cn/minitab/21/help-and-how-to/statistics/basic-statistics/how-to/1-variance/methods-and-formulas/methods-and-formulas/#confidence-intervals-and-bounds-for-the-bonett-method
     *
     *          N      n   / xᵢ - m \⁴
     *  γ =  ————————  Σ  ( ———————— )
     *       (N - 1)² i=1  \    s   /
     *
     *          ————————————
     *         /     N - 3
     *        / γ - ———————
     *  se = /         N
     *      /  ——————————————
     *     √       N - 1
     *
     * @param array $data
     * @param string $N
     * @param string $std
     * @param string $type
     * @param bool $round
     * @access public
     * @return void
     */
    public function se($N, $data, $std)
    {
        if($N <= 5) $m = Describe::mean($data);
        if($N > 5)
        {
            $mPercent = 1 / (2 * sqrt($N - 4));
            $m = Describe::trimMean($data, $mPercent);
        }

        $γ = 0;
        foreach($data as $value)
        {
            $γ += (($value - $m) / $std) ** 4;
        }

        $γ *= ($N / ($N - 1) ** 2);

        $se = sqrt(($γ - (($N - 3) / $N)) / ($N - 1));

        return $se;
    }

    /**
     * Calculate confidence interval of chi-squared.
     * see https://support.minitab.com/zh-cn/minitab/21/help-and-how-to/statistics/basic-statistics/how-to/1-variance/methods-and-formulas/methods-and-formulas/#confidence-intervals-and-bounds-for-the-chi-square-method
     *
     * @param string $N
     * @param string $var
     * @param string $type
     * @param bool $round
     * @access public
     * @return void
     */
    public function chisqCI($N, $var, $type, $round = true)
    {
        $chisq = new ChiSquare($N - 1);
        $α = $this->alterHypo == 'ne' ? $this->α / 2 : $this->α;
        $lowInv = $chisq->inverse(1 - $α);
        $upInv  = $chisq->inverse($α);

        $part = ($N - 1) * $var;

        $lowLimit = $lowInv == 0 ? 0 : $part / $lowInv;
        $upLimit  = $upInv == 0 ? 0 : $part / $upInv;

        return $this->formatCI($lowLimit, $upLimit, $type, $round);
    }

    /**
     * Format CI result to array.
     *
     * @param float $lowLimit
     * @param float $upLimit
     * @param string $type
     * @param bool $round
     * @access public
     * @return void
     */
    public function formatCI($lowLimit, $upLimit, $type, $round)
    {
        if($type != 'variance')
        {
            $lowLimit = sqrt($lowLimit);
            $upLimit  = sqrt($upLimit);
        }

        if($round)
        {
            $lowLimit = $this->roundDigit($lowLimit);
            $upLimit  = $this->roundDigit($upLimit);
        }

        if($this->alterHypo == 'lt')
        {
            return array('-', $upLimit);
        }

        if($this->alterHypo == 'gt')
        {
            return array($lowLimit, '-');
        }

        return array($lowLimit, $upLimit);
    }

    /**
     * Get Test result.
     *
     * @param bool $round
     * @access public
     * @return void
     */
    public function getTest($round = true)
    {
        $result = array();
        foreach($this->test as $test)
        {
            if($this->isContinuous) $result[] = $this->bonettTest($test, $round);
            $result[] = $this->chisqTest($test, $round);
        }

        return $result;
    }

    /**
     * Calculate test result of bonett.
     * see https://support.minitab.com/zh-cn/minitab/21/help-and-how-to/statistics/basic-statistics/how-to/1-variance/methods-and-formulas/methods-and-formulas/#hypothesis-test-for-the-bonett-method
     *
     * @param object $test
     * @param bool $round
     * @access public
     * @return void
     */
    public function bonettTest($test, $round = true)
    {
        $σ₀² = $this->hypo == 'variance' ? $this->hypoVal : $this->hypoVal ** 2;
        $N = $test->N;
        $data = $test->data;
        $std = $test->std;

        $se   = $this->se($N, $data, $std);
        $stdN = new standardNormal();

        $foundLow = false;
        $foundUp  = false;

        $aL = 0;
        $aU = 0;
        for($α = 0.001; $α <= 1; $α += 0.001)
        {
            $z = $stdN->inverse($α);
            $c = $N / ($N - $z);

            $part1 = $c * $std ** 2;
            $part2 = $z * $se;

            $low = $part1 * exp(-1 * $c * $part2);
            $up  = $part1 * exp($c * $part2);

            if(!$foundLow and $σ₀² >= $low)
            {
                $foundLow = true;
                $aL = $α;
            }
            if(!$foundUp and $σ₀² <= $up)
            {
                $foundUp = true;
                $aU = $α;
            }
            if($foundLow and $foundUp) break;
        }

        if($this->alterHypo == 'ne')
        {
            $aL = (1 - $aL) * 2;
            $aU = (1 - $aU) * 2;
            $p = min($aL, $aU);
        }
        if($this->alterHypo == 'lt') $p = $aU == 0 ? 0 : 1 - $aU;
        if($this->alterHypo == 'gt') $p = $aL == 0 ? 0 : 1 - $aL;

        if($p > 1) $p = 2 - $p;

        $obj = clone($test);
        $obj->method  = 'bonett';
        $obj->test    = '-';
        $obj->freedom = '-';
        $obj->p       = $this->roundDigit($p, $round);

        return $obj;
    }

    /**
     * Calculate test result of chi-squared.
     * see https://support.minitab.com/zh-cn/minitab/21/help-and-how-to/statistics/basic-statistics/how-to/1-variance/methods-and-formulas/methods-and-formulas/#hypothesis-test-for-the-chi-square-method
     *
     * @param object $test
     * @param bool $round
     * @access public
     * @return void
     */
    public function chisqTest($test, $round = true)
    {
        $σ₀² = $this->hypo == 'variance' ? $this->hypoVal : $this->hypoVal ** 2;
        $x² = $test->freedom * $test->variance / $σ₀²;

        $chisq = new ChiSquare($test->freedom);

        $cdf = $chisq->cdf($x²);
        $gtP = 1 - $cdf;
        $ltP = $cdf;
        $neP = 2 * min($gtP, $ltP);

        if($this->alterHypo == 'ne') $p = $neP;
        if($this->alterHypo == 'lt') $p = $ltP;
        if($this->alterHypo == 'gt') $p = $gtP;

        $obj         = clone($test);
        $obj->method = 'chisq';
        $obj->test   = $this->roundDigit($x², $round);
        $obj->p      = $this->roundDigit($p, $round);

        return $obj;
    }
}
