%kappamu.m
function [res1,res2,kappa,mu]=kappamu(m,a,b)
s1=0;
s2=0;
for i=1:1:m
s1=s1+(randn(1,1000000)-a).^2;
s2=s2+(randn(1,1000000)-b).^2;
end
s=s1+s2;
[res1 res2]=hist(sqrt(s),100);
kappa=m*(a^2+b^2)/(2*m);
mu=mean(s);
v=var(s);
mu=((1+2*kappa)/(1+kappa)^2)*(mu^2/v);
%plotkappamu.m
[p1,q1,k1,m1]=kappamu(10,1,2);
[p2,q2,k2,m2]=kappamu(1,1,2);
[p3,q3,k3,m3]=kappamu(1,0,0);
[p4,q4,k4,m4]=kappamu(10,0,0);
figure
subplot(2,2,1)
plot(q1,p1)
title(strcat('kappa=',num2str(k1), 'mu=',num2str(m1)))
subplot(2,2,2)
plot(q2,p2,'r')
title(strcat('kappa=',num2str(k2),'mu=',num2str(m2)))
subplot(2,2,3)
plot(q3,p3)
title(strcat('kappa=',num2str(k3),'mu=',num2str(m3)))
subplot(2,2,4)
plot(q4,p4,'r')
title(strcat('kappa=',num2str(k4),'mu=',num2str(m4)))