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