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

class MultipleRegression extends RegressionBase
{

    public $factors;

    public $response;

    public $coefN;

    public $observeN;

    public $isSimple;

    public $factor;

    public $factorX;

    public $factorY;

    public $dataframe;

    public function __construct($dataframe, $factors, $response)
    {
        $this->factors  = $factors;
        $this->response = $response;
        $this->coefN    = count($factors);
        $this->observeN = count($response);

        $this->isSimple = $this->coefN === 1;
        $this->factor   = current($factors);


        parent::__construct($dataframe);
    }

    public function getFourmula($xaxis, $yaxis)
    {
        $params  = $this->coef();
        $formula = array();
        foreach($params as $key => $param)
        {
            if($key == 0)
            {
                $formula[] = $this->dataframe->roundDigit($param);
            }
            else
            {
                $mark  = $param >= 0 ? ' + ' : ' - ';
                $param = abs($param);
                $formula[] = $mark . $this->dataframe->roundDigit($param) . ' ' . $xaxis[$key - 1];
            }
        }
        return '' . $yaxis . ' = ' . implode('', $formula);
    }

    public function checkInvertible()
    {
        if(!$this->isSimple)
        {
            $X = RegressionMatrix::designMatrix($this->factors);
            $Xt = $X->transpose(); // Xᵀ
            return $Xt->multiply($X)->isInvertible();
        }

        return true;
    }

    public function getFactor($round = true)
    {
        $X      = null;
        $Xt     = null;
        $XtXinv = null;
        if(!$this->isSimple)
        {
            $X = RegressionMatrix::designMatrix($this->factors);
            $Xt = $X->transpose(); // Xᵀ
            $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        }

        $ssError = $this->ssError();
        $dfError = $this->dfError();
        $s²  = $this->msError($ssError, $dfError);

        $coef = $this->coef($X, $Xt, $XtXinv);
        $se  = $this->sdError($s², $X, $Xt, $XtXinv);

        $t = $this->tFactor($coef, $se);
        $p = $this->pFactor($dfError, $t);

        $vif = $this->isSimple ? array(1) : $this->vif();
        array_unshift($vif, '');

        $result = new stdclass();
        $result->coef = array();
        $result->se   = array();
        $result->t    = array();
        $result->p    = array();
        $result->vif  = array();

        foreach($coef as $index => $c)
        {
            $result->coef[] = $this->dataframe->roundDigit($c);
            $result->se[]   = $this->dataframe->roundDigit($se[$index]);
            $result->t[]    = $this->dataframe->roundDigit($t[$index]);
            $result->p[]    = $this->dataframe->roundDigit($p[$index]);
            $result->vif[]  = $this->dataframe->roundDigit($vif[$index]);
        }

        return $result;
    }

    public function getModel($round = true)
    {
        $dfError = $this->dfError();
        $dfTotal = $this->dfTotal();
        $ssError = $this->ssError();
        $ssTotal = $this->ssTotal();
        $msError = $this->msError($ssError, $dfError);
        $msTotal = $this->msTotal($ssTotal, $dfTotal);

        $s      = is_numeric($msError) ? sqrt($msError) : '-';
        $rsq    = $this->rsq($ssError, $ssTotal);
        $rsqAdj = $this->rsqAdj($msError, $msTotal);
        $rsqPre = $this->rsqPre();

        if(is_numeric($rsq) and $rsq < 0) $rsq = 0;
        if(is_numeric($rsqAdj) and $rsqAdj < 0) $rsqAdj = 0;

        $s      = $this->dataframe->roundDigit($s, $round);
        $rsq    = $this->dataframe->roundDigit($rsq, $round);
        $rsqAdj = $this->dataframe->roundDigit($rsqAdj, $round);
        $rsqPre = $this->dataframe->roundDigit($rsqPre, $round);

        $result  = new stdclass();

        $result->s      = $s;
        $result->rsq    = is_numeric($rsq) ? $rsq * 100 : $rsq;
        $result->rsqAdj = is_numeric($rsqAdj) ? $rsqAdj * 100 : $rsqAdj;
        $result->rsqPre = is_numeric($rsqPre) ? $rsqPre * 100 : $rsqPre;

        return $result;
    }

