Distributions

[salvomic]

(02-19-2015 11:58 AM)Offroad Wrote:  You could make an app yourself in the meantime if you need them that much. It's part of the fun of owing a Prime. At the end you'll have your very own customized high performance calculator that meets all your needs.

Offroad

yes, you're right!
I love my Prime...

Let's go, but, please, help me to find errors, improve the programs and simplify them

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 ;

My purpose is now standardize these programs also using the same variable naming of the others distributions in the Prime and make controls to suggest and guide the user input. Any help much appreciated!

I hope the programs could be helpful to whom use Prime for statistics.

Salvo