Distributions

[salvomic]

hi all,
in the Prime we have density (and also cumulative and inverse) commands for some distribution: Normal, T (Student), Chi-square, F (Fisher distribution), Binomial and Poisson. Than we have Zeta (for zipf d.)

Sometime (for didactical purpose) I need to treat also other distributions: exponentiald, geometric, negbinomial, uniformd, hypergeometric.

Any hints to implement any of them in a simple way, in a formula definition or program, almost for density and cumulative type?
Especially geometric, hypergeometric and negbinomial.

Thank you!

***
EDIT: I would suggest to insert some of them others in the Prime...

[walter b]

I'd suggest you take a look at the WP 34S Owner's Manual. You find it here: http://sourceforge.net/projects/wp34s/fi...f/download - please turn to App. I therein.

d:-)

[parisse]

BTW, there are available in Xcas.

[salvomic]

(02-19-2015 06:18 AM)walter b Wrote:  I'd suggest you take a look at the WP 34S Owner's Manual. You find it here: http://sourceforge.net/projects/wp34s/fi...f/download - please turn to App. I therein.
d:-)

(02-19-2015 07:35 AM)parisse Wrote:  BTW, there are available in Xcas.

thanks both!

WP 34S is a wonderful calc, and Xcas it's always ok!

However I would like have those distributions also in the Prime

[Offroad]

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

[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

[salvomic]

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 ;

[salvomic]

...the most simple (apparently):

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;

[salvomic]

Another, that requires a bit of attention...
see here for the definition.
See also this thread for the multinomial coefficient.

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;

[salvomic]

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;

[salvomic]

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;

[salvomic]

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;

[DrD]

Nice job! it's been interesting following your progress with the Prime!

From my perspective, its also refreshing to see someone accomplishing productive results, and sharing them in this forum. Sometimes, after wading through all the complaints these days, I'm left to wonder if the grass roots "can do" attitude has vaporized from society. Perhaps its a throwback to my amateur radio heritage, but I admire those who can accomplish useful things with what they have, in spite of idiosyncrasies within their kit of resources.

I don't know your specific background, but superficially, it seems you've exemplified a good example of how learning both the tool, and the source material, can come together to reach something useful and, in which, others might benefit!

-Dale-

[salvomic]

(04-08-2015 01:22 PM)DrD Wrote:  Nice job! it's been interesting following your progress with the Prime!

hi Dale, thank you!

Quote:From my perspective, its also refreshing to see someone accomplishing productive results, and sharing them in this forum. Sometimes, after wading through all the complaints these days, I'm left to wonder if the grass roots "can do" attitude has vaporized from society. Perhaps its a throwback to my amateur radio heritage, but I admire those who can accomplish useful things with what they have, in spite of idiosyncrasies within their kit of resources.

yes, I like sharing in the forum what I learn about the tool, and I hope everybody here do the same. As you know, I'm ham too, and my experience in radio communication has thought to me this behavior...
I hope these pieces of programs (commands) could help someone to expand its knowledge about Maths and Statistics. I would like some of these (with other already present in XCas) can be included in the Prime, if they are seen as "important"

Quote:I don't know your specific background, but superficially, it seems you've exemplified a good example of how learning both the tool, and the source material, can come together to reach something useful and, in which, others might benefit!
-Dale-

I'm learning (again) with the Prime
Really, I'm trying to "squeeze" the Prime as an orange, and I can see that it can do really more than it already do and this for me is very amusing and amazing...
My job is indeed Informatics and Electronics, and my passions are Maths, Statistics, Radio and every form of measure's instruments and calculation machines.

Salvo

[Paul Dale]

Look in the WP34 directory: trunk/xrom/distributions.

The implementations for the ones you want that the 34S supports are there as keystroke programs. They are all pretty straightforward.


Pauli

[Claudio L.]

(04-08-2015 01:22 PM)DrD Wrote:  ... I'm left to wonder if the grass roots "can do" attitude has vaporized from society. Perhaps its a throwback to my amateur radio heritage, but I admire those who can accomplish useful things with what they have, in spite of idiosyncrasies within their kit of resources.

It's a generational change. The new generation had internet since their infancy, therefore they are naturally a "search and download" generation.
The previous generation had personal computers but relatively isolated (no www yet), so they were natural "programmers".
The one before them had only basic calculators and kits to build their own computers, so they were natural "electronics designers, builders, and tinkerers in general".

It's clear to see, even in this forum: some people are designing calculators, soldering iron in hand (older generation). Some people are writing software for their calculators, and some people just ask where to download something, and complain if it doesn't exist. If you ask for the age of each person, you'll see a perfect correlation.

[salvomic]

(04-08-2015 10:04 PM)Claudio L. Wrote:  It's clear to see, even in this forum: some people are designing calculators, soldering iron in hand (older generation). Some people are writing software for their calculators, and some people just ask where to download something, and complain if it doesn't exist. If you ask for the age of each person, you'll see a perfect correlation.

indeed

[Gerald H]

(04-08-2015 01:22 PM)DrD Wrote:  Nice job! it's been interesting following your progress with the Prime!

From my perspective, its also refreshing to see someone accomplishing productive results, and sharing them in this forum. Sometimes, after wading through all the complaints these days, I'm left to wonder if the grass roots "can do" attitude has vaporized from society. Perhaps its a throwback to my amateur radio heritage, but I admire those who can accomplish useful things with what they have, in spite of idiosyncrasies within their kit of resources.

I don't know your specific background, but superficially, it seems you've exemplified a good example of how learning both the tool, and the source material, can come together to reach something useful and, in which, others might benefit!

-Dale-

Yes, it's good to see more programmes being published.

I agree that working within the confines of limited resources sharpens the wits.

I don't agree that the Prime is limited resources, in fact I'd describe it as overpowered combined with arbitrarily handicapped, eg largest integers that can be factorized.

[debrouxl]

In absolute terms, all graphing calculators on the marketplace are laughably underpowered, and sold at very high price tags for such limited CPU power, amount of RAM, amount of Flash memory, limited I/O possibilities, low programmability, etc.
In relative terms, the Prime is one of the two most powerful models of calculators on the marketplace. Compared to the Nspire CX (CAS), the Prime has twice the amount of Flash, a CPU about twice faster, but sadly, only half the RAM.