    public function getAnova($round = true)
    {
        $df      = $this->df();
        $dfError = $this->dfError();
        $dfTotal = $this->dfTotal();

        $yFit = $this->fitY();

        $ssRegr  = $this->ssRegr($yFit);
        $ssError = $this->ssError($yFit);
        $ssTotal = $this->ssTotal();

        $msRegr  = $this->msRegr($ssRegr, $df);
        $msError = $this->msError($ssError, $dfError);
        $msTotal = $this->msTotal($ssTotal, $dfTotal);

        $fRegr = $this->fRegr($msRegr, $msError);
        $pRegr = $this->pRegr($df, $dfError, $fRegr);

        $same  = $this->sameRow();
        $dfPE  = '-';
        $dfLOF = '-';
        $ssPE  = '-';
        $ssLOF = '-';
        $msPE  = '-';
        $msLOF = '-';
        $fLOF  = '-';
        $pLOF  = '-';
        if(!empty($same))
        {
            $dfPE  = $this->dfPE($same);
            $dfLOF = $this->dfLOF($dfError, $dfPE);

            $ssPE  = $this->ssPE($same);
            $ssLOF = $this->ssLOF($ssError, $ssPE);

            $msPE = $this->msPE($ssPE, $dfPE);
            $msLOF = $this->msLOF($ssLOF, $dfLOF);

            if($msPE != 0)
            {
                $fLOF  = $this->fLOF($msLOF, $msPE);
                $pLOF  = $this->pLOF($dfLOF, $dfPE, $fLOF);
            }
        }

        if($this->isSimple)
        {
            $items = array(array($df, $ssRegr, $msRegr, $fRegr, $pRegr));
        }
        else
        {
            $items = $this->itemAnova($ssRegr, $msError, $dfError);
        }

        $result = new stdclass();

        $result->hasSame = !empty($same);

        $result->df      = $df;
        $result->dfError = $dfError;
        $result->dfLOF   = $dfLOF;
        $result->dfPE    = $dfPE;
        $result->dfTotal = $dfTotal;

        $result->ssRegr  = $this->dataframe->roundDigit($ssRegr, $round);
        $result->ssError = $this->dataframe->roundDigit($ssError, $round);
        $result->ssTotal = $this->dataframe->roundDigit($ssTotal, $round);
        $result->ssPE    = $this->dataframe->roundDigit($ssPE, $round);
        $result->ssLOF   = $this->dataframe->roundDigit($ssLOF, $round);

        $result->msRegr  = $this->dataframe->roundDigit($msRegr, $round);
        $result->msError = $this->dataframe->roundDigit($msError, $round);
        $result->msTotal = $this->dataframe->roundDigit($msTotal, $round);
        $result->msPE    = $this->dataframe->roundDigit($msPE, $round);
        $result->msLOF   = $this->dataframe->roundDigit($msLOF, $round);

        $result->fRegr = $this->dataframe->roundDigit($fRegr, $round);
        $result->pRegr = $this->dataframe->roundDigit($pRegr, $round);
        $result->fLOF  = $this->dataframe->roundDigit($fLOF, $round);
        $result->pLOF  = $this->dataframe->roundDigit($pLOF, $round);

        $result->itemDf = array();
        $result->itemSS = array();
        $result->itemMS = array();
        $result->itemF = array();
        $result->itemP = array();
        for($index = 0; $index < $this->coefN; $index++)
        {
            $item = $items[$index];

            $result->itemDf[] = $item[0];
            $result->itemSS[] = $this->dataframe->roundDigit($item[1], $round);
            $result->itemMS[] = $this->dataframe->roundDigit($item[2], $round);
            $result->itemF[]  = $this->dataframe->roundDigit($item[3], $round);
            $result->itemP[]  = $this->dataframe->roundDigit($item[4], $round);
        }

        return $result;
    }

    public function getTest($round = true)
    {
        $y        = $this->response;
        $yFit     = $this->fitY();
        $residual = $this->residual($yFit);
        $h        = RegressionMatrix::hMatrix($this->factors)->getDiagonalElements();
        $dfError  = $this->dfError();
        $ssError  = $this->ssError($yFit);
        $msError  = $this->msError($ssError, $dfError);

        $p_n = 3 * ($this->coefN + 1) / $this->observeN;
        $p_n = min($p_n, 0.99);

        $result = new stdclass();
        $result->observed = array();
        $result->value    = array();
        $result->fit      = array();
        $result->residual = array();
        $result->stdResid = array();
        $result->hasTest  = false;

        if(!is_numeric($ssError) or !is_numeric($msError) or $msError == 0) return $result;

        foreach($residual as $index => $e)
        {
            $setX = false;
            $setR = false;

            $stdResid = $h[$index] == 1 ? '-' : $e / sqrt($msError * (1 - $h[$index]));
            if($h[$index] > $p_n) $setX = true;
            if(abs($stdResid) > 2) $setR = true;

            if($setX || $setR)
            {
                $result->hasTest    = true;
                $result->observed[] = $index + 1;
                $result->value[]    = $y[$index];
                $result->fit[]      = $this->dataframe->roundDigit($yFit[$index], $round);
                $result->residual[] = $this->dataframe->roundDigit($e, $round);
                $stdResid = $this->dataframe->roundDigit($stdResid, $round);
                if($setR) $stdResid .= ' R';
                if($setX) $stdResid .= ' X';
                $result->stdResid[] = $stdResid;
            }
        }

        return $result;
    }

    public function vif()
    {
        $vif = array();
        foreach($this->factors as $index => $factor)
        {
            $arr = $this->factors;
            unset($arr[$index]);
            $this->factorX = $arr;
            $this->factorY = $factor;

            $rsq = $this->vifFactor();
            $vif[] = is_numeric($rsq) ? 1 / (1 - $rsq) : '-';
        }

        return $vif;
    }

    public function vifFactor()
    {
        $X = RegressionMatrix::designMatrix($this->factorX);
        $Xt = $X->transpose(); // Xᵀ
        $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        $Y = RegressionMatrix::yMatrix($this->factorY);
        $coef = $this->coef($X, $Xt, $XtXinv, $Y);
        $B = RegressionMatrix::coefMatrix($coef);

        $ssError = $this->vifSSError($X, $Xt, $XtXinv, $Y, $B);
        $ssTotal = $this->vifSSTotal($Y);

        if(!is_numeric($ssError) or !is_numeric($ssTotal) or $ssTotal === 0) return '-';

        return 1 - $ssError / $ssTotal;
    }

    public function vifSSRegr($X = null, $Xt = null, $XtXinv = null, $Y = null, $B = null)
    {
        if(empty($X)) $X = RegressionMatrix::designMatrix($this->factorX);
        if(empty($Xt)) $Xt = $X->transpose(); // Xᵀ
        if(empty($XtXinv)) $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        if(empty($Y)) $Y = RegressionMatrix::yMatrix($this->factorY);
        if(empty($B))
        {
            $coef = $this->coef($X, $Xt, $XtXinv, $Y);
            $B = RegressionMatrix::coefMatrix($coef);
        }

        $J = RegressionMatrix::oneMatrix($this->observeN, $this->observeN);

        $part1 = $B->transpose()->multiply($Xt)->multiply($Y);
        $part2 = $Y->transpose()->multiply($J)->multiply($Y)->scalarMultiply(1 /  $this->observeN);

        $ss = $part1->subtract($part2);

        return $ss->get(0, 0);
    }

