Statistic Distributions
05-15-2015, 07:48 PM (This post was last modified: 05-15-2015 07:55 PM by salvomic.)
Post: #1
 salvomic Senior Member Posts: 1,366 Joined: Jan 2015
Statistic Distributions
hi all,
with the new firmware 7820 some new statistic's distributions are been added to the Prime, a part of Normal, T (Student), F (Fisher), Chi square (χ2), Binomial, Poisson.
Now also other 8 from XCas: Beta, Cauchy, Exponential, Gamma, Geometric, NegBinomial, Uniform, Weibull. They are in Catalog (as normal form, cdf, icdf).

Other (and some alternative to those already in the Prime) can be found here in my programs.
If some of my functions should have the same name of those in the Prime, you can easily change the function's name in the code...

***
Logistic distribution
Code:
 EXPORT logisticd(m,s,k) // Logistic distribution m=μ location , s=σ >0 scale, k=x var BEGIN local f,g; IF s<=0 THEN RETURN("Use: m, s>0, k"); END; g:= (k-m)/s; f:=e^(-g)/(s*(1+e^(-g))^2) ; return f; END; EXPORT logistic_cdf(m,s,k) BEGIN local f,g; IF s<=0 THEN RETURN("Use: m, s>0, k"); END; g:= (k-m)/s; f := 1/(1+e^(-g)); return f; END ; EXPORT logisticd_icdf(m,s,p) // inverse m=μ, s=σ, p probability BEGIN local f; IF s<=0 THEN RETURN("Use: m, s>0, p"); END; f:= m+s*ln(p/(1-p)); return f; END;

Lognormal
Code:
 EXPORT LgNrm(m,s,k) // LogNormal distribution m=μ, s=σ>0 shape, k=x var BEGIN local f; f:= piecewise(k<=0,0, (1/(k*s*sqrt(2*pi)))*e^(-(ln(k)-m)^2/(2*s^2))); return f; END; EXPORT LgNrm_cdf(m,s,k) // LgNrm(m,s,k) = normal((ln(k)-m)/s) - normal(k) = LgNrm(0,1,e^k) BEGIN local f; f:=(1/2)+(1/2)*erf((ln(k)-m)/(sqrt(2)*s)); return f; END ;

Exponential
Code:
 EXPORT expond(l,n) // exponential distribution l=λ=1/np, n // expond(1/2, n) = chi2(2, n) // expon(λ) = gammad2(1, 1/λ) BEGIN local f; f:= piecewise(n<0,0,  l*e^(-l*n)); return f; END; EXPORT expond_cdf(l,n) BEGIN local f; f := piecewise(n<0, 0, 1-e^(-l*n)); return f; END ;

Geometric
Code:
 EXPORT geometric(p, k) // Geometric distribution, k trial, p probability (kth trial first success) // geometric(p) = Negbinom(1,p) BEGIN local f; IF k<1 THEN RETURN("k must be int >= 1"); END; k:= IP(k); f:= p*(1-p)^(k-1); return f; END; EXPORT geometric2(p, k) // Geometric distribution, k trial, p probability (failure until 1st success) BEGIN local f; IF k<0 THEN RETURN("k must be int >= 0"); END; k:= IP(k); f:= p*(1-p)^(k); return f; END; EXPORT geometric_cdf(p, k) BEGIN local f; IF k<1 THEN RETURN("k must be int >= 1"); END; f := 1-(1-p)^k; return f; END ; EXPORT geometric2_cdf(p, k) BEGIN local f; IF k<0 THEN RETURN("k must be int >= 0"); END; f := 1-(1-p)^(k+1); return f; END ;

