1;

function [x,y]=bg_gauss(a,b)
  ## probability that a single event is triggered in a small unit of time
  ## choose small values for a,b
  x=linspace(-10,10,21);
  y=b+a*exp(-x.*x/16);
endfunction

function sample=sample_rho(x,y,t)
  sample=0*y;
  for isample=1:numel(x)
    for js=1:t
      ar=unifrnd(0,1);
      if (ar<=y(isample))
        sample(isample)+=1;
      endif
    endfor
  endfor
endfunction

function dk=dk_distr(a,b,x,t)
  dk=b+a*exp(-x.*x/16);
  dk=dk*t;
endfunction

function l=log_bayes(a,b,t,sample)
  sample=sample/t;
  [~,y]=bg_gauss(a,b);
  l=sum(sample.*log(y) - y);
endfunction

function lv=log_bayesv(a,b,t,sample)
  for ind=1:numel(a)
    for jnd=1:numel(b)
      lv(ind,jnd)=log_bayes(a(ind),b(jnd),t,sample);
 #     [ind,jnd]
    endfor
  endfor
endfunction

function res=coin_bias(n_heads,n_tails,p_head)
  res=p_head.^n_heads.*(1-p_head).^n_tails;
endfunction
