vcs.maemo.org Git - gnuplot/blob - demo/stat.inc

   1 #
   2 # $Id: stat.inc,v 1.3 2006/06/14 03:24:09 sfeam Exp $
   3 #
   4 # Library of Statistical Functions version 3.0
   5 #
   6 # Copyright (c) 1991, 1992 Jos van der Woude, jvdwoude@hut.nl
   7
   8 # History:
   9 #    -- --- 1992 Jos van der Woude:  1st version
  10 #    06 Jun 2006 Dan Sebald:  Defined PDF/CDF for whole real line or integer
  11 #                             set and range checked all other parameters.
  12
  13 #
  14 # Shortcut for testing if a variable is an integer
  15 #
  16 isint(x)=(int(x)==x)
  17
  18 # Define useful constants
  19 fourinvsqrtpi=4.0/sqrt(pi)
  20 invpi=1.0/pi
  21 invsqrt2pi=1.0/sqrt(2.0*pi)
  22 log2=log(2.0)
  23 sqrt2=sqrt(2.0)
  24 sqrt2invpi=sqrt(2.0/pi)
  25 twopi=2.0*pi
  26
  27 #
  28 #define 1.0/Beta function
  29 #
  30 Binv(p,q)=exp(lgamma(p+q)-lgamma(p)-lgamma(q))
  31
  32 # NOTE:
  33 #
  34 # The stat functions are defined appropriately for the whole real line or
  35 # set of integers (the first input variable).  There are restrictions on
  36 # some of the parameters, and an undefined value will result if the input
  37 # value falls outside the parameter's range.
  38 #
  39 # The discrete PDFs (and some parameters) must have integer (natural number)
  40 # inputs, otherwise an undefined result is produced.  This means the user
  41 # must appropriately supply a discrete input set, perhaps by some form of
  42 # scaling before and after calling the stat desired function.  To plot the
  43 # output of such a discrete data set passed through a stat function, the
  44 # user can make use of plotting features such as "with steps", and so on.
  45
  46 #
  47 #define Probability Density Functions (PDFs)
  48 #
  49
  50 # Arcsin PDF
  51 arcsin(x,r)=r<=0?1/0:abs(x)>r?0.0:invpi/sqrt(r*r-x*x)
  52
  53 # Beta PDF
  54 beta(x,p,q)=p<=0||q<=0?1/0:x<0||x>1?0.0:Binv(p,q)*x**(p-1.0)*(1.0-x)**(q-1.0)
  55
  56 # Binomial PDF
  57 binom(x,n,p)=p<0.0||p>1.0||n<0||!isint(n)?1/0:\
  58   !isint(x)?1/0:x<0||x>n?0.0:exp(lgamma(n+1)-lgamma(n-x+1)-lgamma(x+1)\
  59   +x*log(p)+(n-x)*log(1.0-p))
  60
  61 # Cauchy PDF
  62 # a location parameter, b > 0 scale parameter
  63 cauchy(x,a,b)=b<=0?1/0:b/(pi*(b*b+(x-a)**2))
  64
  65 # Chi-square PDF
  66 chisq(x,k)=k<=0||!isint(k)?1/0:\
  67   x<=0?0.0:exp((0.5*k-1.0)*log(x)-0.5*x-lgamma(0.5*k)-k*0.5*log2)
  68
  69 # Erlang PDF
  70 erlang(x,n,lambda)=n<=0||!isint(n)||lambda<=0?1/0:\
  71   x<0?0.0:x==0?(n==1?real(lambda):0.0):exp(n*log(lambda)+(n-1.0)*log(x)-lgamma(n)-lambda*x)
  72
  73 # Extreme (Gumbel extreme value) PDF
  74 extreme(x,mu,alpha)=alpha<=0?1/0:alpha*(exp(-alpha*(x-mu)-exp(-alpha*(x-mu))))
  75
  76 # F PDF
  77 f(x,d1,d2)=d1<=0||!isint(d1)||d2<=0||!isint(d2)?1/0:\
  78   Binv(0.5*d1,0.5*d2)*(real(d1)/d2)**(0.5*d1)*x**(0.5*d1-1.0)/(1.0+(real(d1)/d2)*x)**(0.5*(d1+d2))
  79
  80 # Gamma PDF
  81 # rho > 0 shape parameter, lambda > 0 inverse scale parameter
  82 gmm(x,rho,lambda)=rho<=0||lambda<=0?1/0:\
  83   x<0?0.0:x==0?(rho>1?0.0:rho==1?real(lambda):1/0):\
  84   exp(rho*log(lambda)+(rho-1.0)*log(x)-lgamma(rho)-lambda*x)
  85
  86 # Geometric PDF
  87 # p probability of success, x number of failures before first success
  88 geometric(x,p)=p<=0||p>1?1/0:\
  89   !isint(x)?1/0:x<0||p==1?(x==0?1.0:0.0):exp(log(p)+x*log(1.0-p))
  90
  91 # Half normal PDF
  92 halfnormal(x,sigma)=sigma<=0?1/0:x<0?0.0:sqrt2invpi/sigma*exp(-0.5*(x/sigma)**2)
  93
  94 # Hypergeometric PDF
  95 # N objects, C of one class (N-C of another), d drawn without
  96 # replacement, x drawn of class C.
  97 hypgeo(x,N,C,d)=N<0||!isint(N)||C<0||C>N||!isint(C)||d<0||d>N||!isint(d)?1/0:\
  98   !isint(x)?1/0:x>d||x>C||x<0||x<d-(N-C)?0.0:exp(lgamma(C+1)-lgamma(C-x+1)-lgamma(x+1)+\
  99   lgamma(N-C+1)-lgamma(d-x+1)-lgamma(N-C-d+x+1)+lgamma(N-d+1)+lgamma(d+1)-lgamma(N+1))
 100
 101 # Laplace PDF
 102 laplace(x,mu,b)=b<=0?1/0:0.5/b*exp(-abs(x-mu)/b)
 103
 104 # Logistic PDF
 105 logistic(x,a,lambda)=lambda<=0?1/0:lambda*exp(-lambda*(x-a))/(1.0+exp(-lambda*(x-a)))**2
 106
 107 # Lognormal PDF
 108 lognormal(x,mu,sigma)=sigma<=0?1/0:\
 109   x<0?0.0:invsqrt2pi/sigma/x*exp(-0.5*((log(x)-mu)/sigma)**2)
 110
 111 # Maxwell PDF
 112 # a sqrt(2) times standard deviation of individual component of normal triple
 113 maxwell(x,a)=a<=0?1/0:x<0?0.0:fourinvsqrtpi*a**3*x*x*exp(-a*a*x*x)
 114
 115 # Negative binomial PDF
 116 # p probability of success, r number of success to complete,
 117 # x failures before success r
 118 negbin(x,r,p)=r<=0||!isint(r)||p<=0||p>1?1/0:\
 119   !isint(x)?1/0:x<0?0.0:p==1?(x==0?1.0:0.0):exp(lgamma(r+x)-lgamma(r)-lgamma(x+1)+\
 120   r*log(p)+x*log(1.0-p))
 121
 122 # Negative exponential PDF
 123 nexp(x,lambda)=lambda<=0?1/0:x<0?0.0:lambda*exp(-lambda*x)
 124
 125 # Normal PDF
 126 normal(x,mu,sigma)=sigma<=0?1/0:invsqrt2pi/sigma*exp(-0.5*((x-mu)/sigma)**2)
 127
 128 # Pareto PDF
 129 pareto(x,a,b)=a<=0||b<0||!isint(b)?1/0:x<a?0:real(b)/x*(real(a)/x)**b
 130
 131 # Poisson PDF
 132 poisson(x,mu)=mu<=0?1/0:!isint(x)?1/0:x<0?0.0:exp(x*log(mu)-lgamma(x+1)-mu)
 133
 134 # Rayleigh PDF
 135 rayleigh(x,lambda)=lambda<=0?1/0:x<0?0.0:lambda*2.0*x*exp(-lambda*x*x)
 136
 137 # Sine PDF
 138 # f frequency, a length
 139 sine(x,f,a)=a<=0?1/0:\
 140   x<0||x>=a?0.0:f==0?1.0/a:2.0/a*sin(f*pi*x/a)**2/(1-sin(twopi*f))
 141
 142 # t (Student's t) PDF
 143 t(x,nu)=nu<0||!isint(nu)?1/0:\
 144   Binv(0.5*nu,0.5)/sqrt(nu)*(1+real(x*x)/nu)**(-0.5*(nu+1.0))
 145
 146 # Triangular PDF
 147 triangular(x,m,g)=g<=0?1/0:x<m-g||x>=m+g?0.0:1.0/g-abs(x-m)/real(g*g)
 148
 149 # Uniform PDF
 150 uniform(x,a,b)=x<(a<b?a:b)||x>=(a>b?a:b)?0.0:1.0/abs(b-a)
 151
 152 # Weibull PDF
 153 weibull(x,a,lambda)=a<=0||lambda<=0?1/0:\
 154   x<0?0.0:x==0?(a>1?0.0:a==1?real(lambda):1/0):\
 155   exp(log(a)+a*log(lambda)+(a-1)*log(x)-(lambda*x)**a)
 156
 157 #
 158 #define Cumulative Distribution Functions (CDFs)
 159 #
 160
 161 # Arcsin CDF
 162 carcsin(x,r)=r<=0?1/0:x<-r?0.0:x>r?1.0:0.5+invpi*asin(x/r)
 163
 164 # incomplete Beta CDF
 165 cbeta(x,p,q)=ibeta(p,q,x)
 166
 167 # Binomial CDF
 168 cbinom(x,n,p)=p<0.0||p>1.0||n<0||!isint(n)?1/0:\
 169   !isint(x)?1/0:x<0?0.0:x>=n?1.0:ibeta(n-x,x+1.0,1.0-p)
 170
 171 # Cauchy CDF
 172 # a location parameter, b > 0 scale parameter
 173 ccauchy(x,a,b)=b<=0?1/0:0.5+invpi*atan((x-a)/b)
 174
 175 # Chi-square CDF
 176 cchisq(x,k)=k<=0||!isint(k)?1/0:x<0?0.0:igamma(0.5*k,0.5*x)
 177
 178 # Erlang CDF
 179 cerlang(x,n,lambda)=n<=0||!isint(n)||lambda<=0?1/0:x<0?0.0:igamma(n,lambda*x)
 180
 181 # Extreme (Gumbel extreme value) CDF
 182 cextreme(x,mu,alpha)=alpha<=0?1/0:exp(-exp(-alpha*(x-mu)))
 183
 184 # F CDF
 185 cf(x,d1,d2)=d1<=0||!isint(d1)||d2<=0||!isint(d2)?1/0:1.0-ibeta(0.5*d2,0.5*d1,d2/(d2+d1*x))
 186
 187 # incomplete Gamma CDF
 188 # rho > 0 shape parameter, lambda > 0 inverse scale parameter
 189 cgmm(x,rho,lambda)=rho<=0||lambda<=0?1/0:x<0?0.0:igamma(rho,x*lambda)
 190
 191 # Geometric CDF
 192 # p probability of success, x number of failures before first success
 193 cgeometric(x,p)=p<=0||p>1?1/0:\
 194   !isint(x)?1/0:x<0||p==0?0.0:p==1?1.0:1.0-exp((x+1)*log(1.0-p))
 195
 196 # Half normal CDF
 197 chalfnormal(x,sigma)=sigma<=0?1/0:x<0?0.0:erf(x/sigma/sqrt2)
 198
 199 # Hypergeometric CDF
 200 # N objects, C of one class (N-C of another), d drawn without
 201 # replacement, x drawn of class C.
 202 chypgeo(x,N,C,d)=N<0||!isint(N)||C<0||C>N||!isint(C)||d<0||d>N||!isint(d)?1/0:\
 203   !isint(x)?1/0:x<0||x<d-(N-C)?0.0:x>d||x>C?1.0:hypgeo(x,N,C,d)+chypgeo(x-1,N,C,d)
 204
 205 # Laplace CDF
 206 claplace(x,mu,b)=b<=0?1/0:(x<mu)?0.5*exp((x-mu)/b):1.0-0.5*exp(-(x-mu)/b)
 207
 208 # Logistic CDF
 209 clogistic(x,a,lambda)=lambda<=0?1/0:1.0/(1+exp(-lambda*(x-a)))
 210
 211 # Lognormal CDF
 212 clognormal(x,mu,sigma)=sigma<=0?1/0:x<=0?0.0:cnormal(log(x),mu,sigma)
 213
 214 # Maxwell CDF
 215 # a sqrt(2) times standard deviation of individual component of normal triple
 216 cmaxwell(x,a)=a<=0?1/0:x<0?0.0:igamma(1.5,a*a*x*x)
 217
 218 # Negative binomial CDF
 219 # p probability of success, r number of success to complete,
 220 # x failures before success r
 221 cnegbin(x,r,p)=r<=0||!isint(r)||p<=0||p>1?1/0:\
 222   !isint(x)?1/0:x<0?0.0:ibeta(r,x+1,p)
 223
 224 # Negative exponential CDF
 225 cnexp(x,lambda)=lambda<=0?1/0:x<0?0.0:1-exp(-lambda*x)
 226
 227 # Normal CDF
 228 cnormal(x,mu,sigma)=sigma<=0?1/0:0.5+0.5*erf((x-mu)/sigma/sqrt2)
 229
 230 # Pareto CDF
 231 cpareto(x,a,b)=a<=0||b<0||!isint(b)?1/0:x<a?0.0:1.0-(real(a)/x)**b
 232
 233 # Poisson CDF
 234 cpoisson(x,mu)=mu<=0?1/0:!isint(x)?1/0:x<0?0.0:1-igamma(x+1.0,mu)
 235
 236 # Rayleigh CDF
 237 crayleigh(x,lambda)=lambda<=0?1/0:x<0?0.0:1.0-exp(-lambda*x*x)
 238
 239 # Sine CDF
 240 # f frequency, a length
 241 csine(x,f,a)=a<=0?1/0:\
 242   x<0?0.0:x>a?1.0:f==0?real(x)/a:(real(x)/a-sin(f*twopi*x/a)/(f*twopi))/(1.0-sin(twopi*f)/(twopi*f))
 243
 244 # t (Student's t) CDF
 245 ct(x,nu)=nu<0||!isint(nu)?1/0:0.5+0.5*sgn(x)*(1-ibeta(0.5*nu,0.5,nu/(nu+x*x)))
 246
 247 # Triangular PDF
 248 ctriangular(x,m,g)=g<=0?1/0:\
 249   x<m-g?0.0:x>=m+g?1.0:0.5+real(x-m)/g-(x-m)*abs(x-m)/(2.0*g*g)
 250
 251 # Uniform CDF
 252 cuniform(x,a,b)=x<(a<b?a:b)?0.0:x>=(a>b?a:b)?1.0:real(x-a)/(b-a)
 253
 254 # Weibull CDF
 255 cweibull(x,a,lambda)=a<=0||lambda<=0?1/0:x<0?0.0:1.0-exp(-(lambda*x)**a)
 256