# _____________________________________________________________
# Kommenttirivi: Tama on octave ohjelmani Opns2.m
# - Pienimman neliosumman sovitus
# _____________________________________________________________
graphics_toolkit('gnuplot')    # Vahentaa grafiikka ongelmia
clear ; clc                    # Poistetaan ... Tyhjennetaan ...
# _____________________________________________________________
function y=ffun(x,beta)  # leasqr tarvitsee taman aliohjelman 
y=beta(1)*cos(x)+beta(2)*sin(x);
endfunction
# _____________________________________________________________
pkg load optim          ; # Added spring 2017
# _______________________________________________________
# - Luodaan n0 simuloitua datan arvoa x0 ja y0 lisaamalla 
# malliin g0 = a0 cos x + b0 sin x satunnaisvirhe e0. 
# Taman normaalijakautuneen virheen keskihajonta on s0. 
# Simuloitu data on y0 = g0 + e0.
# _____________________________________________________________
# - Loopissa n0 ja s0 arvot MUUTTUU!!!!!!!!!!!!!!!!!!!!!!!!!!!!
n =[10 10 50 50] ; s=[0.1 1.0 1.0 5.0];
# - Loopissa subplot komennon z1 ja z2 arvot MUUTTUU!!!!!!!!!!!
z1=[ 3  7 35 39] ; z2=[17 21 49 53];
# _____________________________________________________________
for i=1:4 # LOOPPI ALKAA >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# _____________________________________________________________
n0=n(i);                               # n0 MUUTTUU!!!!!!!!!!!!
x0=2.*pi*((0:n0-1)/(n0-1));            # 0 < x < 2 pi
s0=s(i) ; e0=s0*normrnd(0.,1.,[1 n0]); # e0 MUUTTUU!!!!!!!!!!!!
a0=1.  ; b0=2.;                        # a0 ja b0 
g0=a0*cos(x0)+b0*sin(x0);              # g0
y0=g0+e0;                              # y0
# _____________________________________________________________
# Pienimman neliosumman sovitus
beta0=[1 1];                           # arvattu beta alkuarvo
[yfit beta1 cvg iter] = leasqr(x0, y0, beta0, "ffun");
# yfit=sovitus arvot, cvg=1 (sovitus OK), iter=iterointi maara
# _____________________________________________________________
xx=2.*pi*(0:100)/100.;                    # Sopiva xx
yy=beta1(1)*cos(xx)+beta1(2)*sin(xx);     # Sovitettu yy
gg=     a0*cos(xx)+     b0*sin(xx);       # Simuloitu gg
# _____________________________________________________________
subplot(8,8,[z1(i),z2(i)])            # MUUTTUU!!!!!!!!!!!!!!!!
x1=0.0            ; x2=2.0*pi;        # x-rajat
y1=min(g0)-3.*s0  ; y2=max(g0)+3.*s0; # y-rajat
axis([x1,x2,y1,y2]);                  # xy-rajat
set(gca,'fontsize',5);                # xy-akseli numerojen koko
hold on;                              # plot ei poista edellista
plot(x0,g0,'ob','markersize',3) ; plot(xx,gg,'b') # plotit
plot(x0,y0,'or','markersize',3) ; plot(xx,yy,'r') # plotit
hold off;                             # uusi plot poistaisi ...
# _____________________________________________________________
txt1=sprintf("%3s%3.0f%6s%5.2f",'n=',n0,' ja s0=',s0);
text(x1+0.0*(x2-x1),y1+1.20*(y2-y1),txt1,'fontsize',8);   
txt2=sprintf("%13s%5.2f%6s%5.2f",'Simuloitu A=',a0,' ja B=',b0);
text(x1+0.0*(x2-x1),y1+1.12*(y2-y1),txt2,'fontsize',8);   
txt3=sprintf("%13s%5.2f%6s%5.2f",'Sovitettu A=',beta1(1),' ja B=',beta1(2));
text(x1+0.0*(x2-x1),y1+1.04*(y2-y1),txt3,'fontsize',8);   
# _____________________________________________________________
endfor # LOOPPI LOPPUU >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# _____________________________________________________________
print('Opns2.jpg','-djpg')          # Kuva tiedosto
waitforbuttonpress;                 # Odota: Nappain poistaa
 