    public function vifSSError($X = null, $Xt = null, $XtXinv = null, $Y = null, $B = null)
    {
        if(empty($X)) $X = RegressionMatrix::designMatrix($this->factorX);
        if(empty($Xt)) $Xt = $X->transpose(); // Xᵀ
        if(empty($XtXinv)) $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        if(empty($Y)) $Y = RegressionMatrix::yMatrix($this->factorY);
        if(empty($B))
        {
            $coef = $this->coef($X, $Xt, $XtXinv, $Y);
            $B = RegressionMatrix::coefMatrix($coef);
        }

        $part1 = $Y->transpose()->multiply($Y);
        $part2 = $B->transpose()->multiply($Xt)->multiply($Y);

        $ss = $part1->subtract($part2);

        return $ss->get(0, 0);
    }

    public function vifSSTotal($Y = null)
    {
        if(empty($Y)) $Y = RegressionMatrix::yMatrix($this->factorY);

        $J = RegressionMatrix::oneMatrix($this->observeN, $this->observeN);

        $part1 = $Y->transpose()->multiply($Y);
        $part2 = $Y->transpose()->multiply($J)->multiply($Y)->scalarMultiply(1 /  $this->observeN);

        $ss = $part1->subtract($part2);

        return $ss->get(0, 0);
    }

    public function pFactor($freedom = null, $t = null)
    {
        if(empty($freedom)) $freedom = $this->dfError();
        if(empty($t)) $t = $this->tValue();

        $studentT = new StudentT($freedom);

        $p = array();
        foreach($t as $value)
        {
            if(!is_numeric($value))
            {
                $p[] = '-';
                continue;
            }
            $pi = 2 * (1 - $studentT->cdf($value));
            $p[] = $pi >= 1 ? 2 - $pi : $pi;
        }

        return $p;
    }

    public function tFactor($coef = null, $se = null)
    {
        if(empty($coef)) $coef = $this->coef();
        if(empty($se)) $se = $this->sdError();

        $t = array();
        foreach($coef as $index => $c)
        {
            if(!is_numeric($c) or !is_numeric($se[$index]) or $se[$index] == 0)
            {
                $t[] = '-';
                continue;
            }
            $t[] = $c / $se[$index];
        }

        return $t;
    }

    public function rsq($ssError = null, $ssTotal = null)
    {
        if(empty($ssError)) $ssError = $this->ssError();
        if(empty($ssTotal)) $ssTotal = $this->ssTotal();

        if(!is_numeric($ssError) or !is_numeric($ssTotal) or $ssTotal == 0) return '-';

        return 1 - $ssError / $ssTotal;
    }

    public function rsqAdj($msError = null, $msTotal = null)
    {
        if(empty($msError)) $msError = $this->msError();
        if(empty($msTotal)) $msTotal = $this->msTotal();

        if(!is_numeric($msError) or !is_numeric($msTotal) or $msTotal == 0) return '-';

        return 1 - $msError / $msTotal;
    }

    public function rsqPre()
    {
        $H    = RegressionMatrix::hMatrix($this->factors);
        $e    = $this->residual();
        $y    = $this->response;
        $mean = Describe::mean($y);

        $sum1 = 0;
        $sum2 = 0;
        foreach($y as $index => $value)
        {
            $h = $H[$index][$index];
            if($h == 1) continue;
            $sum1 += pow($e[$index] / (1 - $h), 2);
            $sum2 += pow($value - $mean, 2);
        }

        if($sum2 == 0) return '-';

        $rsqPre = 1 - ($sum1 / $sum2);
        $rsqPre = $rsqPre < 0 ? 0 : $rsqPre;

        return $rsqPre;
    }

