## 26 May 2013 # V9 with(LinearAlgebra): interface(rtablesize=infinity); intM:=proc(K::integer, C::integer) description "Creat interaction terms (matrix form)."; local M:=Matrix(K,C,symbol=mix), tolt,cols,rows; # Set constraints M(..,1):=0; M(1..2,..):=0; tolt:=convert(M,`+`); cols:=MTM:-sum(M,1); rows:=MTM:-sum(M,2); M(2,..):=-cols; M(..,1):=-rows; M(2,1):=tolt; M; end proc; indprob:=proc(x, K::integer, C::integer, tvarphi::truefalse:=false, tvarp::truefalse:=false, phiMix::truefalse:=false ) description "Compute the Probability for given capture history \$x\$."; local # phi,p, # Vectors of phi and p int0, # Matrix of interaction terms ind, # Find indices of where the '1's are fcap, # First capture li, # Last capture ans; # Desired probability global phi,p; # Vectors of phi and p; # take care of phi (nu) # SINGLE CLASS if tvarphi then phi:=Array(1..K,[seq(1/(1+exp(-nu[j])),j=1..K)]); # phi:=Array(1..K,[seq(nu[j],j=1..K)]); else phi:=Array(1..K,[seq(1/(1+exp(-nu)),j=1..K)]); # phi:=Array(1..K,[seq(nu,j=1..K)]); end if; phi(K):=0; # take care of p if tvarp then if phiMix then int0:=intM(K,C); p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+tau[j]+eta[c]+int0[j,c])))); # p:=Array(1..K,1..C,(j,c)->mu+tau[j]+eta[c]+int0[j,c]); else if C=1 then p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+tau[j])))); # p:=Array(1..K,1..C,(j,c)->mu+tau[j]); else p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+tau[j]+eta[c])))); # p:=Array(1..K,1..C,(j,c)->mu+tau[j]+eta[c]); end if; end if; else p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+eta[c])))); # p:=Array(1..K,1..C,(j,c)->mu+eta[c]); end if; p(1,..):=0; p:=subs({tau[2]=0,eta[1]=0},p); ind:=ArrayTools[SearchArray](Array(x)); fcap,li:=ind[1],ind[-1]; ans:= add( w[c]* add( (1-phi[d])* mul(phi[j],j=fcap..(d-1)) * mul(p[j,c]^(x[j])*(1-p[j,c])^(1-x[j]),j=(fcap+1)..d), d=li..nops(x) ), c=1..C); ans; end proc; myP:=proc( K::integer, C::integer, tvarphi::truefalse:=false, tvarp::truefalse:=false, phiMix::truefalse:=false, dosum::truefalse:=false, n::integer:=30) description "Test for Parameter Redundancy (for mixtures of \$p\$)."; local CH, # All Capture Histories from \$K\$ sample kappa, # Exhaustive Summary allpars, # Complete Set of Model Parameters DD1, # Derivative Matrix (Jacobian) of Kappa ## Numerical ## sys, # System of eqns myrand, # Random number generator phitrue, # Independent Set of phi phinum, # Random points phiparM, # Vector of phi ptrue, # Independent Set of p pnum, # Random points pparM, # Vector of p ws, # Random points wsnum, # Vector of w num, # Complete Set of Random Points DD1num, # Numerical Derivative Matrix rS, # Symbolic Rank rN # Numerical Rank ; # Compute Capture Histories CH:=combinat[permute]([seq(s\$K,s=[0,1])],K); CH:=CH[2..-1]; # Get Exhaustive Summary kappa:=map(indprob,CH,K,C,tvarphi,tvarp,phiMix); kappa:=subs({w[C]=1-add(w[i],i=1..(C-1))},kappa); # kappa:=subs({seq(w[i]=1/(1+exp(-wnew[i])),i=1..(C-1))},kappa); kappa:=kappa[2..-1]; # Get Jacobian allpars:=indets(kappa,name); # pars:=allpars minus {seq(w[i],i=1..(C-1))}; DD1:=VectorCalculus[Jacobian](Vector(kappa),convert(allpars,list)); # rS:=Rank(DD1); # Initialise parameters myrand:=evalf[2](rand(10..90)/100); # Take care of phi phitrue:=convert(phi,set) minus {0}; sys:={op(phitrue)=~seq(myrand(),i=1..nops(phitrue))}; phinum:=solve(sys,indets(phitrue,name)); phiparM:=subs(phinum,phi); phiparM:=evalf(phiparM); phiparM:=phiparM[1..-2]; # Take care of p ptrue:=convert(p,set): if tvarp then if phiMix then ptrue:=ptrue minus {0}; else ptrue:=ptrue[2..(K+C-1)]; end if; else ptrue:=ptrue minus {0}; end if; sys:=ptrue=~seq(myrand(),i=1..nops(indets(ptrue,name))); pnum:=solve({sys},indets(ptrue,name)); pparM:=subs(pnum,p); pparM:=evalf(pparM); pparM:=pparM(2..,..); # Take care of weights wsnum:=[seq(w[i]=myrand(),i=1..(C-1))]; ws:=convert(rhs~(wsnum),Vector); wsnum:=convert(wsnum,set); # Get num into DD1 num:=phinum union pnum union wsnum; DD1num:=eval(DD1,num); rN:=Rank(DD1num); # Print Summary # if dosum then printf("Vector for phi \n"); printf("%.4f\n\n",phiparM); printf("Parameter Matrix for p \n"); printf("%.4f\n\n",pparM); printf("With vector of weights \n"); printf("%.4f\n\n",ws); #printf("Number of Parameters = %d, Numerical Rank = %d, defficiency = %d.\n", nops(allpars), rN, nops(allpars)-rN); printf("Number of Parameters = %d.\n",nops(allpars)); printf("Numerical Rank = %d.\n",rN); printf("Defficiency = %d.\n",nops(allpars)-rN); printf("---------------------------------------------------------------\n"); end if; # Print Summary # -End- # return rN; end proc;