function p=dpo(m)
% p=dpo(m) genererar sann.fkn for Po(m)

N=round(max(20,m+10*sqrt(m))); % övre gräns
z=log(m./(1:N));
lp=cumsum([-m z]); % ger logaritmerade sannolikheter
p=exp(lp);

----------------------

function y=randbin(N,n,p)
% y=randbin(N,n,p) ger N st Bin(n,p) slumptal

U=rand(n,N); % matris av slumptal
V=(U<p); % 0-1 matris
y=sum(V);

----------------------

function y=randdiskr1(p) 
% y=randdiskr(p) ger ett slumptal från sannol.fkn p 

P=cumsum(p); % fordelningsfunktionen 
x=rand(1); 
y=min(find(x<P))-1;

----------------------

function y=randdiskr2(n,p)
% y=randdiskr2(n,p) ger n slumptal från sann.fkn p

p=p(:)'; % omvandlar p till radvektor
P=cumsum(p); % fördelningsfunktionen
x=rand(1,n);
[xs Is]=sort(x);
z=[xs P];
[z I]=sort(z);
[z I]=sort(I);
z=I(1:n)-(1:n); % ger sorterade observationer
y=z(Is); % ger slumpmassig sortering

----------------------

function y=randexp(n,m)
% y=randexp(n,m) ger n Exp(m)-slumptal

u=rand(1,n);
y=-m*log(u);

----------------------

function y=randffg1(p)
% y=randffg1(p) ger ett ffg-fördelat slumptal

y=1;
x=rand(1);
while(x>p)
  x=rand(1);
  y=y+1;
end

----------------------

function y=randffg2(n,p)
% y=randffg2(n,p) ger n ffg(p)-slumptal

m=-1/log(1-p);
x=randexp(n,m);
y=ceil(x);

----------------------

function y=randgamma(N,n,a)
% y=randgamma(N,n,a) ger N gamma(n,a) slumptal

x=-log(rand(n,N))/a; % ger nxN matris av exponentialfördelade slumptal
y=sum(x);

----------------------

function y=randhyp(N,n,p)
% y=randhyp(N,n,p) ger ett slumptal fran Hyp(N,n,p)

v=round(N*p); % v skall vara heltal
x=rand(1,N);
[xs I]=sort(x);
y=sum(I(1:n)<=v);

----------------------

function y=randnorm1(m,s)
% y=randnorm1(m,s) ger 2 N(m,s) slumptal

u=rand(1,2);
z(1)=cos(2*pi*u(1))*sqrt(-2*log(u(2)));
z(2)=sin(2*pi*u(1))*sqrt(-2*log(u(2)));
y=m+s*z;

----------------------

function y=randnorm2(m,s)
% y=randnorm2(m,s) ger 2 N(m,s) slumptal

U=2*rand(1)-1;
V=2*rand(1)-1;
while U^2+V^2>1
  U=2*rand(1)-1;
  V=2*rand(1)-1;
end
R=sqrt(-2*log(rand(1)));
x=[U V]*R/sqrt(U^2+V^2);
y=m+s*x;

----------------------

function y=randnorm3(n,m,s)
% y=randnorm3(n,m,s) ger n N(m,s) slumptal

u=rand(12,n);
x=sum(u)-6;
y=m+s*x;

----------------------

function y=randpoiss(m)
% y=randpoiss(m) ger ett Po(m) slumptal

x=randexp(1,1);
y=0;
while x< m
  x=x+randexp(1,1);
  y=y+1;
end

----------------------

function y=randr(n,a,b)
% randr(n,a,b) ger n R(a,b) slumptal

y=(b-a)*rand(1,n)+a;

----------------------

function y=randurval(n,N)
% y=randurval(n,N) ger n slumpmässiga tal ur 1,2,...,N

x=rand(1,N);
[xs I]=sort(x);
y=I(1:n);

----------------------

function y=randwei(n,a,c)
% y=randwei(n,a,c) ger n weibull-slumptal, skalpar a, formpar c

u=rand(1,n);
y=a*(-log(u)).^(1/c);