    public function residual($yFit = null)
    {
        $e    = array();
        $y    = $this->response;
        if(empty($yFit)) $yFit = $this->fitY();

        foreach($yFit as $index => $value)
        {
           $e[] = $y[$index] - $value;
        }

        return $e;
    }

    public function sameRow()
    {
        $x = $this->factors;
        $y = $this->response;

        $sameY = array();
        foreach($y as $index => $yValue)
        {
            $key = '';
            foreach($x as $xValue)
            {
                $key .= $xValue[$index];
            }
            if(!isset($sameY[$key])) $sameY[$key] = array();
            $sameY[$key][] = $yValue;
        }

        foreach($sameY as $index => $same)
        {
            if(count($same) == 1) unset($sameY[$index]);
        }

        return $sameY;
    }

    public function itemAnova($ssRegr = null, $msError = null, $dfError = null)
    {
        if(empty($ssRegr)) $ssREgr = $this->ssRegr();
        if(empty($msError)) $msError = $this->msError();
        if(empty($dfError)) $dfError = $this->dfError();

        $items = array();
        foreach($this->factors as $index => $factor)
        {
            $arr = $this->factors;
            unset($arr[$index]);
            $this->factorX = $arr;
            $this->factorY = $this->response;

            $dfItem = $this->dfItem();
            $ssItem = $this->ssItem($ssRegr);
            $msItem = $this->msItem($ssItem, $dfItem);
            $fItem  = $this->fItem($msItem, $msError);
            $pItem  = $this->pItem($dfError, $dfItem, $fItem);

            $items[] = array($dfItem, $ssItem, $msItem, $fItem, $pItem);
        }

        return $items;

    }

    public function dfItem()
    {
        return 1;
    }

    public function ssItem($ssRegr = null)
    {
        if(empty($ssRegr)) $ssRegr = $this->ssRegr();

        $X      = RegressionMatrix::designMatrix($this->factorX);
        $Xt     = $X->transpose(); // Xᵀ
        $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        $Y      = RegressionMatrix::yMatrix($this->response);
        $coef   = $this->coef($X, $Xt, $XtXinv, $Y);
        $B      = RegressionMatrix::coefMatrix($coef);

        $J = RegressionMatrix::oneMatrix($this->observeN, $this->observeN);

        $part1 = $B->transpose()->multiply($Xt)->multiply($Y);
        $part2 = $Y->transpose()->multiply($J)->multiply($Y)->scalarMultiply(1 /  $this->observeN);

        $ssOther = $part1->subtract($part2)->get(0, 0);

        if(!is_numeric($ssRegr)) return '-';

        return $ssRegr - $ssOther;

    }

    public function msItem($ssItem, $dfItem)
    {
        if(!is_numeric($ssItem) or !is_numeric($dfItem) or $dfItem == 0) return '-';
        return $ssItem / $dfItem;
    }

    public function fItem($msItem, $msError)
    {
        if(!is_numeric($msItem) or !is_numeric($msError) or $msError == 0) return '-';
        return $msItem / $msError;
    }

    public function pItem($dfError, $dfItem, $fItem)
    {
        if(!is_numeric($dfError) or !is_numeric($dfItem) or !is_numeric($fItem)) return '-';
        $fdist = new F($dfItem, $dfError);
        return 1 - $fdist->cdf($fItem);
    }

    public function df()
    {
        return $this->coefN;
    }

    public function ssRegr($yFit = null)
    {
        $yMean = Describe::mean($this->response);
        if(empty($yFit)) $yFit = $this->fitY();

        $sum = 0;
        foreach($yFit as $value)
        {
            $sum += ($value - $yMean) ** 2;
        }
        return $sum;
    }

    public function msRegr($ssRegr = null, $df = null)
    {
        if(empty($ssRegr)) $ssRegr = $this->ssRegr();
        if(empty($df)) $df = $this->df();

        if(!is_numeric($ssRegr) or !is_numeric($df) or $df == 0) return '-';

        return $ssRegr / $df;

    }

