<?php
require_once dirname(__DIR__) . '/base/hypothesis.php';
require_once dirname(dirname(__DIR__)) . '/zendasmath/distribution/stdnormal.php';
require_once dirname(dirname(__DIR__)) . '/zendasmath/distribution/beta.php';
require_once dirname(dirname(__DIR__)) . '/zendasmath/distribution/binomial.php';

/**
 * Proportion1
 */
class proportion1 extends hypothesis
{
    public $dataframe;

    /**
     * Confidence interval ci for the normal method
     * https://support.minitab.com/zh-cn/minitab/21/help-and-how-to/statistics/basic-statistics/how-to/1-proportion/methods-and-formulas/methods-and-formulas/#confidence-interval-ci-for-the-normal-approximation
     *
     * @param int $event
     * @param int $sample
     * @param int $cl
     * @param int $limit up low all
     * @access public
     * @return array
     */
    public function normalCI($event, $sample, $cl = 95, $limit = 'all')
    {
        $standardNormal = new standardNormal();
        $a      = 1 - $cl / 100;
        $p      = $limit == 'all' ? 1 - $a / 2 : 1 - $a;
        $z      = $standardNormal->inverse($p);
        $result = array('-', '-');

        $p = $event / $sample;

        if($event !=0 && $event != $sample)
        {
            $part = $z * sqrt($p * (1 - $p) / $sample);
            $result = array();
            $lower  = $p - $part;
            $lower  = max(0, $lower);
            $upper  = $p + $part;
            $upper  = min(1, $upper);
            $result[] = $limit != 'up' ? round($lower, 6) : '-';
            $result[] = $limit != 'low' ? round($upper, 6) : '-';
        }

        return $result;
    }

    /**
     * normal method test method;
     *
     * @param int $event
     * @param int $sample
     * @param float $p₀
     * @param string $type
     * @access public
     * @return float
     */
    public function normalTest($event, $sample, $p₀, $type, $round = true)
    {
        $p = $event / $sample;
        $z = ($p - $p₀) / sqrt($p₀ * (1 - $p₀) / $sample);
        $standardNormal = new standardNormal();
        $cdf           = $standardNormal->cdf($z);
        if($type == 'gt')
        {
            $result = 1 - $cdf;
        }

        if($type == 'lt')
        {
            $result = $cdf;
        }

        if($type == 'noEqual')
        {
            $result = 2 * (1 - $cdf);
            $result = $result > 1 ? 2 - $result : $result;
        }

        $result = $this->dataframe->roundDigit($result, $round);
        $z      = $this->dataframe->roundDigit($z, $round);
        return array($z, $result);
    }

    /**
     * Confidence interval ci for the exact method
     * https://support.minitab.com/zh-cn/minitab/21/help-and-how-to/statistics/basic-statistics/how-to/1-proportion/methods-and-formulas/methods-and-formulas/#confidence-interval-ci-for-the-exact-method
     *
     * @param int $event
     * @param int $sample
     * @param int $limit up low all
     * @access public
     * @return array
     */
    public function exactCI($event, $sample, $cl, $limit = 'all')
    {
        $α      = (1 - $cl / 100) / 2;
        $result = array();

        if($event != 0)
        {
            if($event == $sample or $limit == 'low') $α = 2 * $α;
            if($event == $sample) $result[1] = '1.000000';

            $x = 1 - $α;
            $a = $sample - $event + 1;
            $b = $event;
            $beta = new Beta($a, $b);

            $LL   = 1 - $beta->inverse($x);
            $LL   = $this->dataframe->toString($LL, 6);

            $result[0] = $LL;
        }

        if($event != $sample)
        {
            if($event == 0 or $limit == 'up') $α = 2 * $α;
            if($event == 0) $result[0] = '0.000000';

            $x = $α;
            $a = $sample - $event;
            $b = $event + 1;
            $beta = new Beta($a, $b);

            $UL   = 1 - $beta->inverse($x);
            $UL   = $this->dataframe->toString($UL, 6);

            $result[1] = $UL;
        }

        if($limit == 'low') $result[1] = '-';
        if($limit == 'up') $result[0] = '-';

        return $result;
    }

    /**
     * Exact method test method;
     *
     * @param int $event
     * @param int $sample
     * @param float $p
     * @param string $type
     * @access public
     * @return float
     */
    public function exactTest($event, $sample, $p0, $type, $round = true)
    {
        $bino   = new Binomial($p0, $sample);
        if($type == 'gt')
        {
            $result = 1 - $bino->cdf($event) + $bino->pdf($event);
        }

        if($type == 'lt')
        {
            $result = $bino->cdf($event);
        }

        if($type == 'noEqual')
        {
            $p = $event / $sample;
            if($p0 == 0.5)
            {
                $y       = min($event, $sample - $event);
                $n_y     = $sample - $y;

                $yp     = $bino->cdf($y);
                $n_yCdf = $bino->cdf($n_y);
                $n_yPdf = $bino->pdf($n_y);
                $n_yp   = 1 - $n_yCdf + $n_yPdf;
                $result = $yp + $n_yp;
            }
            else
            {
                $LRyLog = $this->likelihoodLog($p / $p0, (1 - $p) / (1 - $p0), $event, $sample);
                $xRange = array();
                $last   = null;
                $begin  = null;

                $result = '?';
                for($x = 0; $x <= $sample; $x++)
                {
                    $px     = $x / $sample;
                    $LRxLog = $this->likelihoodLog($px / $p0, (1 - $px) / (1 - $p0), $x, $sample);

                    $isGt = ($LRxLog >= $LRyLog || is_nan($LRxLog));

                    if($x === 0 and $isGt) $begin = $x;

                    if($x !== 0)
                    {
                        if($isGt && !$last) $begin = $x;
                        if(!$isGt && $last) $xRange[] = array($begin, $x - 1);
                    }

                    if($x == $sample)
                    {
                        if($isGt)
                        {
                           if(!$last) $begin = $x;
                           $xRange[] = array($begin, $x);
                        }
                        if(!$isGt && $last) $xRange[] = array($begin, $x);
                    }

                    $last = $isGt;
                }

                $result = 0;
                foreach($xRange as $range)
                {
                    if($range[0] == 0)
                    {
                        $result += $bino->cdf($range[1]);
                    }

                    if($range[1] == $sample && $range[0] !== $range[1])
                    {
                        $result += 1 - $bino->cdf($range[0]);
                        $result += $bino->pdf($range[0]);
                    }
                }
            }
        }
        $result = max(0, $result);
        $result = min(1, $result);

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

        return $result;
    }

    public function likelihood($p, $p0, $x, $n)
    {
        $part1 = $p / $p0; // p / p0
        $part2 = (1 - $p) / (1 - $p0); // (1 - p) / (1 - p0)

        return pow($part1, $x) * pow($part2, $n - $x);
    }

    public function likelihoodLog($p_p0, $p1_p01, $x, $n)
    {
        return $x * log($p_p0) + ($n - $x) * log($p1_p01);
    }
}
