Musings and comments about our common interest
Let’s us finish with the distribution functions for the HP Prime. Remember that we posted two blog issues about it: here and here.
All were taken from user Salvomic in HPmuseum. This is a monumental job and a great help for those that cannot bear to have an approximation when you can have the real thing.
Here are the last ones:
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/(s2))*e(-(t2)/(2*(s2)))));
return f;
END;
EXPORT rayleigh_cdf(s, t)
BEGIN
local f;
f:= piecewise(t<0, 0, 1 - e(-t2/(2*(s2))));
return f;
END;
Pareto distribution
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*xma)/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*La*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-La*k(-a))/(1-(L/H)^a);
return f;
END;
n-Erlang distribution
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, ((lk)*(n(k–1)e^(-ln)))/(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))/((mk)*(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;
We continue with the distributions as promised on last blog issue:
Beta Distribution:
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 distribution
EXPORT gammad(a,l,n)
// Gamma distributio 1st formn a=α>0 shape param
// l=λ>0 rate param, n var
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
BEGIN
local f;
f:=piecewise(n<0,0, (n(k-1)*e(-n/t)) / ((tk) * 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 distribution
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/(Xs), X, 1, k);
// nth generalized armonic number
f:= hks/Zeta(s);
return f;
END;
Laplace distribution:
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
BEGIN
local f;
f:= (1/(2b))e^(-(ABS(n-m))/b);
return f;
END;
EXPORT laplaced_cdf(m,b,n)
BEGIN
local f;
f := piecewise(n<n, (1/2)*e((n-m)/b), 1-(1/2)*e(-(n-m)/b));
return f;
END ;
Uniform Distribution
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;
start: postbit_signature
Multinomial Distribution:
The multinomial distribution is a generalization of the binomial distribution. For n independent trials each of which leads to a success for exactly one of k categories, with each category having a given fixed success probability, the multinomial distribution gives the probability of any particular combination of numbers of successes for the various categories.
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(pk);
END; // if
END;
There is an interesting article about statistical distributions for the HP Prime in hpmuseum.org. As you may know, the HP Prime comes with several functions to provide statistical distributions. The ones that come to my mind now are:
These come in all three flavors: density, cumulative and inverse. So, compared with my beloved HP41CL, it is very well provided. But there are many other distributions that may be needed if you’re heavily into statistics. User Salvomic provides programs for several interesting cases. As he does not claim any copyright to them, and in the interest of the readers, I will post here the programs that create them. For the sake of reader comfort, we’ll post them in several installments, lest the page becomes too long!
Logistic Distribution:
EXPORT logisticd(m,s,k)
// Logistic distribution m=µ, s=s, k=x
BEGIN
local f,g;
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;
g:= (k-m)/s;
f := 1/(1+e^(-g));
return f;
END ;
EXPORT logisticd_icdf(m,s,p)
// inverse m=µ, s=s, p probability
BEGIN
local f;
f:= m+s*ln(p/(1-p));
return f;
END;
Lognormal
EXPORT LgNrm(m,s,k)
// LogNormal distribution m=µ, s=s, k=x
BEGIN
local f;
f:= (1/(kssqrt(2pi)))e(-(ln(k)-m)2/(2*s2));
return f;
END;
EXPORT LgNrm_cdf(m,s,k)
BEGIN
local f;
f := normal((ln(k)-m)/s) ;
return f;
END ;
Exponential
EXPORT expond(l,n)
// exponential distribution l=λ=1/np, n
BEGIN
local f;
f:= piecewise(n<0,0, le^(-ln));
return f;
END;
EXPORT expond_cdf(l,n)
BEGIN
local f;
f := piecewise(n<0, 0, 1-e^(-l*n));
return f;
END ;
Geometric
EXPORT geometric(n,p)
// Geometric distribution, n tries, p prob
BEGIN
local f;
f:= p*(1-p)^n;
return f;
END;
EXPORT geometric_cdf(n,p)
BEGIN
local f;
f := 1-(1-p)^(n+1);
return f;
END ;
Hypergeometric
EXPORT ipergeom(n,m,k,j)
// Hypergeometric distribution, n total, m one kind, k extracted, j var
BEGIN
local f;
f:= (comb(m,j)*comb(n-m,k-j))/(comb(n,k));
return f;
END;
Negative Binomial
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,…
BEGIN
local f;
IF (n<r) THEN return “n must be >= r”; ELSE
f:= comb(n–1, r–1)*(pr)*(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))*(pk)*((1-p)r);
return f;
END;
Gompertz
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:= bhe(b*n)*ehe(-h*e(bn)) ;
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
EXPORT weibull(l, k, t)
// Weibull distribution: l=λ>0 scale param (characteristic lifetime), k>0 shape parameter t var (time)
BEGIN
local f;
f:= piecewise(t<0,0, (k/l)(t/l)^(k–1)e(-(t/l)k));
return f;
END;
EXPORT weibull_cdf(l, k, t)
BEGIN
local f;
f:= piecewise(t<=0,0, 1-e(-(t/l)k)) ;
return f;
END ;
EXPORT weibull_translate(l, k, t, v)
// Weibull distribution tranlated: v location parameter
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(l, k, t, v)
BEGIN
local f;
f:= piecewise(t<=v,0, 1-e(-((t-v)/l)k)) ;
return f;
END ;
Cauchy
EXPORT cauchyd(x0,g,n)
// Cauchy distribution x0 location param, g=γ scale param, n var
BEGIN
local f;
f:= 1/ (pig(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+gtan(pi(p-(1/2)));
return f;
END;