1;

function [x,y]=bg_gauss(a,b)
  x=linspace(-10,10,21);
  y=b+a*exp(-x.*x/16);
endfunction

function sample=sample_rho(x,y,n)
  sample=0*y;
  integ=cumsum(y);
  minteg=integ(end);
  for i=1:n
    ar=unifrnd(0,minteg);
    idx=lookup(integ,ar)+1;
    sample(idx)=sample(idx)+1;
  endfor
endfunction

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

function l=log_bayes(a,b,x,n,sample)
  dk=dk_distr(a,b,x,n);
  l=sum(sample.*log(dk) - dk);
  l=0;
  for ind=1:numel(x)
    l=l+sample(ind)*log(dk(ind))-dk(ind);
  endfor
endfunction

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

