m_qtlmap_analyse_discret_unitrait
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_analyse_discret_unitrait
SYNOPSIS
DESCRIPTION
Analysis module for LA analysis with discrete data and model description
NOTES
SEE ALSO
funct_0qtl_discret_family
[ Top ] [ m_qtlmap_analyse_discret_unitrait ] [ Subroutines ]
NAME
funct_0qtl_discret_family
DESCRIPTION
likelihood function 0 QTL, 1 trait for one family
INPUTS
ip : index sire
jm : idex dam
n : number of parameter of the vector x
x : vector inputs
iuser : user parameters of integer type
user : user parameters of real type
OUTPUTS
f : result of the function
NOTES
MATH_QTLMAP_G01EAF
SOURCE
355 subroutine funct_0qtl_discret_family(ip,jm,n,x,f,iuser,user) 356 integer , intent(in) :: ip,jm,n 357 ! 358 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 359 ! 360 real (kind=dp) ,dimension(n), intent(in) :: x 361 real (kind=dp) ,intent(inout) :: f 362 integer ,dimension(1), intent(in) :: iuser 363 real (kind=dp) ,dimension(1), intent(in) :: user 364 365 ! Tableaux dimensionnes selon nm, le nombre de meres 366 real (kind=dp) ,dimension(nmod(current_ic)) :: threshold 367 368 !modif mkw 369 ! declaration de l comme entier et on a plus besoin de la variable vpf 370 371 integer :: i, kkd, ifail 372 integer ::jnit,neffet,ief,ico,ic 373 real (kind=dp) :: v, vpf,tig 374 !mkw 375 ! 376 neffet=np+nbnivest 377 378 !****************************************************************************** 379 ! 380 ifail=1 381 ic = current_ic 382 383 threshold(1)=x(neffet+1) 384 do i=2,nmod(ic)-1 385 threshold(i)=threshold(i-1)+x(neffet+i) 386 end do 387 ! 388 389 jnit=2 390 if (nbfem.eq.0)jnit=1 391 392 f=0.d0 393 tig=grand 394 if (x(ip)> 0)tig=1.d0/x(ip) 395 ! 396 ! on ne consid�e que les m�res 397 398 do kkd=ndm(jm)+1,ndm(jm+1) 399 if(presentc(ic,kkd)) then 400 v=0.d0 401 if(vecsol(nivdir(kkd,1))) v=v+x(np+corniv(nivdir(kkd,1))) 402 403 if(estime(ic,jm)) then 404 if(vecsol(nivdir(kkd,2))) v=v+x(np+corniv(nivdir(kkd,2))) 405 end if 406 407 do ief=jnit+1,jnit+nbef 408 if(vecsol(nivdir(kkd,ief))) v=v+x(np+corniv(nivdir(kkd,ief))) 409 end do 410 411 do ico=1,nbco 412 if(vecsol(ntnifix+ico)) v=v+covdir(kkd,ico)*x(np+corniv(ntnifix+ico)) 413 end do 414 415 if (ydiscretord(ic,kkd)==1) then 416 vpf=MATH_QTLMAP_G01EAF('L',((threshold(1)-v)*tig),ifail) 417 else 418 if (ydiscretord(ic,kkd)==nmod(ic)) then 419 vpf=1-MATH_QTLMAP_G01EAF('L',(threshold(nmod(ic)-1)-v)*tig,ifail) 420 else 421 vpf=MATH_QTLMAP_G01EAF('L',(threshold(ydiscretord(ic,kkd))-v)*tig,ifail)-& 422 MATH_QTLMAP_G01EAF('L',(threshold(ydiscretord(ic,kkd)-1)-v)*tig,ifail) 423 end if 424 end if 425 426 if (vpf <= 0) then 427 f=INIFINY_REAL_VALUE 428 return 429 else 430 f=f-dlog(vpf) 431 end if 432 433 end if ! si presentc 434 end do ! sur kkd 435 436 end subroutine funct_0qtl_discret_family
funct_0qtl_discret_unitrait
[ Top ] [ m_qtlmap_analyse_discret_unitrait ] [ Subroutines ]
NAME
funct_0qtl_discret_unitrait
DESCRIPTION
likelihood function 0 QTL, 1 trait
INPUTS
n : number of parameter of the vector x
x : vector inputs
iuser : user parameters of integer type
user : user parameters of real type
OUTPUTS
f : result of the function
NOTES
MATH_QTLMAP_G01EAF
SOURCE
235 subroutine funct_0qtl_discret_unitrait(n,x,f,iuser,user) 236 use m_qtlmap_analyse_lin_gen, only : nbnivest, nbfem,nbef,ntnifix,nbco,covdir,nivdir,corniv,vecsol 237 use m_qtlmap_analyse_gen, only : estime 238 239 implicit none 240 integer , intent(in) :: n 241 ! 242 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 243 ! 244 real (kind=dp) ,dimension(n), intent(in) :: x 245 real (kind=dp) ,intent(inout) :: f 246 integer ,dimension(1), intent(in) :: iuser 247 real (kind=dp) ,dimension(1), intent(in) :: user 248 249 ! Tableaux dimensionnes selon nm, le nombre de meres 250 real (kind=dp) ,dimension(nmod(current_ic)) :: threshold 251 252 !modif mkw 253 ! declaration de l comme entier et on a plus besoin de la variable vpf 254 255 integer :: jm, i, kkd, ifail 256 integer ::jnit,ip,neffet,ief,ico,ic 257 real (kind=dp) :: v, vpf,tig 258 !mkw 259 ! 260 neffet=np+nbnivest 261 262 !****************************************************************************** 263 ! 264 ifail=1 265 ic = current_ic 266 267 threshold(1)=x(neffet+1) 268 do i=2,nmod(ic)-1 269 threshold(i)=threshold(i-1)+x(neffet+i) 270 end do 271 ! 272 273 jnit=2 274 if (nbfem.eq.0)jnit=1 275 276 f=0.d0 277 do ip=1,np 278 tig=grand 279 if (x(ip)> 0)tig=1.d0/x(ip) 280 fp0(ip)=0.d0 281 282 do jm=nmp(ip)+1,nmp(ip+1) 283 fm0(jm)=0.d0 284 ! 285 ! on ne considere que les meres 286 287 do kkd=ndm(jm)+1,ndm(jm+1) 288 ! 289 if(presentc(ic,kkd)) then 290 291 v=0.d0 292 if(vecsol(nivdir(kkd,1))) v=v+x(np+corniv(nivdir(kkd,1))) 293 294 if(estime(ic,jm)) then 295 if(vecsol(nivdir(kkd,2))) v=v+x(np+corniv(nivdir(kkd,2))) 296 end if 297 298 do ief=jnit+1,jnit+nbef 299 if(vecsol(nivdir(kkd,ief))) v=v+x(np+corniv(nivdir(kkd,ief))) 300 end do 301 302 do ico=1,nbco 303 if(vecsol(ntnifix+ico)) v=v+covdir(kkd,ico)*x(np+corniv(ntnifix+ico)) 304 end do 305 306 if (ydiscretord(ic,kkd)==1) then 307 vpf=MATH_QTLMAP_G01EAF('L',((threshold(1)-v)*tig),ifail) 308 else 309 if (ydiscretord(ic,kkd)==nmod(ic)) then 310 vpf=1-MATH_QTLMAP_G01EAF('L',(threshold(nmod(ic)-1)-v)*tig,ifail) 311 else 312 vpf=MATH_QTLMAP_G01EAF('L',(threshold(ydiscretord(ic,kkd))-v)*tig,ifail)-& 313 MATH_QTLMAP_G01EAF('L',(threshold(ydiscretord(ic,kkd)-1)-v)*tig,ifail) 314 end if 315 end if 316 317 if (vpf <= 0) then 318 fm0(jm)=INIFINY_REAL_VALUE 319 else 320 fm0(jm)=fm0(jm)-dlog(vpf) 321 end if 322 323 end if ! si presentc 324 end do ! sur kkd 325 326 fp0(ip)=fp0(ip)+fm0(jm) 327 f=f+fm0(jm) 328 329 end do ! sur jm 330 331 end do ! sur ip 332 333 return 334 335 end subroutine funct_0qtl_discret_unitrait
funct_1qtl_discret_family
[ Top ] [ m_qtlmap_analyse_discret_unitrait ] [ Subroutines ]
NAME
funct_1qtl_discret_family
DESCRIPTION
likelihood function 1 QTL, 1 trait for a family
INPUTS
ip : index sire
jm : idex dam
n : number of parameter of the vector x
x : vector inputs
iuser : user parameters of integer type
user : user parameters of real type
OUTPUTS
f : result of the function
NOTES
structure used : vecsol,nivdir,corniv,covdir,presentc,ydiscretord,nmod,ngenom,ndesc,ngend function : estime,MATH_QTLMAP_G01EAF
SOURCE
841 subroutine funct_1qtl_discret_family(ip,jm,n,x,f,iuser,user) 842 843 integer , intent(in) :: ip,jm,n 844 ! 845 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 846 ! 847 real (kind=dp) ,dimension(n), intent(in) :: x 848 real (kind=dp) , intent(inout) :: f 849 integer , dimension(1), intent(in) :: iuser 850 real (kind=dp) ,dimension(1), intent(in) :: user 851 real (kind=dp) ,dimension(nmod(current_ic)) :: threshold 852 853 ! Tableaux dimensionnes selon nm, le nombre de meres 854 real (kind=dp) :: effm 855 856 857 ! 858 ! Divers 859 integer :: init,jnit,nm1,nm2,ngeno1,ngeno2,ig,chr,ic 860 integer :: nd1,nd2,neffet,kd,kkd,ilev,ief,ico,lambda 861 integer :: ntnivmax,neffmax, m, i, j, m1, temp, k, ii, ifail 862 863 864 real (kind=dp) :: sig,var,vmere,vpf,v,wpf,tig 865 866 !****************************************************************************** 867 ! 868 chr =current_chr 869 ic = current_ic 870 neffet=np+nbnivest 871 ifail=1 872 threshold(1)=x(neffet+1) 873 do i=2,nmod(ic)-1 874 threshold(i)=threshold(i-1)+x(neffet+i) 875 end do 876 877 init=3 878 jnit=4 879 if (nbfem.eq.0)then 880 init=2 881 jnit=2 882 end if 883 884 f=0.d0 885 tig=grand 886 if (x(ip)> 0)tig=1.d0/x(ip) 887 vmere=0.d0 888 do ig=ngenom(chr,jm)+1,ngenom(chr,jm+1) 889 vpf=grand 890 do kd=ngend(chr,ig)+1,ngend(chr,ig+1) 891 kkd=ndesc(chr,kd) 892 893 if(presentc(ic,kkd)) then 894 v=0.d0 895 896 if(vecsol(nivdir(kkd,1))) v=v+x(np+corniv(nivdir(kkd,1)))*ppt(kd) 897 898 if(vecsol(nivdir(kkd,init))) v=v+x(np+corniv(nivdir(kkd,init))) 899 900 if(estime(ic,jm)) then 901 if(vecsol(nivdir(kkd,2))) v=v+x(np+corniv(nivdir(kkd,2)))*pmt(kd) 902 if(vecsol(nivdir(kkd,4))) v=v+x(np+corniv(nivdir(kkd,4))) 903 end if 904 905 do ief=jnit+1,jnit+nbef 906 if(vecsol(nivdir(kkd,ief))) v=v+x(np+corniv(nivdir(kkd,ief))) 907 end do 908 909 do ico=1,nbco 910 if(vecsol(ntnifix+ico)) v=v+covdir(kkd,ico)*x(np+corniv(ntnifix+ico)) 911 end do 912 913 if (ydiscretord(ic,kkd)==1) then 914 wpf=MATH_QTLMAP_G01EAF('L',((threshold(1)-v)*tig),ifail) 915 else 916 if (ydiscretord(ic,kkd)==nmod(ic)) then 917 wpf=1.d0-MATH_QTLMAP_G01EAF('L',(threshold(nmod(ic)-1)-v)*tig,ifail) 918 else 919 wpf=MATH_QTLMAP_G01EAF('L',(threshold(ydiscretord(ic,kkd))-v)*tig,ifail)-& 920 MATH_QTLMAP_G01EAF('L',(threshold(ydiscretord(ic,kkd)-1)-v)*tig,ifail) 921 end if 922 end if 923 924 vpf=vpf*wpf 925 926 end if ! presentc 927 end do ! kd 928 vmere=vmere+probg(chr,ig)*vpf 929 end do !ig 930 if (vmere <= 0) then 931 f=INIFINY_REAL_VALUE 932 else 933 f=dlog(grand)-dlog(vmere) 934 end if 935 936 end subroutine funct_1qtl_discret_family
funct_1qtl_discret_unitrait
[ Top ] [ m_qtlmap_analyse_discret_unitrait ] [ Subroutines ]
NAME
funct_1qtl_discret_unitrait
DESCRIPTION
likelihood function 1 QTL, 1 trait
INPUTS
n : number of parameter of the vector x
x : vector inputs
iuser : user parameters of integer type
user : user parameters of real type
OUTPUTS
f : result of the function
NOTES
structure used : vecsol,nivdir,corniv,covdir,presentc,ydiscretord,nmod,ngenom,ndesc,ngend function : estime,MATH_QTLMAP_G01EAF
SOURCE
708 subroutine funct_1qtl_discret_unitrait(n,x,f,iuser,user) 709 use m_qtlmap_analyse_lin_gen, only : nbnivest 710 use m_qtlmap_analyse_gen, only : estime 711 712 implicit none 713 714 integer , intent(in) :: n 715 ! 716 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 717 ! 718 real (kind=dp) ,dimension(n), intent(in) :: x 719 real (kind=dp) , intent(inout) :: f 720 integer , dimension(1), intent(in) :: iuser 721 real (kind=dp) ,dimension(1), intent(in) :: user 722 real (kind=dp) ,dimension(nmod(current_ic)) :: threshold 723 724 ! Tableaux dimensionnes selon nm, le nombre de meres 725 real (kind=dp) ,dimension(nm) :: effm 726 727 728 ! 729 ! Divers 730 integer :: init,jnit,ip,nm1,nm2,jm,ngeno1,ngeno2,ig,chr,ic 731 integer :: nd1,nd2,neffet,kd,kkd,ilev,ief,ico,lambda 732 integer :: ntnivmax,neffmax, m, i, j, m1, temp, k, ii, ifail 733 734 735 real (kind=dp) :: sig,var,vmere,vpf,v,wpf,tig 736 737 !****************************************************************************** 738 ! 739 chr =current_chr 740 ic = current_ic 741 neffet=np+nbnivest 742 ifail=1 743 744 threshold(1)=x(neffet+1) 745 do i=2,nmod(ic)-1 746 threshold(i)=threshold(i-1)+x(neffet+i) 747 end do 748 749 init=3 750 jnit=4 751 if (nbfem.eq.0)then 752 init=2 753 jnit=2 754 end if 755 756 f=0.d0 757 do ip=1,np 758 tig=grand 759 if (x(ip)> 0)tig=1.d0/x(ip) 760 fp1(ip)=0.d0 761 do jm=nmp(ip)+1,nmp(ip+1) 762 vmere=0.d0 763 fm1(jm)=0.d0 764 do ig=ngenom(chr,jm)+1,ngenom(chr,jm+1) 765 vpf=grand 766 do kd=ngend(chr,ig)+1,ngend(chr,ig+1) 767 kkd=ndesc(chr,kd) 768 769 if(presentc(ic,kkd)) then 770 v=0.d0 771 772 if(vecsol(nivdir(kkd,1))) v=v+x(np+corniv(nivdir(kkd,1)))*ppt(kd) 773 774 if(vecsol(nivdir(kkd,init))) v=v+x(np+corniv(nivdir(kkd,init))) 775 776 if(estime(ic,jm)) then 777 if(vecsol(nivdir(kkd,2))) v=v+x(np+corniv(nivdir(kkd,2)))*pmt(kd) 778 if(vecsol(nivdir(kkd,4))) v=v+x(np+corniv(nivdir(kkd,4))) 779 end if 780 781 do ief=jnit+1,jnit+nbef 782 if(vecsol(nivdir(kkd,ief))) v=v+x(np+corniv(nivdir(kkd,ief))) 783 end do 784 785 do ico=1,nbco 786 if(vecsol(ntnifix+ico)) v=v+covdir(kkd,ico)*x(np+corniv(ntnifix+ico)) 787 end do 788 789 if (ydiscretord(ic,kkd)==1) then 790 wpf=MATH_QTLMAP_G01EAF('L',((threshold(1)-v)*tig),ifail) 791 else 792 if (ydiscretord(ic,kkd)==nmod(ic)) then 793 wpf=1.d0-MATH_QTLMAP_G01EAF('L',(threshold(nmod(ic)-1)-v)*tig,ifail) 794 else 795 wpf=MATH_QTLMAP_G01EAF('L',(threshold(ydiscretord(ic,kkd))-v)*tig,ifail)-& 796 MATH_QTLMAP_G01EAF('L',(threshold(ydiscretord(ic,kkd)-1)-v)*tig,ifail) 797 end if 798 end if 799 800 vpf=vpf*wpf 801 802 end if ! presentc 803 end do ! kd 804 vmere=vmere+probg(chr,ig)*vpf 805 end do !ig 806 if (vmere <= 0) then 807 fm1(jm)=INIFINY_REAL_VALUE 808 else 809 fm1(jm)=dlog(grand)-dlog(vmere) 810 end if 811 fp1(ip)=fp1(ip)+fm1(jm) 812 f=f+fm1(jm) 813 814 end do ! jm 815 end do ! ip 816 817 return 818 end subroutine funct_1qtl_discret_unitrait
opti_0qtl_discret_unitrait
[ Top ] [ m_qtlmap_analyse_discret_unitrait ] [ Subroutines ]
NAME
opti_0qtl_discret_unitrait
DESCRIPTION
Compute the likelihood 0 QTL, 1 trait , fixed effect and covaruiate included
INPUTS
ic : index trait
NOTES
contingence,precision,set_filter_optim,minimizing_funct_family,funct_0qtl_discret_family
SOURCE
102 subroutine opti_0qtl_discret_unitrait(ic) 103 use m_qtlmap_analyse_lin_gen, only : contingence,precision,vecsol,precis, & 104 nbnivest,ntniv,par0,precis0,vecsol0,xx 105 106 implicit none 107 108 integer , intent(in) :: ic 109 ! 110 ! Divers 111 112 integer :: iuser(1) 113 real (kind=dp) ,dimension(:),allocatable :: par,borni,borns 114 real (kind=dp) :: user(1) 115 real (kind=dp) :: f,eff,cumul 116 117 integer :: npar,ibound,ideb,ip,ix,ifail,liw,lw,i,indest,ntot 118 integer :: km,jm,liwx,lwx, j, m, m1, temp, k, kkd, ii 119 logical itest 120 logical found 121 logical , dimension(:,:,:),pointer :: filter_inc 122 123 !****************************************************************************** 124 ! 125 ! recherche des elements estimables a partir de la decomposition de choleski 126 ! comptage et recodification de ces elements (on commence par mettre a 0 les 127 ! compteurs destines aux tests des effes fixes 128 ! 129 130 itest=.false. 131 call contingence(ic,0,itest,.false.) 132 call precision(xx,precis) 133 134 current_ic=ic 135 npar=np+nbnivest+nmod(ic)-1 136 allocate (borni(npar)) 137 allocate (borns(npar)) 138 allocate (par(npar)) 139 !MODIF - OPMIZATION 140 allocate (filter_inc(np,nm, npar)) 141 call set_filter_optim(ic,.true.,.false.,ntnivmax,ntniv,vecsol,xinc,filter_inc) 142 !FIN MODIF - OPMIZATION 143 144 ibound=0 145 do ip=1,np 146 borni(ip)=SIG_MIN 147 borns(ip)=SIG_MAX 148 end do 149 do ix=np+1,npar-nmod(ic)+2 150 borni(ix)=XMU_MIN 151 borns(ix)=XMU_MAX 152 end do 153 do ix=npar-nmod(ic)+3,npar 154 borni(ix)=SIG_MIN 155 borns(ix)=XMU_MAX 156 end do 157 ! 158 ! Point de depart 159 do ip=1,np 160 par(ip)=sig0(ip) 161 end do 162 do ix=np+1,np+nbnivest 163 par(ix)=0.d0 164 end do 165 par(np+nbnivest+1)=seuil(ic,1) 166 do ix=2,nmod(ic)-1 167 par(np+nbnivest+ix)=seuil(ic,ix)-seuil(ic,ix-1) 168 end do 169 ! 170 ! 171 ! Optimisation de la vraisemblance 172 ifail=1 173 ! call minimizing_funct(npar,ibound,funct_0qtl_discret_unitrait,borni,borns,par,f,iuser,user,ifail) 174 call minimizing_funct_family(npar,ibound,funct_0qtl_discret_family,filter_inc,fm0,fp0,borni,borns,par,f,iuser,user,ifail) 175 f0=f 176 177 do i=1,npar 178 par0(i)=par(i) 179 end do 180 181 do i=1,ntniv 182 vecsol0(i)=vecsol(i) 183 precis0(i)=precis(i) 184 end do 185 186 do ip = 1,np 187 sig1(ip)=par(ip) 188 end do 189 190 xmu1g=par(np+1) 191 indest=1 192 do ip = 1,np 193 if (vecsol(1+ip)) then 194 indest=indest+1 195 xmu1p(ip)=par(np+indest) 196 end if 197 end do 198 199 ntot=np+1 200 km=0 201 do jm=1,nm 202 if (estime(ic,jm).and.opt_sib.eq.2)then 203 km=km+1 204 if(vecsol(ntot+km)) then 205 indest=indest+1 206 xmu1m(jm)=par(np+indest) 207 end if 208 end if 209 end do 210 211 deallocate (borni) 212 deallocate (borns) 213 deallocate (par) 214 deallocate (filter_inc) 215 216 end subroutine opti_0qtl_discret_unitrait
opti_1qtl_discret_unitrait
[ Top ] [ m_qtlmap_analyse_discret_unitrait ] [ Subroutines ]
NAME
opti_1qtl_discret_unitrait
DESCRIPTION
Compute the test statistic among the chromosome : 1 QTL, 1 trait
INPUTS
ic : index trait
OUTPUTS
lrtsol : likelihood ratio test information (see TYPE_LRT_SOLUTION) fmax : value of the likelihood function at the maximum supnbnivest : number of paramter estimated in the solution (see contingence)
NOTES
prepinc,contingence,precision,set_filter_optim,minimizing_funct_family,funct_1qtl_discret_family
SOURCE
455 subroutine opti_1qtl_discret_unitrait(ic,lrtsol,fmax,supnbnivest) 456 457 implicit none 458 459 integer , intent(in) :: ic 460 type(TYPE_LRT_SOLUTION) , intent(out) :: lrtsol 461 real (kind=dp) , intent(out) :: fmax 462 integer , intent(out) :: supnbnivest 463 ! 464 465 ! Divers 466 integer :: iuser(1) 467 real (kind=dp) ,dimension(:),allocatable :: par,borni,borns 468 real (kind=dp),dimension(:),allocatable :: val 469 real (kind=dp) :: user(1) 470 real (kind=dp) :: f,eff,cumul 471 real (kind=dp) :: xlrt_t,f1 472 logical :: found 473 474 logical itest,stvecsol(size(vecsol1)) 475 integer :: ip,i,ideb,n,ilong,ibound,npar,ix,j,ifail,liw 476 integer :: ii,m, m1, k,temp,chr 477 logical , dimension(:,:,:),pointer :: filter_inc 478 479 480 call new(1,lrtsol) 481 482 current_ic = ic 483 484 ! 485 !****************************************************************************** 486 ! Calcul de la vraisemblance sous H1 487 488 ! 489 ! initialisation 490 ! on utilisera : 491 ! PAR (vecteur des param�tres � optimiser), 492 ! CORNIV (vecteur des positions, parmi NTNIV, des effets extimables), 493 ! SOLVEC (vecteur disant l'estimabilit� des effets) 494 ! VAL le vecteur complet (ntniv positions) des valeurs des niveaux des effets 495 ! STVAL,STVECSOL et STCORNIV les copies de VAL, VECSOL et CORNIV � l'it�ration 496 ! n-1 quand on analyse la position n 497 ! 498 ! 499 !****************************************************************************** 500 !****************************************************************************** 501 ! 502 ! transformation des donn�es discret � des donn�e utilisable par QTLMAP 503 ! 504 ! 505 !***************************************************************************** 506 !****************************************************************************** 507 ! 508 ! 509 ! 510 ! integer ,dimension(:),allocatable :: ydiscretord, indicemod 511 ! integer ::i, j, m, m1, temp, k, nmod, ii 512 ! Point de depart 513 514 allocate ( val( np+ntnivmax +nmod(ic)-1 ) ) 515 allocate ( par( np+ntnivmax +nmod(ic)-1) ) 516 allocate ( borni( np+ntnivmax +nmod(ic)-1) ) 517 allocate ( borns( np+ntnivmax +nmod(ic)-1) ) 518 519 !MODIF - OPMIZATION OFI 520 allocate (filter_inc(np,nm, np+ntnivmax +nmod(ic)-1)) 521 !FIN MODIF - OPMIZATION OFI 522 523 do ip=1,np 524 par(ip)=sig1(ip) 525 end do 526 527 par(np+1)=xmu1g 528 529 do i=np+2,size(par) 530 par(i)=0.d0 531 end do 532 533 do i=np+1,size(val) 534 val(i-np)=par(i) 535 end do 536 ! 537 do i=1,size(vecsol) 538 vecsol(i)=.true. 539 end do 540 541 542 lrtsol%lrtmax(0)=-1.d75 543 544 do chr=1,nchr 545 n=0 546 current_chr = chr 547 ilong=get_ilong(chr) 548 do ix=0,ilong,pas 549 n=n+1 550 ! 551 ! on stocke les conditions en n 552 do i=1,ntniv 553 stvecsol(i)=vecsol(i) 554 end do 555 556 call prepinc(1,n,ic,"LA ") 557 itest=.false. 558 call contingence(ic,1,itest,.false.) 559 npar=np+nbnivest+nmod(ic)-1 560 561 !MODIF - OPMIZATION 562 call set_filter_optim(ic,.true.,.false.,ntnivmax,ntniv,vecsol,xinc,filter_inc) 563 !FIN MODIF - OPMIZATION 564 565 566 ibound=0 567 do ip=1,np 568 borni(ip)=SIG_MIN 569 borns(ip)=SIG_MAX 570 end do 571 do i=np+1,npar-nmod(ic)+2 572 borni(i)=XMU_MIN 573 borns(i)=XMU_MAX 574 end do 575 do i=npar-nmod(ic)+3,npar 576 borni(i)=SIG_MIN 577 borns(i)=XMU_MAX 578 end do 579 ! 580 ! Point de depart 581 do ip=1,np 582 par(ip)=sig1(ip) 583 end do 584 do i=np+1,np+nbnivest 585 par(i)=0.d0 586 end do 587 588 j=np 589 do i=1,ntniv 590 if(vecsol(i))then 591 j=j+1 592 par(j)=0.d0 593 if(stvecsol(i)) par(j)=val(i) 594 end if 595 stvecsol(i)=vecsol(i) 596 end do 597 598 par(np+nbnivest+1)=seuil(ic,1) 599 do i=2,nmod(ic)-1 600 par(np+nbnivest+i)=seuil(ic,i)-seuil(ic,i-1) 601 end do 602 603 604 ! on stocke les conditions en n 605 !call prepinc(n,ic) 606 !itest=.false. 607 !call contingence(ic,1,itest) 608 609 ! Optimisation de la vraisemblance a la position dx 610 611 ifail=1 612 613 !call minimizing_funct(npar,ibound,funct_1qtl_discret_unitrait,borni,borns,par,f1,iuser,user,ifail) 614 call minimizing_funct_family(npar,ibound,funct_1qtl_discret_family,filter_inc,fm1,fp1,borni,borns,par,f1,iuser,user,ifail) 615 616 j=np 617 do i=1,ntniv 618 if(vecsol(i)) then 619 j=j+1 620 val(i)=par(j) 621 else 622 val(i)=9999.d0 623 end if 624 end do 625 ! 626 ! preparation de la matrice d'incidence 627 ! 628 ! 629 ! recherche des elements estimables � partir de la d�composition de choleski 630 ! comptage et recodification de ces elements (on commence par mettre � 0 les 631 ! compteurs destin�s aux tests des effes fix�s 632 633 if ( f1 < INIFINY_REAL_VALUE ) then 634 xlrt_t=-2.d0*(f1-f0) 635 else 636 xlrt_t=0 637 end if 638 639 640 ! on met les profil / prog�niteur dnas les ficheir ad hoc 641 642 do ii=1,np 643 lrtsol%xlrp(chr,ii,n)=-2.d0*(fp1(ii)-fp0(ii)) 644 end do 645 do ii=1,nm 646 lrtsol%xlrm(chr,ii,n)=-2.d0*(fm1(ii)-fm0(ii)) 647 end do 648 649 ! on stocke la position si elle est meilleure que les pr�c�dents 650 ! 651 if(lrtsol%lrtmax(0) < xlrt_t) then 652 653 lrtsol%chrmax(0)=chr 654 lrtsol%nxmax(0)=n 655 lrtsol%lrtmax(0)=xlrt_t 656 fmax=f1 657 supnbnivest=nbnivest 658 659 do i=1,npar 660 par1(i)=par(i) 661 end do 662 663 end if 664 ! 665 ! on garde les valeurs du LRT pour pouvopir dessiner la courbe de vraisemblance 666 ! 667 lrtsol%lrt1(chr,n)=xlrt_t 668 end do ! ix 669 end do ! chr 670 ! 671 ! on calcule la pr�cision des estimation au point correspondant 672 ! au LRT maximum 673 ! 674 call prepinc(lrtsol%chrmax(0),lrtsol%nxmax(0),ic,"LA ") 675 call contingence(ic,1,itest,.false.) 676 call precision(xx,precis) 677 678 do i=1,ntniv 679 precis1(i)=precis(i) 680 vecsol1(i)=vecsol(i) 681 end do 682 683 deallocate ( val ) 684 deallocate ( par ) 685 deallocate ( borni ) 686 deallocate ( borns ) 687 deallocate (filter_inc) 688 689 end subroutine opti_1qtl_discret_unitrait
set_solution_hypothesis0
[ Top ] [ m_qtlmap_analyse_discret_unitrait ] [ Subroutines ]
NAME
set_solution_hypothesis0
DESCRIPTION
fill the solution under H0 in the structure TYPE_INCIDENCE_SOLUTION
INPUTS
ic : index trait
OUTPUTS
incsol : the variable to fill
NOTES
set_solution_hypothesis1
[ Top ] [ m_qtlmap_analyse_discret_unitrait ] [ Subroutines ]
NAME
set_solution_hypothesis1
DESCRIPTION
fill the solution under H1 in the structure TYPE_INCIDENCE_SOLUTION
INPUTS
ic : index trait
OUTPUTS
incsol : the variable to fill
NOTES
test_lin
[ Top ] [ m_qtlmap_analyse_discret_unitrait ] [ Subroutines ]
NAME
test_lin
DESCRIPTION
Test des differents effets de nuisance du modele par une LRT compare a une chi2 and write in the output result file.
INPUTS
chr : chromosome index
ic : index trait
est_moy : if the general mean is estimated
supnbnivest : number of paramter of the solution
fmax : value of the likelihood at the maximum position
nposx : index of the position (see absi)
NOTES