+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP Prime Software Library (/forum-15.html)
+--- Thread: Statistic Distributions (/thread-3850.html)
Statistic Distributions - salvomic - 05-15-2015 07:48 PM
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é
RE: Statistic Distributions - salvomic - 12-23-2015 10:18 PM
Gumbel (and Standard Gumbel) distribution
Code:
EXPORT gumbel(m,b, k)
// Gumbel distribution m=μ location, b=β>0 scale, k=x var
// also log-Weibull
BEGIN
local f, z;
IF b<=0 THEN RETURN("Use: m, b>0, x"); END;
z:= (k-m)/b;
f:= (1/b)*e^-(z+e^-z);
RETURN f;
END;
EXPORT gumbel_cdf(m,b,k)
BEGIN
local f;
IF b<=0 THEN RETURN("Use: m, b>0, x"); END;
f:= e^-(e^(-(k-m)/b));
END;
EXPORT gumbel_std(k)
// Standard Gumbel distribution
// (m=μ location, b=β>0 scale, k=x var)
// μ=0 and β=1
BEGIN
local f;
f:= e^-(k+e^-k);
RETURN f;
END;
EXPORT gumbel_std_cdf(k)
BEGIN
local f;
f:= e^-e^(-k);
END;Fréchet Distribution (also inverse Weibull)
Code:
EXPORT frechet(a,s,m,k)
// Fréchet distribution
// also inverse Weibull
// a=α>0 shape, s>0 scale (default 1)
// m localtion of minimum (default 0)
// k=x var
BEGIN
local f;
IF a<=0 OR s<=0 THEN RETURN("Use: a>0, s>0, m[def 0], x"); END;
f:= (a/s)*((k-m)/s)^(-1-a)*e^(-((k-m)/s)^(-a));
RETURN f;
END;
EXPORT frechet_cdf(a,s,m,k)
BEGIN
local f;
IF a<=0 OR s<=0 THEN RETURN("Use: a>0, s>0, m[def 0], x"); END;
f:= e^(-((k-m)/s)^(-a));
RETURN f;
END;