Annotation of OpenXM_contrib/gnuplot/demo/stat.inc, Revision 1.1
1.1 ! maekawa 1: #
! 2: # $Id: stat.inc,v 1.3 1998/04/14 00:16:54 drd Exp $
! 3: #
! 4: # Library of Statistical Functions version 3.0
! 5: #
! 6: # Permission granted to distribute freely for non-commercial purposes only
! 7: #
! 8: # Copyright (c) 1991, 1992 Jos van der Woude, jvdwoude@hut.nl
! 9:
! 10: # If you don't have the gamma() and/or lgamma() functions in your library,
! 11: # you can use the following recursive definitions. They are correct for all
! 12: # values i / 2 with i = 1, 2, 3, ... This is sufficient for most statistical
! 13: # needs.
! 14: #logsqrtpi = log(sqrt(pi))
! 15: #lgamma(x) = (x<=0.5)?logsqrtpi:((x==1)?0:log(x-1)+lgamma(x-1))
! 16: #gamma(x) = exp(lgamma(x))
! 17:
! 18: # If you have the lgamma() function compiled into gnuplot, you can use
! 19: # alternate definitions for some PDFs. For larger arguments this will result
! 20: # in more efficient evalution. Just uncomment the definitions containing the
! 21: # string `lgamma', while at the same time commenting out the originals.
! 22: # NOTE: In these cases the recursive definition for lgamma() is NOT sufficient!
! 23:
! 24: # Some PDFs have alternate definitions of a recursive nature. I suppose these
! 25: # are not really very efficient, but I find them aesthetically pleasing to the
! 26: # brain.
! 27:
! 28: # Define useful constants
! 29: fourinvsqrtpi=4.0/sqrt(pi)
! 30: invpi=1.0/pi
! 31: invsqrt2pi=1.0/sqrt(2.0*pi)
! 32: log2=log(2.0)
! 33: sqrt2=sqrt(2.0)
! 34: sqrt2invpi=sqrt(2.0/pi)
! 35: twopi=2.0*pi
! 36:
! 37: # define variables plus default values used as parameters in PDFs
! 38: # some are integers, others MUST be reals
! 39: a=1.0
! 40: alpha=0.5
! 41: b=2.0
! 42: df1=1
! 43: df2=1
! 44: g=1.0
! 45: lambda=1.0
! 46: m=0.0
! 47: mm=1
! 48: mu=0.0
! 49: nn=2
! 50: n=2
! 51: p=0.5
! 52: q=0.5
! 53: r=1
! 54: rho=1.0
! 55: sigma=1.0
! 56: u=1.0
! 57:
! 58: #
! 59: #define 1.0/Beta function
! 60: #
! 61: Binv(p,q)=exp(lgamma(p+q)-lgamma(p)-lgamma(q))
! 62:
! 63: #
! 64: #define Probability Density Functions (PDFs)
! 65: #
! 66:
! 67: # NOTE:
! 68: # The discrete PDFs are calulated for all real values, using the int()
! 69: # function to truncate to integers. This is a monumental waste of processing
! 70: # power, but I see no other easy solution. If anyone has any smart ideas
! 71: # about this, I would like to know. Setting the sample size to a larger value
! 72: # makes the discrete PDFs look better, but takes even more time.
! 73:
! 74: # Arcsin PDF
! 75: arcsin(x)=invpi/sqrt(r*r-x*x)
! 76:
! 77: # Beta PDF
! 78: beta(x)=Binv(p,q)*x**(p-1.0)*(1.0-x)**(q-1.0)
! 79:
! 80: # Binomial PDF
! 81: #binom(x)=n!/(n-int(x))!/int(x)!*p**int(x)*(1.0-p)**(n-int(x))
! 82:
! 83: bin_s(x)=n!/(n-int(x))!/int(x)!*p**int(x)*(1.0-p)**(n-int(x))
! 84: bin_l(x)=exp(lgamma(n+1)-lgamma(n-int(x)+1)-lgamma(int(x)+1)\
! 85: +int(x)*log(p)+(n-int(x))*log(1.0-p))
! 86: binom(x)=(n<20)?bin_s(x):bin_l(x)
! 87:
! 88: # Cauchy PDF
! 89: cauchy(x)=b/(pi*(b*b+(x-a)**2))
! 90:
! 91: # Chi-square PDF
! 92: #chi(x)=x**(0.5*df1-1.0)*exp(-0.5*x)/gamma(0.5*df1)/2**(0.5*df1)
! 93: chi(x)=exp((0.5*df1-1.0)*log(x)-0.5*x-lgamma(0.5*df1)-df1*0.5*log2)
! 94:
! 95: # Erlang PDF
! 96: erlang(x)=lambda**n/(n-1)!*x**(n-1)*exp(-lambda*x)
! 97:
! 98: # Extreme (Gumbel extreme value) PDF
! 99: extreme(x)=alpha*(exp(-alpha*(x-u)-exp(-alpha*(x-u))))
! 100:
! 101: # F PDF
! 102: f(x)=Binv(0.5*df1,0.5*df2)*(df1/df2)**(0.5*df1)*x**(0.5*df1-1.0)/\
! 103: (1.0+df1/df2*x)**(0.5*(df1+df2))
! 104:
! 105: # Gamma PDF
! 106: #g(x)=lambda**rho*x**(rho-1.0)*exp(-lambda*x)/gamma(rho)
! 107: g(x)=exp(rho*log(lambda)+(rho-1.0)*log(x)-lgamma(rho)-lambda*x)
! 108:
! 109: # Geometric PDF
! 110: #geometric(x)=p*(1.0-p)**int(x)
! 111: geometric(x)=exp(log(p)+int(x)*log(1.0-p))
! 112:
! 113: # Half normal PDF
! 114: halfnormal(x)=sqrt2invpi/sigma*exp(-0.5*(x/sigma)**2)
! 115:
! 116: # Hypergeometric PDF
! 117: hypgeo(x)=(int(x)>mm||int(x)<mm+n-nn)?0:\
! 118: mm!/(mm-int(x))!/int(x)!*(nn-mm)!/(n-int(x))!/(nn-mm-n+int(x))!*(nn-n)!*n!/nn!
! 119:
! 120: # Laplace PDF
! 121: laplace(x)=0.5/b*exp(-abs(x-a)/b)
! 122:
! 123: # Logistic PDF
! 124: logistic(x)=lambda*exp(-lambda*(x-a))/(1.0+exp(-lambda*(x-a)))**2
! 125:
! 126: # Lognormal PDF
! 127: lognormal(x)=invsqrt2pi/sigma/x*exp(-0.5*((log(x)-mu)/sigma)**2)
! 128:
! 129: # Maxwell PDF
! 130: maxwell(x)=fourinvsqrtpi*a**3*x*x*exp(-a*a*x*x)
! 131:
! 132: # Negative binomial PDF
! 133: #negbin(x)=(r+int(x)-1)!/int(x)!/(r-1)!*p**r*(1.0-p)**int(x)
! 134: negbin(x)=exp(lgamma(r+int(x))-lgamma(r)-lgamma(int(x)+1)+\
! 135: r*log(p)+int(x)*log(1.0-p))
! 136:
! 137: # Negative exponential PDF
! 138: nexp(x)=lambda*exp(-lambda*x)
! 139:
! 140: # Normal PDF
! 141: normal(x)=invsqrt2pi/sigma*exp(-0.5*((x-mu)/sigma)**2)
! 142:
! 143: # Pareto PDF
! 144: pareto(x)=x<a?0:b/x*(a/x)**b
! 145:
! 146: # Poisson PDF
! 147: poisson(x)=mu**int(x)/int(x)!*exp(-mu)
! 148: #poisson(x)=exp(int(x)*log(mu)-lgamma(int(x)+1)-mu)
! 149: #poisson(x)=(x<1)?exp(-mu):mu/int(x)*poisson(x-1)
! 150: #lpoisson(x)=(x<1)?-mu:log(mu)-log(int(x))+lpoisson(x-1)
! 151:
! 152: # Rayleigh PDF
! 153: rayleigh(x)=lambda*2.0*x*exp(-lambda*x*x)
! 154:
! 155: # Sine PDF
! 156: sine(x)=2.0/a*sin(n*pi*x/a)**2
! 157:
! 158: # t (Student's t) PDF
! 159: t(x)=Binv(0.5*df1,0.5)/sqrt(df1)*(1.0+(x*x)/df1)**(-0.5*(df1+1.0))
! 160:
! 161: # Triangular PDF
! 162: triangular(x)=1.0/g-abs(x-m)/(g*g)
! 163:
! 164: # Uniform PDF
! 165: uniform(x)=1.0/(b-a)
! 166:
! 167: # Weibull PDF
! 168: weibull(x)=lambda*n*x**(n-1)*exp(-lambda*x**n)
! 169:
! 170: #
! 171: #define Cumulative Distribution Functions (CDFs)
! 172: #
! 173:
! 174: # Arcsin CDF
! 175: carcsin(x)=0.5+invpi*asin(x/r)
! 176:
! 177: # incomplete Beta CDF
! 178: cbeta(x)=ibeta(p,q,x)
! 179:
! 180: # Binomial CDF
! 181: #cbinom(x)=(x<1)?binom(0):binom(x)+cbinom(x-1)
! 182: cbinom(x)=ibeta(n-x,x+1.0,1.0-p)
! 183:
! 184: # Cauchy CDF
! 185: ccauchy(x)=0.5+invpi*atan((x-a)/b)
! 186:
! 187: # Chi-square CDF
! 188: cchi(x)=igamma(0.5*df1,0.5*x)
! 189:
! 190: # Erlang CDF
! 191: # approximation, using first three terms of expansion
! 192: cerlang(x)=1.0-exp(-lambda*x)*(1.0+lambda*x+0.5*(lambda*x)**2)
! 193:
! 194: # Extreme (Gumbel extreme value) CDF
! 195: cextreme(x)=exp(-exp(-alpha*(x-u)))
! 196:
! 197: # F CDF
! 198: cf(x)=1.0-ibeta(0.5*df2,0.5*df1,df2/(df2+df1*x))
! 199:
! 200: # incomplete Gamma CDF
! 201: cgamma(x)=igamma(rho,x)
! 202:
! 203: # Geometric CDF
! 204: cgeometric(x)=(x<1)?geometric(0):geometric(x)+cgeometric(x-1)
! 205:
! 206: # Half normal CDF
! 207: chalfnormal(x)=erf(x/sigma/sqrt2)
! 208:
! 209: # Hypergeometric CDF
! 210: chypgeo(x)=(x<1)?hypgeo(0):hypgeo(x)+chypgeo(x-1)
! 211:
! 212: # Laplace CDF
! 213: claplace(x)=(x<a)?0.5*exp((x-a)/b):1.0-0.5*exp(-(x-a)/b)
! 214:
! 215: # Logistic CDF
! 216: clogistic(x)=1.0/(1.0+exp(-lambda*(x-a)))
! 217:
! 218: # Lognormal CDF
! 219: clognormal(x)=cnormal(log(x))
! 220:
! 221: # Maxwell CDF
! 222: cmaxwell(x)=igamma(1.5,a*a*x*x)
! 223:
! 224: # Negative binomial CDF
! 225: cnegbin(x)=(x<1)?negbin(0):negbin(x)+cnegbin(x-1)
! 226:
! 227: # Negative exponential CDF
! 228: cnexp(x)=1.0-exp(-lambda*x)
! 229:
! 230: # Normal CDF
! 231: cnormal(x)=0.5+0.5*erf((x-mu)/sigma/sqrt2)
! 232: #cnormal(x)=0.5+((x>mu)?0.5:-0.5)*igamma(0.5,0.5*((x-mu)/sigma)**2)
! 233:
! 234: # Pareto CDF
! 235: cpareto(x)=x<a?0:1.0-(a/x)**b
! 236:
! 237: # Poisson CDF
! 238: #cpoisson(x)=(x<1)?poisson(0):poisson(x)+cpoisson(x-1)
! 239: cpoisson(x)=1.0-igamma(x+1.0,mu)
! 240:
! 241: # Rayleigh CDF
! 242: crayleigh(x)=1.0-exp(-lambda*x*x)
! 243:
! 244: # Sine CDF
! 245: csine(x)=x/a-sin(n*twopi*x/a)/(n*twopi)
! 246:
! 247: # t (Student's t) CDF
! 248: ct(x)=(x<0.0)?0.5*ibeta(0.5*df1,0.5,df1/(df1+x*x)):\
! 249: 1.0-0.5*ibeta(0.5*df1,0.5,df1/(df1+x*x))
! 250:
! 251: # Triangular PDF
! 252: ctriangular(x)=0.5+(x-m)/g-(x-m)*abs(x-m)/(2.0*g*g)
! 253:
! 254: # Uniform CDF
! 255: cuniform(x)=(x-a)/(b-a)
! 256:
! 257: # Weibull CDF
! 258: cweibull(x)=1.0-exp(-lambda*x**n)
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to stat.inc CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gnuplot / demo