Hypergeometric and negative hypergeometric
Code:
 EXPORT hypergeom(N, n, m, k) // Hypergeometric distribution // N population size, m (or K) #successes in population // n number of draws (without replacement), k (or i)  number of successes BEGIN local f; f:= (comb(m,k)*comb(N-m,n-k))/(comb(N,n)); return f; END; EXPORT negHypergeom(r, n, m, k) // Negative Hypergeomtric distribution // r nth special item, N=n+m, n special, m normal, k var BEGIN local f; IF (k<r OR r<0 OR k<0) THEN RETURN("Use: (r, n, m, k) - k must be >= r"); END; f:=(comb(n,r-1)*comb(m,k-r)/comb(n+m,k-1))*((n-r+1)/(n+m-k+1)); return f; END;

Negative Binomial
Code:
 EXPORT NegBin(r, p,n) // Negative Binomial distribution observing until r success, with p probability of success //  n num trials for r success (k failures), n=r, r+1,... // negBinom(1,p) = geometric(p) BEGIN local f; IF (n<r) THEN return "n must be >= r"; ELSE f:= comb(n-1, r-1)*(p^r)*(1-p)^(n-r); return f; END; //if END; EXPORT NegBin_cdf(r, p, n) BEGIN local f, b1, a, b; a:=r; b:= n+1; b1:= int((X^(a-1))*(1-X)^(b-1),X,0,p); // incomplete beta function f:=( b1/Beta(a,b)); return f; END; EXPORT NegBin2(r, p,k) // Negative Binomial observing until r failures, with p probability of success (1-p failure) // n num trials, k = n-r success // i.e. NegBin(5,0.4,15) = NegBin2(5,0.6,10) BEGIN local f; f:= (comb(k+r-1,k))*(p^k)*((1-p)^r); return f; END;

Gompertz
Code:
 EXPORT gompertz(h,b,n) // Gompertz distribution h=η shape param, b scale param, n var BEGIN local f; IF (n<0 OR h<=0 OR b<=0) THEN return "Not defined for η or b < =0 or n <0"; ELSE f:= b*h*e^(b*n)*e^h*e^(-h*e^(b*n)) ; return f; END; // if END; EXPORT gompertz_cdf(h,b,n) BEGIN local f; IF (n<0 OR h<=0 OR b<=0) THEN return "Not defined for η or b < =0 or n <0"; ELSE f := 1-e^(-h*(e^(b*n)-1)); return f; END; END ;

Weibull and Weibull translated

Code:
 EXPORT weibull(k, l, t) // Weibull distribution: l=λ>0 scale param (characteristic lifetime), k>0 shape parameter t var (time) // weibull(λ,1) = exponential(1/λ), weibull(λ,2) = Rayleigh(λ/sqrt(2)) BEGIN local f; f:= piecewise(t<0,0,  (k/l)*(t/l)^(k-1)*e^(-(t/l)^k)); return f; END; EXPORT weibull_cdf(k, l, t) BEGIN local f; f:= piecewise(t<=0,0, 1-e^(-(t/l)^k)) ; return f; END ; EXPORT weibull_translate(k, l, v, t) // Weibull distribution translated: v (or θ theta) location parameter, normally 0 BEGIN local f; f:= piecewise(t<=v,0,  (k/l)*((t-v)/l)^(k-1)*e^(-((t-v)/l)^k)); return f; END; EXPORT weibull_translate_cdf(k, l, v, t) BEGIN local f; f:= piecewise(t<=v,0, 1-e^(-((t-v)/l)^k)) ; return f; END ;

Cauchy
Code:
 EXPORT cauchyd(x0,g,n) // Cauchy distribution x0 location param, g=γ scale param, n var BEGIN local f; f:= 1/ (pi*g*(1+((n-x0)/g)^2)) ; return f; END; EXPORT cauchy_cdf(x0,g,n) BEGIN local f; f:= (1/pi)*atan((n-x0)/g)+1/2 ; return f; END ; EXPORT cauchy_icdf(x0,g,p) // inverse x0, g=γ, p probability BEGIN local f; f:= x0+g*tan(pi*(p-(1/2))); return f; END;

