*
* KUMAC pour montrer l'estimation max L de ligne droite
*
his/cre/1d 100 'y1 (x1=1.) ' 50  0.  5.
his/cre/1d 101 'y2 (x1=2.) ' 50  0.  5.
his/cre/1d 102 'y3 (x1=3.) ' 50  0.  5.
his/cre/1d 103 'y4 (x1=4.) ' 50  0.  5.
his/cre/1d 200 'Intersection a' 50  -1.0 1.0
his/cre/1d 201 'Pente b       ' 50   0.5 1.5
his/cre/2d 202 'a vs  b       ' 50  -1.0 1.0 50 0.5 1.5
his/cre/1d 300 'Intersection a' 50  -1.0 1.0
his/cre/1d 301 'Pente b       ' 50   0.5 1.5
*                             Generer echantillon de mesures ligne droite
sigma r1=rndm(array(1000)) ; sigma r2=rndm(array(1000))
sigma r3=rndm(array(1000)) ; sigma r4=rndm(array(1000))
sigma alpha=0.25 ; sigma beta=1.0 ; sigma s=0.2
sigma x1=1. ; sigma x2=2. ; sigma x3=3. ; sigma x4=4.
sigma y1=sin(2.*pi*r1)*sqrt(-2.*log(r2))*s+alpha+beta*x1
sigma y2=cos(2.*pi*r1)*sqrt(-2.*log(r2))*s+alpha+beta*x2
sigma y3=sin(2.*pi*r3)*sqrt(-2.*log(r4))*s+alpha+beta*x3
sigma y4=cos(2.*pi*r3)*sqrt(-2.*log(r4))*s+alpha+beta*x4
vec/hfill y1 100 ; vec/hfill y2 101 ; vec/hfill y3 102 ; vec/hfill y4 103
*                             Calculer intersection et pente
sigma xy=.25*(x1*y1+x2*y2+x3*y3+x4*y4)
sigma xx=.25*(x1*x1+x2*x2+x3*x3+x4*x4)
sigma xb=.25*(x1+x2+x3+x4)
sigma yb=.25*(y1+y2+y3+y4)
sigma b = (xb*yb-xy)/(xb**2-xx)
sigma a = yb - b*xb
vec/hfill  a 200
vec/hfill  b 201
do i=1,1000
   call hfill(202,a([i]),b([i]),1.)
enddo
*                             Moyennes, variances, correllation echantillon
sigma am=vsum(a)/1000
sigma bm=vsum(b)/1000
sigma sa=sqrt(vsum((a-am)**2)/1000)
sigma sb=sqrt(vsum((b-bm)**2)/1000)
sigma ro=(vsum((a-am)*(b-bm))/1000)/sa/sb
*                             Variances, correllation distribution
sigma sae=(s/2.)*sqrt(xx/(xx-xb**2))
sigma sbe=(s/2.)*sqrt(1./(xx-xb**2))
sigma roe=(s/2.)**2*xb/(xb**2-xx)/sae/sbe
*                             Pente vs intersection avec son
*                             ellipse de concentration (variance predite)
zone 1 1
his/plot 202
symb a b -1 1.
vec/cre ela(100)
vec/cre elb(100)
sigma ela=array(100,-2.*sae#+2.*sae)
sigma elb=beta+sbe*sqrt(1.-roe**2)*sqrt(4.-ela**2/sae**2)+sbe*roe*ela/sae
sigma ela=ela+alpha
pline 100 ela elb
sigma ela=ela-alpha
sigma elb=beta-sbe*sqrt(1.-roe**2)*sqrt(4.-ela**2/sae**2)+sbe*roe*ela/sae
sigma ela=ela+alpha
pline 100 ela elb
wait
*                             Distributions marginales avec predictions
zone 1 2
his/plot 200
sigma xx=array(50,-.98#.98)
sigma yy=(1000.*.04/sqrt(2.*pi)/sae)*exp(-(xx-alpha)**2/(2.*sae**2))
his/put/cont 300 yy
his/plot 300 'cs'
his/plot 201
sigma xx=array(50,.51#1.49)
sigma yy=(1000.*.02/sqrt(2.*pi)/sbe)*exp(-(xx-beta)**2/(2.*sbe**2))
his/put/cont 301 yy
his/plot 301 'cs'