    public function fRegr($msRegr = null, $msError = null)
    {
        if(empty($msRegr)) $msRegr = $this->msRegr();
        if(empty($msError)) $msError = $this->msError();

        if(!is_numeric($msRegr) or !is_numeric($msError) or $msError == 0) return '-';

        return $msRegr / $msError;
    }

    public function pRegr($df = null, $dfError = null, $fRegr = null)
    {
        if(empty($df)) $df = $this->df();
        if(empty($dfError)) $dfError = $this->dfError();
        if(empty($fRegr)) $fRegr = $this->fRegr();

        if(!is_numeric($dfError) or !is_numeric($df) or !is_numeric($fRegr)) return '-';

        $fdist = new F($df, $dfError);
        return 1 - $fdist->cdf($fRegr);
    }

    public function dfError()
    {
        return $this->observeN - $this->coefN - 1;
    }

    public function ssError($yFit = null)
    {
        $y = $this->response;
        if(empty($yFit)) $yFit = $this->fitY();

        $sum = 0;
        foreach($y as $index => $value)
        {
            $sum += ($value - $yFit[$index]) ** 2;
        }

        return $sum;
    }

    public function msError($ssError = null, $dfError = null)
    {
        if(empty($ssError)) $ssError = $this->ssError();
        if(empty($dfError)) $dfError = $this->dfError();

        if(!is_numeric($ssError) or !is_numeric($dfError) or $dfError === 0) return '-';

        return $ssError / $dfError;
    }

    public function dfTotal()
    {
        return $this->observeN - 1;
    }

    public function ssTotal($vif = false)
    {
        $yMean = Describe::mean($this->response);

        $sum = 0;
        foreach($this->response as $value)
        {
            $sum += ($value - $yMean) ** 2;
        }

        return $sum;
    }

    public function msTotal($ssTotal = null, $dfTotal = null)
    {
        if(empty($ssTotal)) $ssTotal = $this->ssTotal();
        if(empty($dfTotal)) $dfTotal = $this->dfTotal();

        if(!is_numeric($ssTotal) or !is_numeric($dfTotal) or $dfTotal === 0) return '-';

        return $ssTotal / $dfTotal;
    }

    public function dfLOF($dfError = null, $dfPE = null)
    {
        if(empty($dfError)) $dfError = $this->dfError();
        if(empty($dfPE)) $dfPE = $this->dfPE();

        if(!is_numeric($dfError) or !is_numeric($dfPE)) return '-';

        return $dfError - $dfPE;
    }

    public function ssLOF($ssError = null, $ssPE = null)
    {
        if(empty($ssError)) $ssError = $this->ssError();
        if(empty($ssPE)) $ssPE = $this->ssPE();

        if(!is_numeric($ssError) or !is_numeric($ssPE)) return '-';

        return $ssError - $ssPE;
    }

    public function msLOF($ssLOF = null, $dfLOF = null)
    {
        if(empty($ssLOF)) $ssLOF = $this->ssLOF();
        if(empty($dfLOF)) $dfLOF = $this->dfLOF();

        if(!is_numeric($ssLOF) or !is_numeric($dfLOF) or $dfLOF === 0) return '-';

        return $ssLOF / $dfLOF;
    }

    public function fLOF($msLOF = null, $msPE = null)
    {
        if(empty($msLOF)) $msLOF = $this->msLOF();
        if(empty($msPE)) $msPE = $this->msPE();

        if(!is_numeric($msLOF) or !is_numeric($msPE) or $msPE === 0) return '-';

        return $msLOF / $msPE;
    }

    public function pLOF($dfLOF = null, $dfPE = null, $fLOF = null)
    {
        if(empty($dfLOF)) $dfLOF = $this->dfLOF();
        if(empty($dfPE)) $dfPE = $this->dfPE();
        if(empty($fLOF)) $fLOF = $this->fLOF();

        if(!is_numeric($dfLOF) or !is_numeric($dfPE) or !is_numeric($fLOF)) return '-';

        $fdist = new F($dfLOF, $dfPE);
        return 1 - $fdist->cdf($fLOF);
    }