Beta
Code:
 EXPORT betad(a, b, n) // Beta distribution: a=α>0, b=β>0 shape param, n var (0<=n<=1) BEGIN local f; f:= piecewise(n<0 ,0, n>=1, 0, (1/Beta(a,b))*(n^(a-1))*(1-n)^(b-1)); return f; END; EXPORT betad_cdf(a, b, n) BEGIN local f, b1; b1:= int((X^(a-1))*(1-X)^(b-1),X,0,n); // incomplete beta function f:=piecewise(n<0,0,n>=1, 0,  b1/Beta(a,b)); return f; END;

Gamma
Code:
 EXPORT gammad(a,l,n) // Gamma distribution 1st form a=α>0 shape param //  l=λ>0 rate param, n var // gammad::gamma2d α=k, β=1/θ // gamma(1,1/λ) = expon(λ), gamma(n/2,1/2) = chi2(n) BEGIN local f; f:= piecewise( n<0,0, (l*e^(-l*n)*(l*n)^(a-1))/Gamma(a)  ); END; EXPORT gammad_cdf(a,l,n) BEGIN local f; f:= int(X^(a-1)*e^(-X),X,0,l*n)/Gamma(a); return f; END ; EXPORT gammad2(k,t,n) // Gamma distribution 2nd form  k>0 shape param, // t=θ>0 scale param, n var // if k is N (natural) -> Erlang distribution (k=1 -> exponential) BEGIN local f; f:=piecewise(n<0,0,  (n^(k-1)*e^(-n/t)) / ((t^k) * Gamma(k)) ); return f; END; EXPORT gammad2_cdf(k,t,n) BEGIN local f; f:= int(X^(k-1)*e^(-X),X,0,n/t)/Gamma(k); return f; END ;

Zeta (Zpif)
Code:
 EXPORT Zetazipf(s, k) // Zeta (Zipf) distribution, not defined in s=1 BEGIN local f; IF s=1 THEN return "Not defined in s=1";  ELSE f:= (k^(-s))/Zeta(s); return f; END;  //if END; EXPORT Zetazipf_cdf(s,k) BEGIN local f, hks; hks:= sum(1/(X^s), X, 1, k); // nth generalized armonic number f:= hks/Zeta(s); return f; END;

Laplace
Code:
 EXPORT laplaced(m,b,n) // Laplace distribution m=μ location param, b scale param,  n var // if m=0 and b=1 -> expond scaled by 1/2 (λ=1/b) BEGIN local f; f:= (1/(2*b))*e^(-(ABS(n-m))/b); return f; END; EXPORT laplaced_cdf(m,b,n) BEGIN local f; f := piecewise(n<m, (1/2)*e^((n-m)/b), 1-(1/2)*e^(-(n-m)/b)); return f; END ;

Uniform distribution
Code:
 EXPORT uniformd(a,b,n) // Uniform, distribution [a-b], n var BEGIN local f:= piecewise(n>=a AND n<=b, 1/(b-a), 0); return f; END; EXPORT uniformd_cdf(a,b,n) BEGIN local f:= piecewise(n<a, 0, n≥b, 1, (n-a)/(b-a)); return f; END;

Multinomial distribution (requires n = total number of items, list of k items of a kind, list of their probabilities...)

Code:
 EXPORT Multinomd(n, k, p) // Multinomial distribution (n>0,{list k_i}, {list p_i}) BEGIN IF ((type(k) ≠ 6) OR (type(p) ≠ 6)) THEN  return "2nd and 3th argument must be a list"; ELSE IF n<1 THEN return "n must be >0"; END; n:= ip(n); //n must be integer IF (size(k) ≠ size(p) ) THEN return "items in k must be = those in p"; END; IF (ΣList(k) ≠ n) THEN return "ΣList(k) must be = n"; END; IF (ΣList(p) > 1) THEN return "ΣList(p) must be <= 1"; END; // ΣList(p) should be = 1 but we can use in list k only values with no 0 //so items in p could be less than those in k... return  (n!/ΠLIST(k!))*ΠLIST(p^k); END; // if  END;

n-Erlang distribution (see also Gamma and exponential)

