Imprimir

# Calculator blog

Musings and comments about our common interest

## Final post on statistical distributions on HP Prime

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;

Leer mensaje completo
Más sobre:

## Statistical distributions for the HP Prime - 2

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

// 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;

BEGIN

local f;

f:= int(X(a-1)*e(-X),X,0,l*n)/Gamma(a);

return f;

END ;

// 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;

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;

Leer mensaje completo
Más sobre:

## Statistical Distributions in HP Prime

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:

• Normal,
• T (Student),
• Chi-square,
• F (Fisher distribution),
• Binomial and
• Poisson.

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;

Leer mensaje completo