    public function dfPE($same = null)
    {
        if(empty($same)) $same = $this->sameRow();

        $k = count($same);
        $s = 0;
        foreach($same as $values)
        {
            $s += count($values);
        }
        $m = $this->observeN - $s + $k;

        return $this->observeN - $m;
    }

    public function ssPE($same = null)
    {
        if(empty($same)) $same = $this->sameRow();
        $sum = 0;
        foreach($same as $values)
        {
            $mean = Describe::mean($values);
            foreach($values as $value)
            {
                $sum += ($value - $mean) ** 2;
            }
        }

        return $sum;
    }

    public function msPE($ssPE = null, $dfPE = null)
    {
        if(empty($ssPE)) $ssPE = $this->ssPE();
        if(empty($dfPE)) $dfPE = $this->dfPE();

        if(!is_numeric($ssPE) or !is_numeric($dfPE) or $dfPE === 0) return '-';

        return $ssPE / $dfPE;
    }

    public function fitY($coef = null)
    {
        if(empty($coef)) $coef = $this->coef();
        $common = array_fill(0, $this->observeN, 1);
        $x = $this->factors;
        array_unshift($x, $common);

        $y = array();
        for($row = 0; $row < $this->observeN; $row++)
        {
            $sum = 0;
            foreach($x as $index => $value)
            {
                $sum += $value[$row] * $coef[$index];
            }
            $y[] = $sum;
        }

        return $y;
    }

    public function sdError($s² = null, $X = null, $Xt = null, $XtXinv = null)
    {
        return $this->multipleSdError($s², $X, $Xt, $XtXinv);
    }

    public function multipleSdError($s² = null, $X = null, $Xt = null, $XtXinv = null)
    {
        if(empty($s²)) $s² = $this->msError();
        if(empty($X)) $X = RegressionMatrix::designMatrix($this->factors);
        if(empty($Xt)) $Xt = $X->transpose(); // Xᵀ
        if(empty($XtXinv)) $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹

        if(!is_numeric($s²)) return array_fill(0, $this->coefN + 1, '-');

        $sde = $XtXinv->scalarMultiply($s²);

        $diagonal = $sde->getDiagonalElements();

        $sqrtDiago = array();
        foreach($diagonal as $value)
        {
            $sqrtDiago[] = sqrt($value);
        }

        return $sqrtDiago;
    }

    public function coef($X = null, $Xt = null, $XtXinv = null, $Y = null)
    {
        if($this->isSimple) return $this->simpleCoef();

        return $this->multipleCoef($X, $Xt, $XtXinv, $Y);
    }

    public function multipleCoef($X = null, $Xt = null, $XtXinv = null, $Y = null)
    {
        if(empty($X)) $X = RegressionMatrix::designMatrix($this->factors);
        if(empty($Xt)) $Xt = $X->transpose(); // Xᵀ
        if(empty($XtXinv)) $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        if(empty($Y)) $Y = RegressionMatrix::yMatrix($this->response);

        return $XtXinv->multiply($Xt)->multiply($Y)->getColumn(0);
    }

    public function simpleCoef()
    {
        $x     = $this->factor;
        $y     = $this->response;
        $xMean = Describe::mean($x);
        $yMean = Describe::mean($y);

        $sum1 = 0;
        $sum2 = 0;
        foreach($x as $index => $xValue)
        {
            $yValue = $y[$index];
            $sum1 += ($xValue - $xMean) * ($yValue - $yMean);
            $sum2 += ($xValue - $xMean) ** 2;
        }

        $b₁ = $sum1 / $sum2;
        $b₀ = $yMean - $b₁ * $xMean;

        return array($b₀, $b₁);
    }
}