Code:
 EXPORT erlang(k,l,n) // n-Erlang distribution k shape parameter, l=λ >=0 rate parameter // from Gamma d.; if k=1 -> erlang(1,l,n) = exponential(l,n) // erlang(k,λ) = gamma(k,1/λ) BEGIN local f; f:=piecewise(n<0,0, l<0, 0, ((l^k)*(n^(k-1)*e^(-l*n)))/(k-1)! ); END; EXPORT erlang2(k, m, n) //k shape parameter,  m=μ=1/λ >=0 scale parameter // if μ=2 -> chi2 with 2k degree of freedom BEGIN local f; f:=piecewise(n<0,0, m<0, 0, (n^(k-1)*e^(-n/m))/((m^k)*(k-1)!)); END; EXPORT erlang_cdf(k, l, n) // k shape parameter, l=λ >=0 rate parameter (μ=1/λ) BEGIN local f; f:= 1- sum((1/X!)*(e^(-l*n))*(l*n)^X,X,0,k-1); END;

Rayleigh distribution

Code:
 EXPORT rayleigh(s, t) // Rayleigh distribution s=σ scale parameter, t>=0 var // rayleigh(2σ^2) = weibull(σ*sqrt(2),2) BEGIN local f; f:= piecewise(t<0, 0, ((t/(s^2))*e^(-(t^2)/(2*(s^2))))); return f; END; EXPORT rayleigh_cdf(s, t) BEGIN local f; f:= piecewise(t<0, 0, 1 - e^(-t^2/(2*(s^2)))); return f; END;

Pareto distribution

Code:
 EXPORT paretod(xm, a, k) // Pareto distribution x_m >0 scale, a=α > 0 shape BEGIN local f; IF ((a<=0) OR (xm <= 0) OR (k<xm)) THEN RETURN("Use: xm > 0, a >0, k >= xm"); END; f:=(a*xm^a)/k^(a+1); return f; END; EXPORT paretod_cdf(xm, a, k) BEGIN local f; IF ((a<=0) OR (xm <= 0) OR (k<xm)) THEN RETURN("Use: xm > 0, a >0, k >= xm"); END; f:=1-(xm/k)^a; return f; END; EXPORT paretod_bound(a, L, H, k) // Bounded Pareto distribution // a=α >0 shape, L>0, H>L location, k var BEGIN local f; IF (a<=0 OR L<=0 OR H<=L) THEN RETURN("Use a>0, L>0, H>L"); END; f:= (a*L^a*k^(-a-1))/(1-(L/H)^a); return f; END; EXPORT paretod_bnd_cdf(a,L,H,k) BEGIN local f; IF (a<=0 OR L<=0 OR H<=L) THEN RETURN("Use a>0, L>0, H>L"); END; f:= (1-L^a*k^(-a))/(1-(L/H)^a); return f; END;

The important and classic Maxwell-Boltzmann Distribution

Code:
 EXPORT maxwell_boltzmannd(a, n) // Maxwell-Boltzmann distribution // scale parameter a = SQRT(kT/m), k=Boltzmann constant // chi2 with 3 DF BEGIN local f; f:=piecewise(a<0, 0, sqrt(2/PI)*((n^2)*e^((-n^2)/(2*a^2)))/(a^3)  ); return f; END; EXPORT maxwell_boltzman_cdf(a, n) // Maxwell-Boltzmann distribution // scale parameter a = SQRT(kT/m), k=Boltzmann constant // chi2 with 3 DF BEGIN local f; f:=piecewise(a<0, 0, erf(n/(sqrt(2)*a) )- (sqrt(2/PI))*(n*e^((-n^2)/(2*a^2)))/a  ); return f; END;

Enjoy with Statistic in the Prime!

Salvo Micciché

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
 « Next Oldest | Next Newest »

 Messages In This Thread Statistic Distributions - salvomic - 05-15-2015 07:48 PM RE: Statistic Distributions - salvomic - 12-23-2015, 10:18 PM

User(s) browsing this thread: 1 Guest(s)