lib/MLIB/Graph/external/jpgraph/src/jpgraph_regstat.php

   1 <?php
   2 /*=======================================================================
   3 // File:        JPGRAPH_REGSTAT.PHP
   4 // Description: Regression and statistical analysis helper classes
   5 // Created:     2002-12-01
   6 // Author:      Johan Persson (johanp@aditus.nu)
   7 // Ver:         $Id: jpgraph_regstat.php,v 1.2 2003/03/08 11:29:21 aditus Exp $
   8 //
   9 // License:     This code is released under QPL
  10 // Copyright (C) 2002 Johan Persson
  11 //========================================================================
  12 */
  13
  14 //------------------------------------------------------------------------
  15 // CLASS Spline
  16 // Create a new data array from an existing data array but with more points.
  17 // The new points are interpolated using a cubic spline algorithm
  18 //------------------------------------------------------------------------
  19 class Spline {
  20     // 3:rd degree polynom approximation
  21
  22     var $xdata,$ydata;   // Data vectors
  23     var $y2;             // 2:nd derivate of ydata
  24     var $n=0;
  25
  26     function Spline($xdata,$ydata) {
  27         $this->y2 = array();
  28         $this->xdata = $xdata;
  29         $this->ydata = $ydata;
  30
  31         $n = count($ydata);
  32         $this->n = $n;
  33
  34         // Natural spline 2:derivate == 0 at endpoints
  35         $this->y2[0]    = 0.0;
  36         $this->y2[$n-1] = 0.0;
  37         $delta[0] = 0.0;
  38
  39         // Calculate 2:nd derivate
  40         for($i=1; $i < $n-1; ++$i) {
  41             $d = ($xdata[$i+1]-$xdata[$i-1]);
  42             if( $d == 0  ) {
  43                 JpGraphError::Raise('Invalid input data for spline. Two or more consecutive input X-values are equal. Each input X-value must differ since from a mathematical point of view it must be a one-to-one mapping, i.e. each X-value must correspond to exactly one Y-value.');
  44             }
  45             $s = ($xdata[$i]-$xdata[$i-1])/$d;
  46             $p = $s*$this->y2[$i-1]+2.0;
  47             $this->y2[$i] = ($s-1.0)/$p;
  48             $delta[$i] = ($ydata[$i+1]-$ydata[$i])/($xdata[$i+1]-$xdata[$i]) -
  49                          ($ydata[$i]-$ydata[$i-1])/($xdata[$i]-$xdata[$i-1]);
  50             $delta[$i] = (6.0*$delta[$i]/($xdata[$i+1]-$xdata[$i-1])-$s*$delta[$i-1])/$p;
  51         }
  52
  53         // Backward substitution
  54         for( $j=$n-2; $j >= 0; --$j ) {
  55             $this->y2[$j] = $this->y2[$j]*$this->y2[$j+1] + $delta[$j];
  56         }
  57     }
  58
  59     // Return the two new data vectors
  60     function Get($num=50) {
  61         $n = $this->n ;
  62         $step = ($this->xdata[$n-1]-$this->xdata[0]) / ($num-1);
  63         $xnew=array();
  64         $ynew=array();
  65         $xnew[0] = $this->xdata[0];
  66         $ynew[0] = $this->ydata[0];
  67         for( $j=1; $j < $num; ++$j ) {
  68             $xnew[$j] = $xnew[0]+$j*$step;
  69             $ynew[$j] = $this->Interpolate($xnew[$j]);
  70         }
  71         return array($xnew,$ynew);
  72     }
  73
  74     // Return a single interpolated Y-value from an x value
  75     function Interpolate($xpoint) {
  76
  77         $max = $this->n-1;
  78         $min = 0;
  79
  80         // Binary search to find interval
  81         while( $max-$min > 1 ) {
  82             $k = ($max+$min) / 2;
  83             if( $this->xdata[$k] > $xpoint )
  84                 $max=$k;
  85             else
  86                 $min=$k;
  87         }
  88
  89         // Each interval is interpolated by a 3:degree polynom function
  90         $h = $this->xdata[$max]-$this->xdata[$min];
  91
  92         if( $h == 0  ) {
  93             JpGraphError::Raise('Invalid input data for spline. Two or more consecutive input X-values are equal. Each input X-value must differ since from a mathematical point of view it must be a one-to-one mapping, i.e. each X-value must correspond to exactly one Y-value.');
  94         }
  95
  96
  97         $a = ($this->xdata[$max]-$xpoint)/$h;
  98         $b = ($xpoint-$this->xdata[$min])/$h;
  99         return $a*$this->ydata[$min]+$b*$this->ydata[$max]+
 100              (($a*$a*$a-$a)*$this->y2[$min]+($b*$b*$b-$b)*$this->y2[$max])*($h*$h)/6.0;
 101     }
 102 }
 103
 104 // EOF
 105 ?>