option nolet
open #1:screen 0,0.5,0,1
open #2:screen .5,1,.5,1
window #1
! serie :analytical chemistry
! filename :trimodalfitfinalversion1
! rem : seg =3 (tri)modalfit
print "*********************************************"
print "* 'version 1.0',trimodal fit *"
print "* first seg : 1 - > p *"
print "* second seg : p+1 - > q *"
print " third seg : q+1 - > n (n=number of data) *"
print "* use:conductometric tit.,spectrofotometric *"
print "* data ,from lab experiment. *"
print "* algo: least square computation,trimodalfit*"
print "*********************************************"
print "* [email protected],25/10/2016 *"
print "*********************************************"
print
print "************"
print "* given = *"
print "************"
!**********************************
! info to display data on the graph
!**********************************
dim xp(2),yp(2)
tb_x = 0 ! xmin
te_x = 20 ! xmax
nt_x = 10 ! tick marks x-axis
tb_y = 1 ! ymin
te_y = 10 ! ymax
nt_y = 12 ! tick marks y-axis
!*********************************************************
! info : graph,title,numbering,label axes
!*********************************************************
title$ = "graph:'x data versus y data'"
x_txt$ = "x-data"
y_txt$ = "y-data"
txt = 2 ! distance from right part of the window
!read data statement, x data, y data
dim xy(2,1)
!read x data
n=16 ! number x data = number y data
mat redim xy(2,n)
! x data
for i=1 to n
read xy(1,i)
next i
print " number of data = ";n
print
! y data
for i=1 to n
read xy(2,i)
next i
print " x-data = '=> "
print
for i=1 to n
print xy(1,i);
next i
print
print " y-data = ' => "
print
for i=1 to n
print xy(2,i);
next i
print
print
print "*******************************"
print " data range: x,y;tick x,y *"
print "*******************************"
print "xmin = ";tb_x,"xmax = ";te_x," number of tick mark,x = ";nt_x
print "ymin = ";tb_y,"ymax = ";te_y," number of tick mark,y = ";nt_y
print "************"
print "* result = *"
print "************"
!sum x,average x,ii=1;sum y,average y,ii=2
call sum_val_average(1,n,xy(,),somvalx,avgx)
call sum_val_average(2,n,xy(,),somvaly,avgy)
print " sum x_values,sum y_values = ";somvalx;",";somvaly
print " mean x_values,mean y_values = ";avgx;",";avgy
! sum prod of x_values*y_values
call sum_prod_x_y(n,xy(,),somprodxy)
print " sum(prod(x,y)) = ";somprodxy
! sum prod of x_values*x_values
call sum_prod_x_x(1,n,xy(,),somprodxx)
! sum prod of y_values*y_values
call sum_prod_x_x(2,n,xy(,),somprodyy)
print " sum(prod(x,x)) = ";somprodxx," sum(prod(y,y)) = ";somprodyy
! calculation regression line ' y = aa*x+bb'
spxx = somprodxx
spyy = somprodyy
spxy = somprodxy
sx = somvalx
sy = somvaly
call regr_line_aa(n,spxy,spxx,spyy,sx,sy,aa)
call regr_line_bb(n,avgx,avgy,aa,bb)
print
print "regression line,blue => y = ";aa;" * x(+)";bb
call correlation_coef(spxy,spxx,spyy,r_cor)
print
print "correlation coef. = ";r_cor
print
!*******************************************
! algo:
!*******************************************
! find minimum qp : with respect to aa1,bb1
! (1 to p),aa2,bb2 (p+1 to q)
! ,aa3,bb3 (q+1 to n),y1,y2,y3
!*******************************************
dim data_b(1,1)
! pp,qq,aa1,bb1,aa2,bb2,aa3,bb3 => num_col
num = 0
!*******************************************
seg = 3 ! number of segments
!*********************************************
! math :
!*********************************************
!*********************************************
!number of segments =>'n = number of(x,y)data'
!*********************************** number***
! seg = 1 => y1
!*********************************************
! seg = 2 => ,for loop (1)
! ( 1 -> p,p+1-> n) ,y1,y2
! loop 1: for p=2 to n-2=n-(seg-1)*2 => p
!*********************************************
! seg = 3 => ,for loop (2)
! ( 1->p,p+1->q,q+1->n) ,y1,y2,y3
! loop 1: for p=2 to n-4=n-(seg-1)*2 => p
! loop 2: for q=p+2 to n-2=n-(seg-2)*2 => q
!********************************************
! seg = 4 => ,for loop (3)
! ( 1->p,p+1->q,q+1->r,r+1->n) ,y1,y2,y3,y4
! loop 1: for p=2 to n-6=n-(seg-1)*2 => p
! loop 2: for q=p+2 to n-4=n-(seg-2)*2 => q
! loop 3: for r=q+2 to n-2=n-(seg-3)*2 => r
! etc....
!********************************************
!********************************************
num_col = 3*seg
!********************************************
mat redim data_b(1,num_col)
for pp = 2 to n-(seg-1)*2 !n-4
for qq = pp+2 to n-(seg-2)*2 !n-2
!**********************************
! first : lin equation y=aa1*x+bb1
!**********************************
call x_mean(1,pp,xy(,),xm1)
call y_mean(1,pp,xy(,),ym1)
call s_xy(1,pp,xm1,ym1,xy(,),sxy1)
call s_xx(1,pp,xm1,ym1,xy(,),sxx1)
aa1 = sxy1/sxx1
bb1 = ym1-aa1*xm1
!**********************************
! second : lin equation y=aa2*x+bb2
!**********************************
call x_mean(pp+1,qq,xy(,),xm2)
call y_mean(pp+1,qq,xy(,),ym2)
call s_xy(pp+1,qq,xm2,ym2,xy(,),sxy2)
call s_xx(pp+1,qq,xm2,ym2,xy(,),sxx2)
aa2 = sxy2/sxx2
bb2 = ym2-aa2*xm2
!**********************************
!third : lin equation y=aa3*x+bb3
!**********************************
call x_mean(qq+1,n,xy(,),xm3)
call y_mean(qq+1,n,xy(,),ym3)
call s_xy(qq+1,n,xm3,ym3,xy(,),sxy3)
call s_xx(qq+1,n,xm3,ym3,xy(,),sxx3)
aa3 = sxy3/sxx3
bb3 = ym3-aa3*xm3
!******************************************
! find the minimum : qp,h1 = 1,for sort qp
!******************************************
num=num+1
mat redim data_b(num,num_col)
!********************
! least square method
!********************
call q_p_q_var(n,pp,qq,aa1,aa2,aa3,bb1,bb2,bb3,xy(,),qp)
data_b(num,1)=qp
data_b(num,2)=aa1
data_b(num,3)=bb1
data_b(num,4)=aa2
data_b(num,5)=bb2
data_b(num,6)=aa3
data_b(num,7)=bb3
data_b(num,8)=pp
data_b(num,9)=qq
next qq
next pp
h1 =1 ! first column for sort the rest of columns
nsort = num
for i=1 to nsort
!*************************************************
!after sort : use first row for values of 'data_b'
!*************************************************
call sort_qp(nsort,h1,data_b(,),num_col)
next i
! data from first row data_b
aa1=data_b(h1,2)
bb1=data_b(h1,3)
aa2=data_b(h1,4)
bb2=data_b(h1,5)
aa3=data_b(h1,6)
bb3=data_b(h1,7)
!************** end algo**************************
print
print "regression line,green => y1 = ";aa1;" * x(+)";bb1
print "regression line,red => y2 = ";aa2;" * x(+)";bb2
print "regression line,blue => y3 = ";aa3;" * x(+)";bb3
print
!******************************************
! find from y1=y2 x=xx_1
! aa1*x+bb1 = aa2*x+bb2 =>(x1,y1)
!******************************************
xx_1 = (bb1-bb2)/(aa2-aa1)
person1 = xx_1
tsolve1 = aa1*person1+bb1
print "**************************"
print "*inter section (y1,y2) =>*"
print "**************************"
print
print " x1-intersect = ";person1
print " y1-intersect = ";tsolve1
print
!******************************************
! find from y2=y3 x=xx_2
! aa2*x+bb2 = aa3*x+bb3 =>(x2,y2)
!******************************************
print
xx_2 = (bb2-bb3)/(aa3-aa2)
person2 = xx_2
tsolve2 = aa2*person2+bb2
print "**************************"
print "*inter section (y2,y3) =>*"
print "**************************"
print
print " x2-intersect = ";person2
print " y2-intersect = ";tsolve2
!*******************************************
window #2 ! draw graph right part of screen
!*******************************************
ntx = nt_x ! tick marks x-axis
nty = nt_y ! tick marks y-axis
xp(1)=tb_x ! xmin,cm
xp(2)=te_x ! xmax,cm
yp(1)=tb_y ! ymin,kg
yp(2)=te_y ! ymax,kg
ask pixels px,py
asp_r = px/py !1 !aspect ratio correction(rectangular window)
set window xp(1)*asp_r,xp(2)*asp_r,yp(1),yp(2)
call draw_axes(xp(1)*asp_r,xp(2)*asp_r,yp(1),yp(2),ntx,nty)
set color "black"
!**********************************************************
call num_txt_axes(title$,x_txt$,y_txt$,xp(),yp(),txt,ntx,nty)
!**********************************************************
! draw graph in the window x-axes ,y-axes(right part)
!**********************************************************
z3$ = "red" ! color 'circle'
set color z3$ ! color of the graph
!*************************************************
! data x,y,draw circles
!*************************************************
s1=0.15 ! radius
! you can change this to get circles
asp = py/px ! to draw circle in the graph for data xy
! draw n circles for data xy
for i=1 to n
call smallcir_dr(s1,xy(1,i),xy(2,i),scale,asp,z3$)
flood xy(1,i),xy(2,i) ! fill circle with color 'red'
!*****************************
! put a number for every point
!*****************************
plot text,at xy(1,i)+s1+s1/2,xy(2,i)+s1+s1/2:str$(i)
next i
!*************************************************
! axes ,transformation x=X-avgx,y=Y-avgy
!*************************************************
set color "green"
plot lines :avgx,yp(1);avgx,yp(2)
plot lines :xp(1),avgy;xp(2),avgy
plot text,at avgx,avgy:"(mean(x),mean(y))"
!*************************************************
! draw regression line : y = aa*x+bb,2 points
!*************************************************
set color "blue"
! first point xx1,yy1
xx1 = tb_x
yy1 = aa*xx1+bb
! second point xx2,yy2
xx2 = te_x
yy2 = aa*xx2+bb
plot lines :xx1,yy1;xx2,yy2
!*************************************************
! to find : intersection y1,y2,y3,trimodal fit
!*************************************************
! draw regression line : y1 = aa1*x+bb1,2 points
!*************************************************
set color "green"
! first point xx1,yy1
xx1 = tb_x
yy1 = aa1*xx1+bb1
! second point xx2,yy2
xx2 = te_x
yy2 = aa1*xx2+bb1
plot lines :xx1,yy1;xx2,yy2
!*************************************************
! draw regression line : y2 = aa2*x+bb2,2 points
!*************************************************
set color "red"
! first point xx1,yy1
xx1 = tb_x
yy1 = aa2*xx1+bb2
! second point xx2,yy2
xx2 = te_x
yy2 = aa2*xx2+bb2
plot lines :xx1,yy1;xx2,yy2
!*************************************************
! draw regression line : y3 = aa3*x+bb3,2 points
!*************************************************
set color "blue"
! first point xx1,yy1
xx1 = tb_x
yy1 = aa3*xx1+bb3
! second point xx2,yy2
xx2 = te_x
yy2 = aa3*xx2+bb3
plot lines :xx1,yy1;xx2,yy2
!*************************************************
! 1e)intersection (y1,y2),2e)intersection (y2,y3)
!*************************************************
set color "magenta"
plot lines :person1,yp(1);person1,tsolve1
plot lines :xp(1),tsolve1;person1,tsolve1
cord1$ = "(" & str$(person1) & "," & str$(tsolve1)&")"
plot text,at person1,tsolve1:cord1$
set color "magenta"
plot lines :person2,yp(1);person2,tsolve2
plot lines :xp(1),tsolve2;person2,tsolve2
cord2$ = "(" & str$(person2) & "," & str$(tsolve2)&")"
plot text,at person2,tsolve2:cord2$
!*********************************************************
! data : x,y
!*********************************************************
!*********************************************************
! data statement, x values
!*********************************************************
data 1.2,2.8,5.2,7.4,10,11.6,12.6,14,15.4,16.2,17,18.6,19.1
data 19.3,19.5,20.3
!*********************************************************
! data statement, y values * 1000
!*********************************************************
data 6.523,5.683,4.722,3.841,2.759,2.201,1.943,1.864,2.155
data 2.336,2.493,2.822,3.351,4.821,5.256,6.231,7.235,7.961
!*********************************************************
end
!******************************************
! y = aa*x+bb 'regression line
!******************************************
sub regr_line_aa(n,spxy,spxx,spyy,sx,sy,aa)
tel1 = n*spxy-sx*sy
noem1= n*spxx - (sx)^2
aa = tel1/noem1
end sub
sub regr_line_bb(n,avgx,avgy,aa,bb)
bb = avgy - aa*avgx
end sub
!*****************************************
!*****************************************
sub sum_val_average(ii,n,xy(,),somval,avg)
som1 = 0
for i= 1 to n
som1 = som1 + xy(ii,i)
next i
somval=som1
avg = (somval)/n
end sub
sub sum_prod_x_y(n,xy(,),somprodxy)
som2 = 0
for i=1 to n
som2=som2+xy(1,i)*xy(2,i)
next i
somprodxy = som2
end sub
sub sum_prod_x_x(ii,n,xy(,),somprodxx)
som3 = 0
for i=1 to n
som3=som3+xy(ii,i)*xy(ii,i)
next i
somprodxx = som3
end sub
sub correlation_coef(spxy,spxx,spyy,r_cor)
tel2 = spxy
noem2 = sqr(spxx*spyy)
r_cor = tel2/noem2
end sub
!**********************************************
!**** trimodal fit*****************************
!**********************************************
!**********************************************
!1e) mean x,y => !
!**********************************************
sub x_mean(t1,t2,xy(,),xm)
som4 = 0
for i = t1 to t2
som4 = som4 + xy(1,i)
next i
xm = som4/(t2-t1+1)
end sub
sub y_mean(t1,t2,xy(,),ym)
som5 = 0
for i = t1 to t2
som5 = som5 + xy(2,i)
next i
ym = som5/(t2-t1+1)
end sub
!*********************************************
! components needed : regression line
!*********************************************
sub s_xy(v1,v2,xm,ym,xy(,),sxy)
som6=0
for i=v1 to v2
som6 = som6 +(xy(1,i)-xm)*(xy(2,i)-ym)
next i
sxy = som6
end sub
sub s_xx(v1,v2,xm,ym,xy(,),sxx)
som8 = 0
for i=v1 to v2
som8 = som8 +(xy(1,i)-xm)^2
next i
sxx = som8
end sub
sub s_yy(v1,v2,xm,ym,xy(,),syy)
som9 = 0
for i=v1 to v2
som9 = som9 +(xy(2,i)-ym)^2
next i
syy = som9
end sub
!**********************************************
! coef :regression line => y = a_x*x + b_c
!**********************************************
sub regres_a(sxx,syy,sxy,a_x)
!*****************
! y = a_x*x + b_c :regression line
!*****************
a_x = sxy/sxx
end sub
sub regres_b(a_x,xm,ym,b_c)
!****************
! y = a_x*x + b_c :regression line
!****************
b_c = xm - a_x*ym
end sub
!***********************************************
!***********************************************
!*********************************************
! find minimum from all qp , pp=2..n-2
! 2 : reason (you can't draw a line through
! one point.)
!*********************************************
sub q_p_q_var(n,pp,qq,aa_1,aa_2,aa_3,bb_1,bb_2,bb_3,xy(,),qp)
som10 = 0
som11 = 0
som12 = 0
for j=1 to pp
som10 = som10 + (xy(2,j)-aa_1*xy(1,j)-bb_1)^2
next j
for j=pp+1 to qq
som11 = som11 + (xy(2,j)-aa_2*xy(1,j)-bb_2)^2
next j
for j=qq+1 to n
som12 = som12 + (xy(2,j)-aa_3*xy(1,j)-bb_3)^2
next j
qp = som10 + som11 + som12
end sub
!*********************************************
! sort values : qp find minimum
!*********************************************
sub sort_qp(ns,h1,aq(,),num_col)
!****************************************
!****************************************
! modification : 2D array
!****************************************
! 'small to big' sorting
!****************************************
! see also : website,computer programma's
!****************************************
!****************************************
flag=0
do
flag=1 ! otherwise exit do
for i=1 to ns-1
if aq(i,h1)> aq(i+1,h1) then !
!***********************************
! sort on column h1
!***********************************
call swap(aq(,),i,h1)
!**********************************
! columns update for the rest
!**********************************
for j=1 to num_col
if j <> h1 then
call swap(aq(,),i,j)
end if
next j
flag=0 ! no more swap
end if
next i
loop while (flag=0)
!****************************************
end sub
sub swap(aq(,),i,h1)
!===swap=======
temp =aq(i,h1)
aq(i,h1) =aq(i+1,h1)
aq(i+1,h1)=temp
!==============
end sub
!*********************************************
!*********************************************
sub smallcir_dr(s1,xf,yf,scale,asp,z3$)
set color z3$
for theta=0 to 2*3.14 step 0.0001
plot points:s1*cos(theta)+xf,s1*sin(theta)*asp+yf
next theta
end sub
sub num_txt_axes(title$,x_txt$,y_txt$,xp(),yp(),txt,ntx,nty)
mx = xp(2)-xp(1)
my = yp(2)-yp(1)
plot text,at xp(1)+mx/txt,yp(2)-my/(txt*txt):title$
!*****************************************
!**** range x axis 0 to xmax *******
!*****************************************
for i=2 to ntx
m1 = xp(1)+mx/ntx*(i-1)
plot text,at m1,yp(1):str$(m1)
next i
plot text,at xp(2)-txt,yp(1)+txt/4:x_txt$
!*****************************************
!**** range y axis 0 to ymax ************
!*****************************************
for i=2 to nty
m2 = yp(1)+my/nty*(i-1)
plot text,at xp(1),m2:str$(m2)
next i
plot text,at xp(1)+txt,yp(2)-txt/4:y_txt$
end sub
SUB draw_axes(xmin,xmax,ymin,ymax,ntx,nty)
! axis :for distribution function
! see extra,website
! distance between tick marks on x-axis
LET dx = (xmax - xmin)/ntx
! distance between tick marks on y-axis
LET dy = (ymax - ymin)/nty
SET WINDOW xmin ,xmax ,ymin ,ymax
LET x0 = max(0,xmin)
LET y0 = max(0,ymin)
IF ymin*ymax < 0 then
LET y0 = 0
ELSE
LET y0 = ymin
END IF
PLOT LINES: xmin,y0;xmax,y0 ! horizontal axis
PLOT LINES: x0,ymin;x0,ymax ! vertical axis
LET tx = 0.1*dy ! size of tick mark
LET ty = 0.1*dx
FOR itick = 0 to ntx
LET x = xmin + itick*dx
PLOT LINES: x,y0 - tx; x,y0 + tx ! draw ticks on x axis
next itick
FOR itick = 0 to nty
LET y = ymin + itick*dy
PLOT LINES: x0 - ty,y; x0 + ty,y ! draw ticks on y axis
NEXT itick
END SUB