with(plots):
with(stats):
with(describe):
with(random):

average:= (t) -> mean(t):
N:=(x,m,sig)-> statevalf[cdf, normald[m,sqrt(sig)]](x):
B:=(n,p,k)->(n!*(p^k)*(1-p)^(n-k))/(k!*(n-k)!):
P:=(m,k)->(m^k*exp(-m))/k!:
y:=(a,b)->1/(b-a):
UNI:=(a,b)->plot([[a-sqrt(b^2),0], [a,0], [a,y(a,b)], [b,y(a,b)], [b,0], [b+sqrt(a^2),0]], style=LINE):
g:=(x,m,v)->(exp(-(x-m)^2/(2*v)))/sqrt(2*Pi*v):
phi:=(z)->exp(-(z^2)/2)/sqrt(2*Pi):

MINLS:=proc(ls)
    local i,mini;
    i:='i': mini:=infinity:
    for i from 1 to nops(ls) do mini:=min(mini,ls[i]): od:
    mini;
end:

MAXLS:=proc(ls)  
    local i, maxi;
    i:='i': maxi:=-infinity:
    for i from 1 to nops(ls) do maxi:=max(maxi,ls[i]): od:
    maxi;
end:
                  
FREQHIST:=proc(ls)
    local range,wx,x,m,i,n,w,c,opls,p,q,npls;
    range:=MAXLS(ls)-MINLS(ls):
    wx:=round(range*100):
    x:=array(1..1,0..wx):
    m:='m': i:='i': n:='n':
    for i from 0 to wx do x[1,i]:=0: od:
    for m from 1 to nops(ls) do n:=ls[m]: w:=round(100*(n-MINLS(ls))); c:=x[1,w]: x[1,w]:=c+1: od:
    n:='n': opls:=[]:
    for n from 0 to wx do p:=n/100+MINLS(ls): q:=x[1,n]: npls:=[op(opls),[p,0],[p,q],[p,0]]: opls:=npls: od:
    plot(npls,style=LINE);
end:

#R:= poisson[10.0](generator[15]):
R:=rand(10..20);

SAMPLE := proc (c)
     local j, samp, no;
     j:='j':
     samp:=[]:
     for j from 1 to c do no:=R(): samp:=[op(samp),no]: od;
end:

MEAN:=proc(t)
     local me;
     me:=evalf(average(t));
end:                       

VEST:=proc(u)
     local vest;
     vest:=evalf(variance[0](u));
end:

VARIANCE:=proc(g)
      local n,va;
      n:=nops(g):
      va:=evalf(variance[0](g)*(n-1)/n);
end:

SDEV:=proc(h)         
      local s;
      s:=evalf(sqrt(VARIANCE(h)));
end:

SAMPLEMEANS:=proc (d,f)
     local k, samp, smpm, spm;
     k:='k':
     smpm:=[]:
     for k from 1 to f do samp:=SAMPLE(d): spm:=MEAN(samp): smpm:=[op(smpm),spm]: od;
     smpm;
end:

SAMPLESDEVS:=proc(g,h)
     local k,samp,smsd,spsd;
     k:='k': smsd:=[]:
     for k from 1 to h do samp:=SAMPLE(g): spsd:=SDEV(samp): smsd:=[op(smsd),spsd]: od;
     smsd;
end:

samp1:=SAMPLE(15):

MEANS:=proc(o,p,q,r)
    local ms,i,s;
    ms:=[]:
    i:='i':
    for i from 1 to p do s:=(q-r[i])*o*p: ms:=[op(ms),s]: od;
end:

SAMPMHIST:=proc(c,d,f)
    local msm,means;
    msm:=average(f):
    means:=MEANS(c,d,msm,f):
    HISTGM(means,`Sample means`);
end:

CHISQ:=proc(X,v)
    local ans1,ans2,ans;
    ans1:=int(x^(.5*v-1)*exp(-.5*x),x=X..infinity):
    ans2:=int(x^(.5*v-1)*exp(-.5*x),x=0..infinity):
    ans:=100*ans1/ans2;
end: