Statistic Distributions - Printable Version +- 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"); 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"; 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=a AND n<=b, 1/(b-a), 0); return f; END; EXPORT uniformd_cdf(a,b,n) BEGIN local f:= piecewise(n0,{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 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 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;```