From bd2d83202887555d9c28d98f2c45835d09df5a5e Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Sun, 29 Mar 2015 12:06:50 -0700 Subject: [PATCH 01/19] Added qdev parameters --- src/admb-code/iscam.tpl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index ad4a550f..6268d647 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -1556,6 +1556,12 @@ PARAMETER_SECTION //init_bounded_vector log_m_nodes(1,n_m_devs,-5.0,5.0,m_dev_phz); init_bounded_vector log_m_nodes(1,nMdev,-5.0,5.0,Mdev_phz); + // |---------------------------------------------------------------------------------| + // | DEVIATIONS IN CATCHABILITY COEFFICIENTS ASSUMING A RANDOM WALK | + // |---------------------------------------------------------------------------------| + // | + init_bounded_matrix log_q_devs(1,nItNobs,1,n_it_nobs,-5.0,5.0,-1); + // |---------------------------------------------------------------------------------| // | CORRELATION COEFFICIENTS FOR AGE COMPOSITION DATA USED IN LOGISTIC NORMAL | // |---------------------------------------------------------------------------------| @@ -1630,7 +1636,7 @@ PARAMETER_SECTION vector varphi(1,ngroup); vector sig(1,ngroup); vector tau(1,ngroup); - vector sigma_r(1,ngroup); + vector sigma_r(1,ngroup); // |---------------------------------------------------------------------------------| // | MATRIX OBJECTS From 2723e76f1933c7ec67ece868bb887aa462e6f673 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Mon, 30 Mar 2015 10:03:59 -0700 Subject: [PATCH 02/19] Implemented time-varying q with random walk parameters estimated, revamped survey q controls --- fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl | 26 +-- fba/HALIBUT/DATA/Halibut2014.dat | 210 +++++++++++----------- src/admb-code/iscam.tpl | 134 ++++++++------ 3 files changed, 202 insertions(+), 168 deletions(-) diff --git a/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl b/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl index 94eaed55..3bb8bcc6 100644 --- a/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl +++ b/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl @@ -83,7 +83,7 @@ 1888 1888 1888 -## +## ## ———————————————————————————————————————————————————————————————————————————————————— ## ## TIME VARYING NATURAL MORTALIIY RATES ## ## ———————————————————————————————————————————————————————————————————————————————————— ## @@ -104,18 +104,20 @@ ## ———————————————————————————————————————————————————————————————————————————————————— ## ## ## -## ------------------------------------------------------------------------- ## -## PRIORS FOR SURVEY Q ## -## Prior type: ## -## 0 - uninformative prior ## -## 1 - normal prior density for log(q) ## -## 2 - random walk in q ## -## ------------------------------------------------------------------------- ## -2 # -number of surveys (nits) -0 0 # -prior type (see legend above +## ———————————————————————————————————————————————————————————————————————————————————— ## +## ABUNDANCE OBSERVATION MODELS +## ———————————————————————————————————————————————————————————————————————————————————— ## +## QTYPE: +## 0 = FIXED SURVEY Q (specify log(mean) for prior log(mean)) +## 1 = CONSTANT Q (use MLE for q and optional informative prior) +## 2 = RANDOM WALK Q (use prior mean & sd for penalized random walk) +## ———————————————————————————————————————————————————————————————————————————————————— ## +2 # -number of surveys (n_it_nobs) +2 2 # -QTYPE (see legend above) 0 0 # -prior log(mean) -0 0 # -prior sd -## ------------------------------------------------------------------------- ## +0 0 # -prior sd (set to 0 for uniformative prior) +2 2 # -Estimation phase +## ———————————————————————————————————————————————————————————————————————————————————— ## ## ## ------------------------------------------------------------------------- ## diff --git a/fba/HALIBUT/DATA/Halibut2014.dat b/fba/HALIBUT/DATA/Halibut2014.dat index 42cd972d..c52a7dbb 100644 --- a/fba/HALIBUT/DATA/Halibut2014.dat +++ b/fba/HALIBUT/DATA/Halibut2014.dat @@ -384,111 +384,111 @@ ## commercial_data: ## survey ## iyr index(it) gear area group sex log(SE) log(PE) timing - 1907 280 1 1 1 0 0.100 0.010 0.5 - 1910 271 1 1 1 0 0.100 0.010 0.5 - 1911 237 1 1 1 0 0.100 0.010 0.5 - 1912 176 1 1 1 0 0.100 0.010 0.5 - 1913 129 1 1 1 0 0.100 0.010 0.5 - 1914 124 1 1 1 0 0.100 0.010 0.5 - 1915 118 1 1 1 0 0.100 0.010 0.5 - 1916 137 1 1 1 0 0.100 0.010 0.5 - 1917 98 1 1 1 0 0.100 0.010 0.5 - 1918 96 1 1 1 0 0.100 0.010 0.5 - 1919 93 1 1 1 0 0.100 0.010 0.5 - 1920 96 1 1 1 0 0.100 0.010 0.5 - 1921 88 1 1 1 0 0.100 0.010 0.5 - 1922 73 1 1 1 0 0.100 0.010 0.5 - 1923 78 1 1 1 0 0.100 0.010 0.5 - 1924 74 1 1 1 0 0.100 0.010 0.5 - 1925 68 1 1 1 0 0.100 0.010 0.5 - 1926 67 1 1 1 0 0.100 0.010 0.5 - 1927 65 1 1 1 0 0.100 0.010 0.5 - 1928 58 1 1 1 0 0.100 0.010 0.5 - 1929 51 1 1 1 0 0.100 0.010 0.5 - 1930 46 1 1 1 0 0.100 0.010 0.5 - 1931 50 1 1 1 0 0.100 0.010 0.5 - 1932 60 1 1 1 0 0.100 0.010 0.5 - 1933 63 1 1 1 0 0.100 0.010 0.5 - 1934 62 1 1 1 0 0.100 0.010 0.5 - 1935 76 1 1 1 0 0.100 0.010 0.5 - 1936 71 1 1 1 0 0.100 0.010 0.5 - 1937 80 1 1 1 0 0.100 0.010 0.5 - 1938 88 1 1 1 0 0.100 0.010 0.5 - 1939 80 1 1 1 0 0.100 0.010 0.5 - 1940 81 1 1 1 0 0.100 0.010 0.5 - 1941 85 1 1 1 0 0.100 0.010 0.5 - 1942 90 1 1 1 0 0.100 0.010 0.5 - 1943 95 1 1 1 0 0.100 0.010 0.5 - 1944 110 1 1 1 0 0.100 0.010 0.5 - 1945 102 1 1 1 0 0.100 0.010 0.5 - 1946 101 1 1 1 0 0.100 0.010 0.5 - 1947 99 1 1 1 0 0.100 0.010 0.5 - 1948 99 1 1 1 0 0.100 0.010 0.5 - 1949 95 1 1 1 0 0.100 0.010 0.5 - 1950 95 1 1 1 0 0.100 0.010 0.5 - 1951 96 1 1 1 0 0.100 0.010 0.5 - 1952 110 1 1 1 0 0.100 0.010 0.5 - 1953 131 1 1 1 0 0.100 0.010 0.5 - 1954 133 1 1 1 0 0.100 0.010 0.5 - 1955 119 1 1 1 0 0.100 0.010 0.5 - 1956 129 1 1 1 0 0.100 0.010 0.5 - 1957 110 1 1 1 0 0.100 0.010 0.5 - 1958 121 1 1 1 0 0.100 0.010 0.5 - 1959 129 1 1 1 0 0.100 0.010 0.5 - 1960 132 1 1 1 0 0.100 0.010 0.5 - 1961 127 1 1 1 0 0.100 0.010 0.5 - 1962 115 1 1 1 0 0.100 0.010 0.5 - 1963 105 1 1 1 0 0.100 0.010 0.5 - 1964 100 1 1 1 0 0.100 0.010 0.5 - 1965 99 1 1 1 0 0.100 0.010 0.5 - 1966 100 1 1 1 0 0.100 0.010 0.5 - 1967 101 1 1 1 0 0.100 0.010 0.5 - 1968 103 1 1 1 0 0.100 0.010 0.5 - 1969 95 1 1 1 0 0.100 0.010 0.5 - 1970 91 1 1 1 0 0.100 0.010 0.5 - 1971 89 1 1 1 0 0.100 0.010 0.5 - 1972 78 1 1 1 0 0.100 0.010 0.5 - 1973 63 1 1 1 0 0.100 0.010 0.5 - 1974 61 1 1 1 0 0.100 0.010 0.5 - 1975 61 1 1 1 0 0.100 0.010 0.5 - 1976 55 1 1 1 0 0.100 0.010 0.5 - 1977 63 1 1 1 0 0.100 0.010 0.5 - 1978 71 1 1 1 0 0.100 0.010 0.5 - 1979 75 1 1 1 0 0.100 0.010 0.5 - 1980 94 1 1 1 0 0.100 0.010 0.5 - 1981 111 1 1 1 0 0.100 0.010 0.5 - 1982 127 1 1 1 0 0.100 0.010 0.5 - 1984 291 1 1 1 0 0.100 0.010 0.5 - 1985 357 1 1 1 0 0.036 0.010 0.5 - 1986 320 1 1 1 0 0.041 0.010 0.5 - 1987 321 1 1 1 0 0.040 0.010 0.5 - 1988 368 1 1 1 0 0.035 0.010 0.5 - 1989 358 1 1 1 0 0.025 0.010 0.5 - 1990 319 1 1 1 0 0.028 0.010 0.5 - 1991 318 1 1 1 0 0.038 0.010 0.5 - 1992 319 1 1 1 0 0.034 0.010 0.5 - 1993 373 1 1 1 0 0.099 0.010 0.5 - 1994 306 1 1 1 0 0.072 0.010 0.5 - 1995 330 1 1 1 0 0.036 0.010 0.5 - 1996 392 1 1 1 0 0.038 0.010 0.5 - 1997 405 1 1 1 0 0.025 0.010 0.5 - 1998 408 1 1 1 0 0.025 0.010 0.5 - 1999 393 1 1 1 0 0.023 0.010 0.5 - 2000 403 1 1 1 0 0.022 0.010 0.5 - 2001 362 1 1 1 0 0.041 0.010 0.5 - 2002 360 1 1 1 0 0.019 0.010 0.5 - 2003 329 1 1 1 0 0.018 0.010 0.5 - 2004 319 1 1 1 0 0.019 0.010 0.5 - 2005 296 1 1 1 0 0.017 0.010 0.5 - 2006 270 1 1 1 0 0.019 0.010 0.5 - 2007 252 1 1 1 0 0.020 0.010 0.5 - 2008 232 1 1 1 0 0.017 0.010 0.5 - 2009 222 1 1 1 0 0.018 0.010 0.5 - 2010 204 1 1 1 0 0.020 0.010 0.5 - 2011 198 1 1 1 0 0.015 0.010 0.5 - 2012 195 1 1 1 0 0.021 0.010 0.5 - 2013 179 1 1 1 0 0.017 0.010 0.5 - 2014 191 1 1 1 0 0.052 0.010 0.5 + 1907 280 1 1 1 0 0.100 0.050 0.5 + 1910 271 1 1 1 0 0.100 0.050 0.5 + 1911 237 1 1 1 0 0.100 0.050 0.5 + 1912 176 1 1 1 0 0.100 0.050 0.5 + 1913 129 1 1 1 0 0.100 0.050 0.5 + 1914 124 1 1 1 0 0.100 0.050 0.5 + 1915 118 1 1 1 0 0.100 0.050 0.5 + 1916 137 1 1 1 0 0.100 0.050 0.5 + 1917 98 1 1 1 0 0.100 0.050 0.5 + 1918 96 1 1 1 0 0.100 0.050 0.5 + 1919 93 1 1 1 0 0.100 0.050 0.5 + 1920 96 1 1 1 0 0.100 0.050 0.5 + 1921 88 1 1 1 0 0.100 0.050 0.5 + 1922 73 1 1 1 0 0.100 0.050 0.5 + 1923 78 1 1 1 0 0.100 0.050 0.5 + 1924 74 1 1 1 0 0.100 0.050 0.5 + 1925 68 1 1 1 0 0.100 0.050 0.5 + 1926 67 1 1 1 0 0.100 0.050 0.5 + 1927 65 1 1 1 0 0.100 0.050 0.5 + 1928 58 1 1 1 0 0.100 0.050 0.5 + 1929 51 1 1 1 0 0.100 0.050 0.5 + 1930 46 1 1 1 0 0.100 0.050 0.5 + 1931 50 1 1 1 0 0.100 0.050 0.5 + 1932 60 1 1 1 0 0.100 0.050 0.5 + 1933 63 1 1 1 0 0.100 0.050 0.5 + 1934 62 1 1 1 0 0.100 0.050 0.5 + 1935 76 1 1 1 0 0.100 0.050 0.5 + 1936 71 1 1 1 0 0.100 0.050 0.5 + 1937 80 1 1 1 0 0.100 0.050 0.5 + 1938 88 1 1 1 0 0.100 0.050 0.5 + 1939 80 1 1 1 0 0.100 0.050 0.5 + 1940 81 1 1 1 0 0.100 0.050 0.5 + 1941 85 1 1 1 0 0.100 0.050 0.5 + 1942 90 1 1 1 0 0.100 0.050 0.5 + 1943 95 1 1 1 0 0.100 0.050 0.5 + 1944 110 1 1 1 0 0.100 0.050 0.5 + 1945 102 1 1 1 0 0.100 0.050 0.5 + 1946 101 1 1 1 0 0.100 0.050 0.5 + 1947 99 1 1 1 0 0.100 0.050 0.5 + 1948 99 1 1 1 0 0.100 0.050 0.5 + 1949 95 1 1 1 0 0.100 0.050 0.5 + 1950 95 1 1 1 0 0.100 0.050 0.5 + 1951 96 1 1 1 0 0.100 0.050 0.5 + 1952 110 1 1 1 0 0.100 0.050 0.5 + 1953 131 1 1 1 0 0.100 0.050 0.5 + 1954 133 1 1 1 0 0.100 0.050 0.5 + 1955 119 1 1 1 0 0.100 0.050 0.5 + 1956 129 1 1 1 0 0.100 0.050 0.5 + 1957 110 1 1 1 0 0.100 0.050 0.5 + 1958 121 1 1 1 0 0.100 0.050 0.5 + 1959 129 1 1 1 0 0.100 0.050 0.5 + 1960 132 1 1 1 0 0.100 0.050 0.5 + 1961 127 1 1 1 0 0.100 0.050 0.5 + 1962 115 1 1 1 0 0.100 0.050 0.5 + 1963 105 1 1 1 0 0.100 0.050 0.5 + 1964 100 1 1 1 0 0.100 0.050 0.5 + 1965 99 1 1 1 0 0.100 0.050 0.5 + 1966 100 1 1 1 0 0.100 0.050 0.5 + 1967 101 1 1 1 0 0.100 0.050 0.5 + 1968 103 1 1 1 0 0.100 0.050 0.5 + 1969 95 1 1 1 0 0.100 0.050 0.5 + 1970 91 1 1 1 0 0.100 0.050 0.5 + 1971 89 1 1 1 0 0.100 0.050 0.5 + 1972 78 1 1 1 0 0.100 0.050 0.5 + 1973 63 1 1 1 0 0.100 0.050 0.5 + 1974 61 1 1 1 0 0.100 0.050 0.5 + 1975 61 1 1 1 0 0.100 0.050 0.5 + 1976 55 1 1 1 0 0.100 0.050 0.5 + 1977 63 1 1 1 0 0.100 0.050 0.5 + 1978 71 1 1 1 0 0.100 0.050 0.5 + 1979 75 1 1 1 0 0.100 0.050 0.5 + 1980 94 1 1 1 0 0.100 0.050 0.5 + 1981 111 1 1 1 0 0.100 0.050 0.5 + 1982 127 1 1 1 0 0.100 0.050 0.5 + 1984 291 1 1 1 0 0.100 1.050 0.5 + 1985 357 1 1 1 0 0.036 0.050 0.5 + 1986 320 1 1 1 0 0.041 0.050 0.5 + 1987 321 1 1 1 0 0.040 0.050 0.5 + 1988 368 1 1 1 0 0.035 0.050 0.5 + 1989 358 1 1 1 0 0.025 0.050 0.5 + 1990 319 1 1 1 0 0.028 0.050 0.5 + 1991 318 1 1 1 0 0.038 0.050 0.5 + 1992 319 1 1 1 0 0.034 0.050 0.5 + 1993 373 1 1 1 0 0.099 0.050 0.5 + 1994 306 1 1 1 0 0.072 0.050 0.5 + 1995 330 1 1 1 0 0.036 0.050 0.5 + 1996 392 1 1 1 0 0.038 0.050 0.5 + 1997 405 1 1 1 0 0.025 0.050 0.5 + 1998 408 1 1 1 0 0.025 0.050 0.5 + 1999 393 1 1 1 0 0.023 0.050 0.5 + 2000 403 1 1 1 0 0.022 0.050 0.5 + 2001 362 1 1 1 0 0.041 0.050 0.5 + 2002 360 1 1 1 0 0.019 0.050 0.5 + 2003 329 1 1 1 0 0.018 0.050 0.5 + 2004 319 1 1 1 0 0.019 0.050 0.5 + 2005 296 1 1 1 0 0.017 0.050 0.5 + 2006 270 1 1 1 0 0.019 0.050 0.5 + 2007 252 1 1 1 0 0.020 0.050 0.5 + 2008 232 1 1 1 0 0.017 0.050 0.5 + 2009 222 1 1 1 0 0.018 0.050 0.5 + 2010 204 1 1 1 0 0.020 0.050 0.5 + 2011 198 1 1 1 0 0.015 0.050 0.5 + 2012 195 1 1 1 0 0.021 0.050 0.5 + 2013 179 1 1 1 0 0.017 0.050 0.5 + 2014 191 1 1 1 0 0.052 0.050 0.5 ## survey_data: 1977 1.4713 6 1 1 0 0.100 0.010 0.5 1978 1.1112 6 1 1 0 0.100 0.010 0.5 diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index 6268d647..e95a5e1c 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -1233,6 +1233,7 @@ DATA_SECTION init_ivector q_prior(1,nits); init_vector mu_log_q(1,nits); init_vector sd_log_q(1,nits); + init_ivector q_phz(1,nits); // |---------------------------------------------------------------------------------| // | Miscellaneous controls | @@ -1560,7 +1561,7 @@ PARAMETER_SECTION // | DEVIATIONS IN CATCHABILITY COEFFICIENTS ASSUMING A RANDOM WALK | // |---------------------------------------------------------------------------------| // | - init_bounded_matrix log_q_devs(1,nItNobs,1,n_it_nobs,-5.0,5.0,-1); + init_bounded_vector_vector log_q_devs(1,nItNobs,1,n_it_nobs,-5.1,5.0,q_phz); // |---------------------------------------------------------------------------------| // | CORRELATION COEFFICIENTS FOR AGE COMPOSITION DATA USED IN LOGISTIC NORMAL | @@ -1876,7 +1877,7 @@ FUNCTION void initParameters() steepness = theta(2); m = mfexp(theta(3)); rho = theta(6); - sigma_r = theta(7); + sigma_r = theta(7); //varphi = sqrt(1.0/theta(7)); sig = elem_prod(sqrt(rho) , varphi); //tau = elem_prod(sqrt(1.0-rho) , varphi); @@ -1966,23 +1967,7 @@ FUNCTION dvector cubic_spline(const dvector& spline_coffs, const dvector& la) //return(1.0*la); } -// FUNCTION dvar_matrix cubic_spline_matrix(const dvar_matrix& spline_coffs) -// { -// RETURN_ARRAYS_INCREMENT(); -// int nodes= spline_coffs.colmax()-spline_coffs.colmin()+1; -// int rmin = spline_coffs.rowmin(); -// int rmax = spline_coffs.rowmax(); - -// dvector ia(1,nodes); -// dvector fa(sage,nage); -// ia.fill_seqadd(0,1./(nodes-1)); -// //fa.fill_seqadd(sage,1); -// fa.fill_seqadd(0,1./(nage-sage)); -// vcubic_spline_function_array fna(rmin,rmax,ia,spline_coffs); -// RETURN_ARRAYS_DECREMENT(); -// return(fna(fa)); - -// } + /** @@ -2767,6 +2752,7 @@ FUNCTION calcTotalCatch fishery that would remove potential spawn that would not be surveyed. - d3_survey_data: (iyr index(it) gear area group sex wt timing) - for MLE of survey q, using weighted mean of zt to calculate q. + - March 30, 2015. Added deviation in q for random walk. TODO list: [] - add capability to accomodate priors for survey q's. @@ -2779,11 +2765,12 @@ FUNCTION calcSurveyObservations int ii,kk,ig,nz; double di; - dvariable ftmp; + dvariable ftmp,zbar; dvar_vector Na(sage,nage); dvar_vector va(sage,nage); dvar_vector sa(sage,nage); epsilon.initialize(); + qt.initialize(); xi.initialize(); it_hat.initialize(); @@ -2845,42 +2832,87 @@ FUNCTION calcSurveyObservations dvar_vector t1 = rowsum(V); - dvar_vector zt = log(it) - log(t1(iz,nz)); - dvariable zbar = sum(elem_prod(zt,wt)); - q(kk) = mfexp(zbar); - - // | survey residuals - epsilon(kk).sub(iz,nz) = elem_div(zt - zbar,it_log_se(kk)(iz,nz)); - it_hat(kk).sub(iz,nz) = q(kk) * t1(iz,nz); + dvar_vector zt = log(it); + - // | SPECIAL CASE: penalized random walk in q process error only. - if( q_prior(kk)==2 ) + // | March 30, 2015. Issue #37. + // | Added switch for 3 q_prior(kk) options. + // | q_prior(kk) = 0 = constant fixed Q + // | q_prior(kk) = 1 = constant MLE Q + // | q_prior(kk) = 2 = penalized random walk in Q + switch( q_prior(kk) ) { - epsilon(kk).initialize(); - dvar_vector fd_zt = first_difference(zt); - dvariable zw_bar = sum(elem_prod(fd_zt,wt(iz,nz-1))); - epsilon(kk).sub(iz,nz-1) = elem_div(fd_zt - zw_bar,it_log_se(kk)(iz,nz-1)); - qt(kk)(iz) = exp(zt(iz)); - for(ii=iz+1;ii<=nz;ii++) - { - qt(kk)(ii) = qt(kk)(ii-1) * exp(fd_zt(ii-1)); - } - it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); + case 0: // Constant fixed Q + q(kk) = exp(mu_log_q(kk)); + it_hat(kk).sub(iz,nz) = q(kk) * t1(iz,nz); + zt -= log(it_hat(kk).sub(iz,nz)); + break; + + case 1: // Constant MLE Q + zt -= log(t1(iz,nz)); + zbar = sum(elem_prod(zt,wt)); + q(kk) = mfexp(zbar); + + // | survey residuals + it_hat(kk).sub(iz,nz) = q(kk) * t1(iz,nz); + zt -= zbar; + //epsilon(kk).sub(iz,nz) = elem_div(zt,it_log_se(kk)(iz,nz)); + break; + + case 2: // Penalized random walk in Q + zt -= log(t1(iz,nz)); + qt(kk)(iz) = exp( zt(iz) + log_q_devs(kk)(iz) ); + for(ii=iz+1; ii<=nz; ii++) + { + COUT(ii); + qt(kk)(ii) = qt(kk)(ii-1) * exp(log_q_devs(kk)(ii)); + } + it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); + zt -= log(qt(kk)(iz,nz)); + //epsilon(kk).sub(iz,nz)= elem_div(zt,it_log_se(kk)(iz,nz)); + break; } - // | MIXED ERROR MODEL for random walk in q - if( q_prior(kk)==3 ) + // Standardized observation error residuals. + COUT(qt); + epsilon(kk).sub(iz,nz) = elem_div(zt,it_log_se(kk)(iz,nz)); + + // Standardized process error residuals. + if(active(log_q_devs(kk))) { - dvar_vector proerr = zt - zbar; - qt(kk)(ii) = exp(zbar + proerr(iz)); - for(ii=iz+1;ii<=nz;ii++) - { - proerr(ii) = zt(ii) - zt(ii-1); - qt(kk)(ii) = qt(kk)(ii-1) * exp(proerr(ii)); - } - xi(kk).sub(iz,nz) = elem_div(proerr,it_log_pe(kk)(iz,nz)); - it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); + dvar_vector fd_qt = first_difference( log(qt(kk)(iz,nz)) ); + xi(kk).sub(iz,nz-1) = elem_div(fd_qt,it_log_pe(kk)(iz,nz-1)); } + +// TO BE DEPRECATED +// // | SPECIAL CASE: penalized random walk in q process error only. +// if( q_prior(kk)==2 ) +// { +// epsilon(kk).initialize(); +// dvar_vector fd_zt = first_difference(zt); +// dvariable zw_bar = sum(elem_prod(fd_zt,wt(iz,nz-1))); +// epsilon(kk).sub(iz,nz-1) = elem_div(fd_zt - zw_bar,it_log_se(kk)(iz,nz-1)); +// qt(kk)(iz) = exp(zt(iz)); +// for(ii=iz+1;ii<=nz;ii++) +// { +// qt(kk)(ii) = qt(kk)(ii-1) * exp(fd_zt(ii-1)); +// } +// it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); +// } +// +// // | MIXED ERROR MODEL for random walk in q +// if( q_prior(kk)==3 ) +// { +// dvar_vector proerr = zt - zbar; +// qt(kk)(ii) = exp(zbar + proerr(iz)); +// for(ii=iz+1;ii<=nz;ii++) +// { +// proerr(ii) = zt(ii) - zt(ii-1); +// qt(kk)(ii) = qt(kk)(ii-1) * exp(proerr(ii)); +// } +// xi(kk).sub(iz,nz) = elem_div(proerr,it_log_pe(kk)(iz,nz)); +// it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); +// } } if(verbose)cout<<"**** Ok after calcSurveyObservations ****"< Date: Mon, 30 Mar 2015 10:25:03 -0700 Subject: [PATCH 03/19] q_dev parameters now estimable --- fba/HALIBUT/DATA/Halibut2014.dat | 58 ++++++++++++++++---------------- src/admb-code/iscam.tpl | 8 ++--- 2 files changed, 32 insertions(+), 34 deletions(-) diff --git a/fba/HALIBUT/DATA/Halibut2014.dat b/fba/HALIBUT/DATA/Halibut2014.dat index c52a7dbb..e831c7da 100644 --- a/fba/HALIBUT/DATA/Halibut2014.dat +++ b/fba/HALIBUT/DATA/Halibut2014.dat @@ -458,7 +458,7 @@ 1980 94 1 1 1 0 0.100 0.050 0.5 1981 111 1 1 1 0 0.100 0.050 0.5 1982 127 1 1 1 0 0.100 0.050 0.5 - 1984 291 1 1 1 0 0.100 1.050 0.5 + 1984 291 1 1 1 0 0.100 1.000 0.5 1985 357 1 1 1 0 0.036 0.050 0.5 1986 320 1 1 1 0 0.041 0.050 0.5 1987 321 1 1 1 0 0.040 0.050 0.5 @@ -490,34 +490,34 @@ 2013 179 1 1 1 0 0.017 0.050 0.5 2014 191 1 1 1 0 0.052 0.050 0.5 ## survey_data: - 1977 1.4713 6 1 1 0 0.100 0.010 0.5 - 1978 1.1112 6 1 1 0 0.100 0.010 0.5 - 1980 2.0090 6 1 1 0 0.100 0.010 0.5 - 1981 2.6670 6 1 1 0 0.100 0.010 0.5 - 1982 2.8739 6 1 1 0 0.100 0.010 0.5 - 1983 2.8816 6 1 1 0 0.100 0.010 0.5 - 1984 9.3014 6 1 1 0 0.100 0.010 0.5 - 1985 8.9419 6 1 1 0 0.100 0.010 0.5 - 1986 6.2568 6 1 1 0 0.100 0.010 0.5 - 1996 12.8858 6 1 1 0 0.100 0.010 0.5 - 1997 7.7841 6 1 1 0 0.034 0.010 0.5 - 1998 6.9727 6 1 1 0 0.034 0.010 0.5 - 1999 6.2649 6 1 1 0 0.034 0.010 0.5 - 2000 6.6265 6 1 1 0 0.034 0.010 0.5 - 2001 6.3842 6 1 1 0 0.028 0.010 0.5 - 2002 6.3893 6 1 1 0 0.026 0.010 0.5 - 2003 5.7876 6 1 1 0 0.027 0.010 0.5 - 2004 6.4579 6 1 1 0 0.026 0.010 0.5 - 2005 5.9282 6 1 1 0 0.027 0.010 0.5 - 2006 5.4306 6 1 1 0 0.027 0.010 0.5 - 2007 5.7267 6 1 1 0 0.027 0.010 0.5 - 2008 5.6536 6 1 1 0 0.026 0.010 0.5 - 2009 5.4952 6 1 1 0 0.027 0.010 0.5 - 2010 5.1239 6 1 1 0 0.029 0.010 0.5 - 2011 5.0435 6 1 1 0 0.029 0.010 0.5 - 2012 5.5606 6 1 1 0 0.030 0.010 0.5 - 2013 4.7301 6 1 1 0 0.034 0.010 0.5 - 2014 5.1645 6 1 1 0 0.032 0.010 0.5 + 1977 1.4713 6 1 1 0 0.100 0.050 0.5 + 1978 1.1112 6 1 1 0 0.100 0.050 0.5 + 1980 2.0090 6 1 1 0 0.100 0.050 0.5 + 1981 2.6670 6 1 1 0 0.100 0.050 0.5 + 1982 2.8739 6 1 1 0 0.100 1.000 0.5 + 1983 2.8816 6 1 1 0 0.100 0.050 0.5 + 1984 9.3014 6 1 1 0 0.100 1.000 0.5 + 1985 8.9419 6 1 1 0 0.100 0.050 0.5 + 1986 6.2568 6 1 1 0 0.100 0.050 0.5 + 1996 12.8858 6 1 1 0 0.100 0.050 0.5 + 1997 7.7841 6 1 1 0 0.034 1.000 0.5 + 1998 6.9727 6 1 1 0 0.034 0.050 0.5 + 1999 6.2649 6 1 1 0 0.034 0.050 0.5 + 2000 6.6265 6 1 1 0 0.034 0.050 0.5 + 2001 6.3842 6 1 1 0 0.028 0.050 0.5 + 2002 6.3893 6 1 1 0 0.026 0.050 0.5 + 2003 5.7876 6 1 1 0 0.027 0.050 0.5 + 2004 6.4579 6 1 1 0 0.026 0.050 0.5 + 2005 5.9282 6 1 1 0 0.027 0.050 0.5 + 2006 5.4306 6 1 1 0 0.027 0.050 0.5 + 2007 5.7267 6 1 1 0 0.027 0.050 0.5 + 2008 5.6536 6 1 1 0 0.026 0.050 0.5 + 2009 5.4952 6 1 1 0 0.027 0.050 0.5 + 2010 5.1239 6 1 1 0 0.029 0.050 0.5 + 2011 5.0435 6 1 1 0 0.029 0.050 0.5 + 2012 5.5606 6 1 1 0 0.030 0.050 0.5 + 2013 4.7301 6 1 1 0 0.034 0.100 0.5 + 2014 5.1645 6 1 1 0 0.032 0.100 0.5 ## ## ------------------------------------------------------------------------- ## ## AGE COMPOSITION DATA (ROW YEAR, COL=AGE) Ragged object ## diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index e95a5e1c..d63697e1 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -1561,7 +1561,7 @@ PARAMETER_SECTION // | DEVIATIONS IN CATCHABILITY COEFFICIENTS ASSUMING A RANDOM WALK | // |---------------------------------------------------------------------------------| // | - init_bounded_vector_vector log_q_devs(1,nItNobs,1,n_it_nobs,-5.1,5.0,q_phz); + init_bounded_vector_vector log_q_devs(1,nItNobs,1,n_it_nobs,-5.0,5.0,q_phz); // |---------------------------------------------------------------------------------| // | CORRELATION COEFFICIENTS FOR AGE COMPOSITION DATA USED IN LOGISTIC NORMAL | @@ -2864,7 +2864,6 @@ FUNCTION calcSurveyObservations qt(kk)(iz) = exp( zt(iz) + log_q_devs(kk)(iz) ); for(ii=iz+1; ii<=nz; ii++) { - COUT(ii); qt(kk)(ii) = qt(kk)(ii-1) * exp(log_q_devs(kk)(ii)); } it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); @@ -2874,13 +2873,12 @@ FUNCTION calcSurveyObservations } // Standardized observation error residuals. - COUT(qt); epsilon(kk).sub(iz,nz) = elem_div(zt,it_log_se(kk)(iz,nz)); // Standardized process error residuals. if(active(log_q_devs(kk))) { - dvar_vector fd_qt = first_difference( log(qt(kk)(iz,nz)) ); + dvar_vector fd_qt = first_difference( log_q_devs(kk)(iz,nz) ); xi(kk).sub(iz,nz-1) = elem_div(fd_qt,it_log_pe(kk)(iz,nz-1)); } @@ -3139,7 +3137,7 @@ FUNCTION calcObjectiveFunction // } // nlvec(2,k)=dnorm(epsilon(k),sig_it); nlvec(2,k)=dnorm(epsilon(k),1.0); - nlvec(2,k)=dnorm(xi(k),1.0); + nlvec(2,k)+=dnorm(xi(k),1.0); } // |---------------------------------------------------------------------------------| From 9cef2267940db6879099bd9126eda5f87d9aff5d Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Mon, 30 Mar 2015 10:28:14 -0700 Subject: [PATCH 04/19] Extended nlvec to accomodate process errors in catchability --- src/admb-code/iscam.tpl | 39 ++++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index d63697e1..0466a5d8 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -1656,7 +1656,7 @@ PARAMETER_SECTION // | - delta -> residuals between estimated R and R from S-R curve (process err) // | matrix log_rt(1,n_ag,syr-nage+sage,nyr); - matrix nlvec(1,7,1,ilvec); + matrix nlvec(1,8,1,ilvec); matrix epsilon(1,nItNobs,1,n_it_nobs); matrix xi(1,nItNobs,1,n_it_nobs); matrix it_hat(1,nItNobs,1,n_it_nobs); @@ -3079,10 +3079,11 @@ FUNCTION calcObjectiveFunction Likelihoods (nlvec): -1) likelihood of the catch data -2) likelihood of the survey abundance index - -3) likelihood of age composition data - -4) likelihood for stock-recruitment relationship - -5) penalized likelihood for fishery selectivities + -3) likelihood component for random walk in q. + -4) likelihood of age composition data + -5) likelihood for stock-recruitment relationship -6) penalized likelihood for fishery selectivities + -7) penalized likelihood for fishery selectivities TODO list: @@ -3137,7 +3138,7 @@ FUNCTION calcObjectiveFunction // } // nlvec(2,k)=dnorm(epsilon(k),sig_it); nlvec(2,k)=dnorm(epsilon(k),1.0); - nlvec(2,k)+=dnorm(xi(k),1.0); + nlvec(3,k)=dnorm(xi(k),1.0); } // |---------------------------------------------------------------------------------| @@ -3196,19 +3197,19 @@ FUNCTION calcObjectiveFunction switch( int(nCompLikelihood(k)) ) { case 1: - nlvec(3,k) = dmvlogistic(O,P,nu,age_tau2(k),dMinP(k)); + nlvec(4,k) = dmvlogistic(O,P,nu,age_tau2(k),dMinP(k)); break; case 2: - nlvec(3,k) = dmultinom(O,P,nu,age_tau2(k),dMinP(k)); + nlvec(4,k) = dmultinom(O,P,nu,age_tau2(k),dMinP(k)); break; case 3: if( !active(log_age_tau2(k)) ) // LN1 Model { - nlvec(3,k) = cLN_Age(); + nlvec(4,k) = cLN_Age(); } else { - nlvec(3,k) = cLN_Age( exp(log_age_tau2(k)) ); + nlvec(4,k) = cLN_Age( exp(log_age_tau2(k)) ); } // Residual @@ -3223,11 +3224,11 @@ FUNCTION calcObjectiveFunction //logistic_normal cLN_Age( O,P,dMinP(k),dEps(k) ); if( active(phi1(k)) && !active(phi2(k)) ) // LN2 Model { - nlvec(3,k) = cLN_Age(exp(log_age_tau2(k)),phi1(k)); + nlvec(4,k) = cLN_Age(exp(log_age_tau2(k)),phi1(k)); } if( active(phi1(k)) && active(phi2(k)) ) // LN3 Model { - nlvec(3,k) = cLN_Age(exp(log_age_tau2(k)),phi1(k),phi2(k)); + nlvec(4,k) = cLN_Age(exp(log_age_tau2(k)),phi1(k),phi2(k)); } // Residual @@ -3242,11 +3243,11 @@ FUNCTION calcObjectiveFunction case 5: // Logistic-normal with student-t if( !active(log_degrees_of_freedom(k)) ) { - nlvec(3,k) = cLST_Age(); + nlvec(4,k) = cLST_Age(); } else { - nlvec(3,k) = cLST_Age(exp(log_degrees_of_freedom(k))); + nlvec(4,k) = cLST_Age(exp(log_degrees_of_freedom(k))); } // Residual @@ -3257,10 +3258,10 @@ FUNCTION calcObjectiveFunction } break; case 6: // Multinomial with estimated effective sample size. - nlvec(3,k) = mult_likelihood(O,P,nu,log_degrees_of_freedom(k)); + nlvec(4,k) = mult_likelihood(O,P,nu,log_degrees_of_freedom(k)); break; case 7: // Multivariate-t - nlvec(3,k) = multivariate_t_likelihood(O,P,log_age_tau2(k), + nlvec(4,k) = multivariate_t_likelihood(O,P,log_age_tau2(k), log_degrees_of_freedom(k), phi1(k),nu); age_tau2(k) = exp(value(log_age_tau2(k))); @@ -3292,7 +3293,7 @@ FUNCTION calcObjectiveFunction { for(g=1;g<=ngroup;g++) { - nlvec(4,g) = dnorm(delta(g),sigma_r(g)); + nlvec(5,g) = dnorm(delta(g),sigma_r(g)); } } @@ -3322,12 +3323,12 @@ FUNCTION calcObjectiveFunction { //curvature in selectivity parameters dvar_vector df2 = first_difference(first_difference(log_sel(k)(ig)(i))); - nlvec(5,k) += lambda_1(k)/(nage-sage+1)*df2*df2; + nlvec(6,k) += lambda_1(k)/(nage-sage+1)*df2*df2; //penalty for dome-shapeness for(j=sage;j<=nage-1;j++) if(log_sel(k,ig,i,j)>log_sel(k,ig,i,j+1)) - nlvec(6,k)+=lambda_2(k) + nlvec(7,k)+=lambda_2(k) *square( log_sel(k,ig,i,j)-log_sel(k,ig,i,j+1) ); } } @@ -3348,7 +3349,7 @@ FUNCTION calcObjectiveFunction for(j=sage;j<=nage;j++) { dvar_vector df2 = first_difference(first_difference(trans_log_sel(j))); - nlvec(7,k) += lambda_3(k)/(nage-sage+1)*norm2(df2); + nlvec(8,k) += lambda_3(k)/(nage-sage+1)*norm2(df2); } } } From dad1c07035d17d66474c61ee9f3e88c718e1fb46 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Mon, 30 Mar 2015 10:33:11 -0700 Subject: [PATCH 05/19] Time varying q, close #37 --- README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/README.md b/README.md index a9e55825..618bc7a4 100644 --- a/README.md +++ b/README.md @@ -279,3 +279,6 @@ _____________________________________________________________ Last changed on: Source code: https://github.com/smartell/iSCAM _____________________________________________________________ + +## Code Changes + - March 30, 2015. Added time varying catchbility options (random walk). From 5dd148fe23ce747951e5fecf63fe54728771ef31 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Mon, 30 Mar 2015 10:50:42 -0700 Subject: [PATCH 06/19] Fixed array out of bounds error --- fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl | 8 ++-- fba/HALIBUT/DATA/Halibut2014.dat | 50 +++++++++++------------ src/admb-code/iscam.tpl | 31 ++++---------- 3 files changed, 38 insertions(+), 51 deletions(-) diff --git a/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl b/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl index 3bb8bcc6..7507c5ff 100644 --- a/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl +++ b/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl @@ -113,10 +113,10 @@ ## 2 = RANDOM WALK Q (use prior mean & sd for penalized random walk) ## ———————————————————————————————————————————————————————————————————————————————————— ## 2 # -number of surveys (n_it_nobs) -2 2 # -QTYPE (see legend above) -0 0 # -prior log(mean) -0 0 # -prior sd (set to 0 for uniformative prior) -2 2 # -Estimation phase + 1 1 # -QTYPE (see legend above) + 0 0 # -prior log(mean) + 0 0 # -prior sd (set to 0 for uniformative prior) +-2 -2 # -Estimation phase ## ———————————————————————————————————————————————————————————————————————————————————— ## ## diff --git a/fba/HALIBUT/DATA/Halibut2014.dat b/fba/HALIBUT/DATA/Halibut2014.dat index e831c7da..9715bca7 100644 --- a/fba/HALIBUT/DATA/Halibut2014.dat +++ b/fba/HALIBUT/DATA/Halibut2014.dat @@ -490,34 +490,34 @@ 2013 179 1 1 1 0 0.017 0.050 0.5 2014 191 1 1 1 0 0.052 0.050 0.5 ## survey_data: - 1977 1.4713 6 1 1 0 0.100 0.050 0.5 - 1978 1.1112 6 1 1 0 0.100 0.050 0.5 - 1980 2.0090 6 1 1 0 0.100 0.050 0.5 - 1981 2.6670 6 1 1 0 0.100 0.050 0.5 + 1977 1.4713 6 1 1 0 0.100 0.010 0.5 + 1978 1.1112 6 1 1 0 0.100 0.010 0.5 + 1980 2.0090 6 1 1 0 0.100 0.010 0.5 + 1981 2.6670 6 1 1 0 0.100 0.010 0.5 1982 2.8739 6 1 1 0 0.100 1.000 0.5 - 1983 2.8816 6 1 1 0 0.100 0.050 0.5 + 1983 2.8816 6 1 1 0 0.100 0.010 0.5 1984 9.3014 6 1 1 0 0.100 1.000 0.5 - 1985 8.9419 6 1 1 0 0.100 0.050 0.5 - 1986 6.2568 6 1 1 0 0.100 0.050 0.5 - 1996 12.8858 6 1 1 0 0.100 0.050 0.5 + 1985 8.9419 6 1 1 0 0.100 0.010 0.5 + 1986 6.2568 6 1 1 0 0.100 0.010 0.5 + 1996 12.8858 6 1 1 0 0.100 0.010 0.5 1997 7.7841 6 1 1 0 0.034 1.000 0.5 - 1998 6.9727 6 1 1 0 0.034 0.050 0.5 - 1999 6.2649 6 1 1 0 0.034 0.050 0.5 - 2000 6.6265 6 1 1 0 0.034 0.050 0.5 - 2001 6.3842 6 1 1 0 0.028 0.050 0.5 - 2002 6.3893 6 1 1 0 0.026 0.050 0.5 - 2003 5.7876 6 1 1 0 0.027 0.050 0.5 - 2004 6.4579 6 1 1 0 0.026 0.050 0.5 - 2005 5.9282 6 1 1 0 0.027 0.050 0.5 - 2006 5.4306 6 1 1 0 0.027 0.050 0.5 - 2007 5.7267 6 1 1 0 0.027 0.050 0.5 - 2008 5.6536 6 1 1 0 0.026 0.050 0.5 - 2009 5.4952 6 1 1 0 0.027 0.050 0.5 - 2010 5.1239 6 1 1 0 0.029 0.050 0.5 - 2011 5.0435 6 1 1 0 0.029 0.050 0.5 - 2012 5.5606 6 1 1 0 0.030 0.050 0.5 - 2013 4.7301 6 1 1 0 0.034 0.100 0.5 - 2014 5.1645 6 1 1 0 0.032 0.100 0.5 + 1998 6.9727 6 1 1 0 0.034 0.010 0.5 + 1999 6.2649 6 1 1 0 0.034 0.010 0.5 + 2000 6.6265 6 1 1 0 0.034 0.010 0.5 + 2001 6.3842 6 1 1 0 0.028 0.010 0.5 + 2002 6.3893 6 1 1 0 0.026 0.010 0.5 + 2003 5.7876 6 1 1 0 0.027 0.010 0.5 + 2004 6.4579 6 1 1 0 0.026 0.010 0.5 + 2005 5.9282 6 1 1 0 0.027 0.010 0.5 + 2006 5.4306 6 1 1 0 0.027 0.010 0.5 + 2007 5.7267 6 1 1 0 0.027 0.010 0.5 + 2008 5.6536 6 1 1 0 0.026 0.010 0.5 + 2009 5.4952 6 1 1 0 0.027 0.010 0.5 + 2010 5.1239 6 1 1 0 0.029 0.010 0.5 + 2011 5.0435 6 1 1 0 0.029 0.010 0.5 + 2012 5.5606 6 1 1 0 0.030 0.010 0.5 + 2013 4.7301 6 1 1 0 0.034 0.010 0.5 + 2014 5.1645 6 1 1 0 0.032 0.010 0.5 ## ## ------------------------------------------------------------------------- ## ## AGE COMPOSITION DATA (ROW YEAR, COL=AGE) Ragged object ## diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index 0466a5d8..66a6aef1 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -1292,12 +1292,13 @@ DATA_SECTION // | ilvec[2] -> number of surveys (nItNobs) // | ilvec[3] -> number of age-compisition data sets (nAgears) // | ilvec[4] -> container for recruitment deviations. - ivector ilvec(1,7); + ivector ilvec(1,8); !! ilvec = ngear; !! ilvec(1) = 1; !! ilvec(2) = nItNobs; - !! ilvec(3) = nAgears; - !! ilvec(4) = ngroup; + !! ilvec(3) = nItNobs; + !! ilvec(4) = nAgears; + !! ilvec(5) = ngroup; // |---------------------------------------------------------------------------------| @@ -1409,29 +1410,16 @@ DATA_SECTION // | MANAGEMENT STRATEGY EVALUATION INPUTS // |---------------------------------------------------------------------------------| // | - - - //LOC_CALCS - // ifstream ifile("Halibut2012.mse"); - // if(ifile) - // { - // cout<<"Vader is happy"< Date: Mon, 30 Mar 2015 13:06:18 -0700 Subject: [PATCH 07/19] Control file change --- fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl b/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl index 7507c5ff..b92d97ab 100644 --- a/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl +++ b/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl @@ -113,10 +113,10 @@ ## 2 = RANDOM WALK Q (use prior mean & sd for penalized random walk) ## ———————————————————————————————————————————————————————————————————————————————————— ## 2 # -number of surveys (n_it_nobs) - 1 1 # -QTYPE (see legend above) + 2 2 # -QTYPE (see legend above) 0 0 # -prior log(mean) 0 0 # -prior sd (set to 0 for uniformative prior) --2 -2 # -Estimation phase + 1 1 # -Estimation phase ## ———————————————————————————————————————————————————————————————————————————————————— ## ## From ad781e430dc81ce46b362fabffabdb9786de0b46 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Mon, 30 Mar 2015 14:14:56 -0700 Subject: [PATCH 08/19] added compositionLikelihoods.hpp --- .../include/compositionLikelihoods.hpp | 37 +++++++++++++++++++ src/admb-code/include/lib_iscam.h | 1 + src/admb-code/iscam.tpl | 9 +++-- 3 files changed, 44 insertions(+), 3 deletions(-) create mode 100644 src/admb-code/include/compositionLikelihoods.hpp diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp new file mode 100644 index 00000000..d3f2ed77 --- /dev/null +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -0,0 +1,37 @@ +// compositionLikelihoods.hpp +#include + +#ifndef COMPOSITIONLIKELIHOODS_H +#define COMPOSITIONLIKELIHOODS_H + +/** + * @file compositionLikelihoods.hpp + * @defgroup Likelihoods + * @author Steven Martell + * @namespace acl + * @date Feb 10, 2014 + * @title Selectivity functions + * @details Uses abstract base class for computing negative loglikelihoods + */ +namespace acl +{ + + template + class compositionLikelihoods + { + private: + T m_x; + public: + virtual const T nloglike(const T & x) const = 0; + virtual const T residual(const T & x) const = 0; + virtual ~compositionLikelihoods(){}; + + void set_x(T & x) { this -> m_x = x; } + T get_x() const{ return m_x; } + + }; + + +} + +#endif \ No newline at end of file diff --git a/src/admb-code/include/lib_iscam.h b/src/admb-code/include/lib_iscam.h index 650b961a..600579d5 100644 --- a/src/admb-code/include/lib_iscam.h +++ b/src/admb-code/include/lib_iscam.h @@ -8,6 +8,7 @@ #include "growthModels.hpp" #include "multinomial.h" #include "selex.hpp" +#include "compositionLikelihoods.hpp" #endif \ No newline at end of file diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index 66a6aef1..64736891 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -2558,9 +2558,9 @@ FUNCTION calcComposition dvar_vector pred_ca = ca * age_age(e); A_hat(kk)(ii) = pred_ca(n_A_sage(kk),n_A_nage(kk)); if( n_A_nage(kk) < nage ) - { - A_hat(kk)(ii)(n_A_nage(kk)) += sum( pred_ca(n_A_nage(kk)+1,nage) ); - } + { + A_hat(kk)(ii)(n_A_nage(kk)) += sum( pred_ca(n_A_nage(kk)+1,nage) ); + } } else { @@ -3167,6 +3167,9 @@ FUNCTION calcObjectiveFunction { O(ii) = d3_A_obs(k)(i).sub(n_A_sage(k),n_A_nage(k)); P(ii) = A_hat(k)(i).sub(n_A_sage(k),n_A_nage(k)); + COUT(O(ii)); + COUT(P(ii)); + exit(1); ii ++; } //if( iyr <= nyr ) naa++; From f262f639f89411434fe5617a126ecc8a61088c97 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Mon, 30 Mar 2015 15:29:06 -0700 Subject: [PATCH 09/19] Working on compositionLikelihoods.hpp --- .../include/compositionLikelihoods.hpp | 58 +++++++++++++++++-- 1 file changed, 53 insertions(+), 5 deletions(-) diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index d3f2ed77..86f6e2ec 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -20,18 +20,66 @@ namespace acl class compositionLikelihoods { private: - T m_x; + int r1,r2; /// index for first and last row. + int c1,c2; /// index for rist and last column. + ivector m_jmin; /// index for ragged start columns. + ivector m_jmax; /// index for ragged end columns. + + T m_O; /// observed composition object + public: - virtual const T nloglike(const T & x) const = 0; - virtual const T residual(const T & x) const = 0; + virtual const T nloglike(const T & _O) const = 0; + virtual const T residual(const T & _O) const = 0; virtual ~compositionLikelihoods(){}; - void set_x(T & x) { this -> m_x = x; } - T get_x() const{ return m_x; } + void set_O(T & O) { this -> m_O = O; } + T get_O() const{ return m_O; } + void tail_compression_indexes(); }; + void acl::test(); + + + + // |--------------------------------------------------------------------------------| + // | MULTIVARIATE LOGISTIC NEGATIVE LOGLIKELIHOOD derived class | + // |--------------------------------------------------------------------------------| + template + inline + const T dmvlogistic(const T& O, const T& P) + { + + } + + /** + * @brief Negative loglikelihood for compostion data using multivaritae logistic + * @author Steve Martell + * @param T usually a dmatrix + */ + template + class multivariteLogistic: public compositionLikelihoods + { + private: + T m_P; /// predicted composition object. + + public: + multivariteLogistic(T _P): m_P(_P) {} + + const T nloglike(const T &O) const + { + return acl::dmvlogistic(O, this->get_P()); + } + + T get_P() const {return m_P; } + void set_P(T P) {this->m_P=P;} + + }; + } + + + #endif \ No newline at end of file From 4e0ef7ace34d19292cefebc294cea1e5162c3b79 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Mon, 30 Mar 2015 15:40:43 -0700 Subject: [PATCH 10/19] Still working out methods for abstract base class --- .../include/compositionLikelihoods.hpp | 31 +++++++++++++++++-- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index 86f6e2ec..13047e44 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -16,6 +16,12 @@ namespace acl { + /** + * @brief Base class for composition likelihoods. + * @details Virtual methods for negative loglikelihood and standardized residuals. + * + * @tparam T [description] + */ template class compositionLikelihoods { @@ -28,17 +34,35 @@ namespace acl T m_O; /// observed composition object public: + // constructors + compositionLikelihoods(){} + compositionLikelihoods(const T& _O) + :m_O(_O) + { + r1 = m_O.rowmin(); + r2 = m_O.rowmax(); + c1 = m_O.colmin(); + c2 = m_O.colmax(); + } + + // virtual methods + virtual ~compositionLikelihoods(){}; virtual const T nloglike(const T & _O) const = 0; virtual const T residual(const T & _O) const = 0; - virtual ~compositionLikelihoods(){}; + // getters & setters void set_O(T & O) { this -> m_O = O; } T get_O() const{ return m_O; } + // methods void tail_compression_indexes(); + }; - void acl::test(); + // void comtail_compression_indexes() + // { + + // } @@ -64,7 +88,8 @@ namespace acl T m_P; /// predicted composition object. public: - multivariteLogistic(T _P): m_P(_P) {} + multivariteLogistic(T &_O, T &_P) + :compositionLikelihoods(_O),m_P(_P) {} const T nloglike(const T &O) const { From a1113a872c02386ae018e18f92314425d39e6e77 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Tue, 31 Mar 2015 09:36:21 -0700 Subject: [PATCH 11/19] Creating DATA and DVAR typenames for composition likelihoods --- .../include/compositionLikelihoods.hpp | 78 ++++++++++++++++--- src/admb-code/iscam.tpl | 14 +++- 2 files changed, 77 insertions(+), 15 deletions(-) diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index 13047e44..d00ff35e 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -16,6 +16,53 @@ namespace acl { + template + inline + const T tailCompression(const T &_M, const double pmin=0.0) + { + int r1,r2,c1,c2; + r1 = _M.rowmin(); + r2 = _M.rowmax(); + c1 = _M.colmin(); + c2 = _M.colmax(); + + ivector jmin(r1,r2); + ivector jmax(r1,r2); + for(int i = r1; i <= r2; i++ ) + { + auto o = _M(i); + + jmin(i) = c1+1; // index for min column + jmax(i) = c2; // index for max column + auto cumsum = o/sum(o); + for(int j = c1+1; j <= c2; j++ ) + { + cumsum(j) += cumsum(j-1); + cumsum(j) <= pmin ? jmin(i)++ : NULL; + j != c2 ? 1.0 - cumsum(j) < pmin ? jmax(i)-- : NULL : NULL; + } + } + + // Now compress the matrix. + T M_; + T M = _M; + M_.allocate(r1,r2,jmin,jmax); + M_.initialize(); + + // fill ragged array M_ + for(int i = r1; i <= r2; i++ ) + { + M(i) /= sum(M(i)); + M_(i)(jmin(i),jmax(i)) = M(i)(jmin(i),jmax(i)); + + // add cumulative sum to tails. + M_(i)(jmin(i)) = sum(M(i)(c1,jmin(i))); + M_(i)(jmax(i)) = sum(M(i)(jmax(i),c2)); + } + return M_; + + } + /** * @brief Base class for composition likelihoods. * @details Virtual methods for negative loglikelihood and standardized residuals. @@ -53,19 +100,11 @@ namespace acl // getters & setters void set_O(T & O) { this -> m_O = O; } T get_O() const{ return m_O; } - - // methods - void tail_compression_indexes(); - + }; - // void comtail_compression_indexes() - // { - - // } - // |--------------------------------------------------------------------------------| // | MULTIVARIATE LOGISTIC NEGATIVE LOGLIKELIHOOD derived class | // |--------------------------------------------------------------------------------| @@ -73,7 +112,16 @@ namespace acl inline const T dmvlogistic(const T& O, const T& P) { + T nll; + return nll; + } + template + inline + const T dmvlogisticResidual(const T& O, const T& P) + { + T nu; + return nu; } /** @@ -81,20 +129,26 @@ namespace acl * @author Steve Martell * @param T usually a dmatrix */ - template + template class multivariteLogistic: public compositionLikelihoods { private: T m_P; /// predicted composition object. public: - multivariteLogistic(T &_O, T &_P) - :compositionLikelihoods(_O),m_P(_P) {} + multivariteLogistic(const DATA &_O, const DVAR &_P) + :compositionLikelihoods(_O),m_P(_P) {} const T nloglike(const T &O) const { + T rO = tailCompression(O); return acl::dmvlogistic(O, this->get_P()); } + + const T residual(const T &O) const + { + return acl::dmvlogisticResidual(O, this->get_P()); + } T get_P() const {return m_P; } void set_P(T P) {this->m_P=P;} diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index 64736891..dccf4b96 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -3144,6 +3144,15 @@ FUNCTION calcObjectiveFunction // | TODO: // | [ ] - change A_nu to data-type variable, does not need to be differentiable. // | [ ] - issue 29. Fix submatrix O, P for prospective analysis & sex/area/group. + + // Testing new abstract compositionLikelihood class. + + acl::compositionLikelihoods *ptr_AgeCompLike; + + + // delete *ptr_AgeComeLike; + + A_nu.initialize(); for(k=1;k<=nAgears;k++) { @@ -3167,14 +3176,13 @@ FUNCTION calcObjectiveFunction { O(ii) = d3_A_obs(k)(i).sub(n_A_sage(k),n_A_nage(k)); P(ii) = A_hat(k)(i).sub(n_A_sage(k),n_A_nage(k)); - COUT(O(ii)); - COUT(P(ii)); - exit(1); ii ++; } //if( iyr <= nyr ) naa++; //if( iyr < syr ) iaa++; } + ptr_AgeCompLike = new acl::multivariteLogistic(O,P); + dvar_matrix ell = ptr_AgeCompLike->nloglike(P); //dmatrix O = trans(trans(d3_A_obs(k)).sub(n_A_sage(k),n_A_nage(k))).sub(iaa,naa); //dvar_matrix P = trans(trans(A_hat(k)).sub(n_A_sage(k),n_A_nage(k))).sub(iaa,naa); From 72e61cdb58278e1b845c982b2693952a2ccede92 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Tue, 31 Mar 2015 11:25:49 -0700 Subject: [PATCH 12/19] Still developing negLogLikelihood abstract class --- .../include/compositionLikelihoods.hpp | 220 +++++++++++------- src/admb-code/iscam.tpl | 6 +- 2 files changed, 145 insertions(+), 81 deletions(-) diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index d00ff35e..d6946a35 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -1,11 +1,11 @@ -// compositionLikelihoods.hpp +// negLogLikelihood.hpp #include -#ifndef COMPOSITIONLIKELIHOODS_H -#define COMPOSITIONLIKELIHOODS_H +#ifndef negLogLikelihood_H +#define negLogLikelihood_H /** - * @file compositionLikelihoods.hpp + * @file negLogLikelihood.hpp * @defgroup Likelihoods * @author Steven Martell * @namespace acl @@ -16,61 +16,62 @@ namespace acl { - template - inline - const T tailCompression(const T &_M, const double pmin=0.0) - { - int r1,r2,c1,c2; - r1 = _M.rowmin(); - r2 = _M.rowmax(); - c1 = _M.colmin(); - c2 = _M.colmax(); - - ivector jmin(r1,r2); - ivector jmax(r1,r2); - for(int i = r1; i <= r2; i++ ) - { - auto o = _M(i); - - jmin(i) = c1+1; // index for min column - jmax(i) = c2; // index for max column - auto cumsum = o/sum(o); - for(int j = c1+1; j <= c2; j++ ) - { - cumsum(j) += cumsum(j-1); - cumsum(j) <= pmin ? jmin(i)++ : NULL; - j != c2 ? 1.0 - cumsum(j) < pmin ? jmax(i)-- : NULL : NULL; - } - } - - // Now compress the matrix. - T M_; - T M = _M; - M_.allocate(r1,r2,jmin,jmax); - M_.initialize(); - - // fill ragged array M_ - for(int i = r1; i <= r2; i++ ) - { - M(i) /= sum(M(i)); - M_(i)(jmin(i),jmax(i)) = M(i)(jmin(i),jmax(i)); - - // add cumulative sum to tails. - M_(i)(jmin(i)) = sum(M(i)(c1,jmin(i))); - M_(i)(jmax(i)) = sum(M(i)(jmax(i),c2)); - } - return M_; - - } + // template + // inline + // const T tailCompression(const T &_M, const double pmin=0.0) + // { + // int r1,r2,c1,c2; + // r1 = _M.rowmin(); + // r2 = _M.rowmax(); + // c1 = _M.colmin(); + // c2 = _M.colmax(); + + // ivector jmin(r1,r2); + // ivector jmax(r1,r2); + // for(int i = r1; i <= r2; i++ ) + // { + // auto o = _M(i); + + // jmin(i) = c1+1; // index for min column + // jmax(i) = c2; // index for max column + // auto cumsum = o/sum(o); + // for(int j = c1+1; j <= c2; j++ ) + // { + // cumsum(j) += cumsum(j-1); + // cumsum(j) <= pmin ? jmin(i)++ : NULL; + // j != c2 ? 1.0 - cumsum(j) < pmin ? jmax(i)-- : NULL : NULL; + // } + // } + + // // Now compress the matrix. + // T M_; + // T M = _M; + // M_.allocate(r1,r2,jmin,jmax); + // M_.initialize(); + + // // fill ragged array M_ + // for(int i = r1; i <= r2; i++ ) + // { + // M(i) /= sum(M(i)); + // M_(i)(jmin(i),jmax(i)) = M(i)(jmin(i),jmax(i)); + + // // add cumulative sum to tails. + // M_(i)(jmin(i)) = sum(M(i)(c1,jmin(i))); + // M_(i)(jmax(i)) = sum(M(i)(jmax(i),c2)); + // } + // return M_; + + // } /** * @brief Base class for composition likelihoods. * @details Virtual methods for negative loglikelihood and standardized residuals. * - * @tparam T [description] + * @tparam DATA Data type argument. + * @tparam DVAR Variable type argument. */ - template - class compositionLikelihoods + template + class negLogLikelihood { private: int r1,r2; /// index for first and last row. @@ -78,49 +79,111 @@ namespace acl ivector m_jmin; /// index for ragged start columns. ivector m_jmax; /// index for ragged end columns. - T m_O; /// observed composition object + DATA m_O; /// observed composition object public: // constructors - compositionLikelihoods(){} - compositionLikelihoods(const T& _O) + negLogLikelihood(){} + negLogLikelihood(const DATA& _O) :m_O(_O) { r1 = m_O.rowmin(); r2 = m_O.rowmax(); c1 = m_O.colmin(); c2 = m_O.colmax(); + + getRaggedVectors(); } - // virtual methods - virtual ~compositionLikelihoods(){}; - virtual const T nloglike(const T & _O) const = 0; - virtual const T residual(const T & _O) const = 0; + // pure virtual methods + virtual ~negLogLikelihood(){}; + virtual const DVAR nloglike(const DVAR & _P) const = 0; + virtual const DVAR residual(const DVAR & _P) const = 0; + // virtual methods + inline + const DATA tailCompression(const DATA &_M) const; + + // methods + inline + void getRaggedVectors(double pmin = 0.0); + // getters & setters - void set_O(T & O) { this -> m_O = O; } - T get_O() const{ return m_O; } + void set_O(DATA & O) { this -> m_O = O; } + DATA get_O() const{ return m_O; } }; + template + inline + const DATA acl::negLogLikelihood::tailCompression(const DATA &_M) const + { + DATA R; + DATA M; + R.allocate(r1,r2,m_jmin,m_jmax); + R.initialize(); + // fill ragged array R + for(int i = r1; i <= r2; i++ ) + { + M(i) /= sum(M(i)); + R(i)(m_jmin(i),m_jmax(i)) = M(i)(m_jmin(i),m_jmax(i)); + + // add cumulative sum to tails. + R(i)(m_jmin(i)) = sum(M(i)(c1,m_jmin(i))); + R(i)(m_jmax(i)) = sum(M(i)(m_jmax(i),c2)); + } + return R; + } + + /** + * @brief Get column indicies for each row to compress matrix into ragged array. + * @details Determine the start and end positions of each row in which to compute + * the likelihoods where values < pmin are pooled into the tails of the distribution. + * + * @param pmin value to assume for the minium proportion. + * @return NULL Modifies private member variables: m_jmin and m_jmax. + */ + template + inline + void acl::negLogLikelihood::getRaggedVectors(double pmin) + { + m_jmin.allocate(r1,r2); + m_jmax.allocate(r1,r2); + + for(int i = r1; i <= r2; i++ ) + { + dvector o = m_O(i); + + m_jmin(i) = c1+1; // index for min column + m_jmax(i) = c2; // index for max column + dvector cumsum = o/sum(o); + for(int j = c1+1; j <= c2; j++ ) + { + cumsum(j) += cumsum(j-1); + cumsum(j) <= pmin ? m_jmin(i)++ : NULL; + j != c2 ? 1.0 - cumsum(j) < pmin ? m_jmax(i)-- : NULL : NULL; + } + } + + } // |--------------------------------------------------------------------------------| // | MULTIVARIATE LOGISTIC NEGATIVE LOGLIKELIHOOD derived class | // |--------------------------------------------------------------------------------| - template + template inline - const T dmvlogistic(const T& O, const T& P) + const DVAR dmvlogistic(const DATA& O, const DVAR& P) { - T nll; + DVAR nll; return nll; } - template + template inline - const T dmvlogisticResidual(const T& O, const T& P) + const DVAR dmvlogisticResidual(const DATA& O, const DVAR& P) { - T nu; + DVAR nu; return nu; } @@ -130,28 +193,29 @@ namespace acl * @param T usually a dmatrix */ template - class multivariteLogistic: public compositionLikelihoods + class multivariteLogistic: public negLogLikelihood { private: - T m_P; /// predicted composition object. + DVAR m_P; /// predicted composition object. public: + // constructor multivariteLogistic(const DATA &_O, const DVAR &_P) - :compositionLikelihoods(_O),m_P(_P) {} + :negLogLikelihood(_O),m_P(_P) {} - const T nloglike(const T &O) const + const DVAR nloglike(const DVAR &P) const { - T rO = tailCompression(O); - return acl::dmvlogistic(O, this->get_P()); + //DATA rO = tailCompression(this->get_O); + return acl::dmvlogistic(this->get_O(), this->get_P()); } - const T residual(const T &O) const + const DVAR residual(const DVAR &P) const { - return acl::dmvlogisticResidual(O, this->get_P()); + return acl::dmvlogisticResidual(this->get_O(), this->get_P()); } - T get_P() const {return m_P; } - void set_P(T P) {this->m_P=P;} + DVAR get_P() const {return m_P; } + void set_P(DVAR P) {this->m_P=P;} }; diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index dccf4b96..4154e837 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -3147,7 +3147,7 @@ FUNCTION calcObjectiveFunction // Testing new abstract compositionLikelihood class. - acl::compositionLikelihoods *ptr_AgeCompLike; + acl::negLogLikelihood *ptr_AgeCompLike; // delete *ptr_AgeComeLike; @@ -3181,8 +3181,8 @@ FUNCTION calcObjectiveFunction //if( iyr <= nyr ) naa++; //if( iyr < syr ) iaa++; } - ptr_AgeCompLike = new acl::multivariteLogistic(O,P); - dvar_matrix ell = ptr_AgeCompLike->nloglike(P); + ptr_AgeCompLike = new acl::multivariteLogistic(O,P); + //dvar_matrix ell = ptr_AgeCompLike->nloglike(P); //dmatrix O = trans(trans(d3_A_obs(k)).sub(n_A_sage(k),n_A_nage(k))).sub(iaa,naa); //dvar_matrix P = trans(trans(A_hat(k)).sub(n_A_sage(k),n_A_nage(k))).sub(iaa,naa); From 6602b5577bdbc1d559cd076362fd481a4f4e1697 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Tue, 31 Mar 2015 21:16:07 -0700 Subject: [PATCH 13/19] Fixed minor error in getRaggedVectors --- scripts/Template.ctl | 27 +++-- .../include/compositionLikelihoods.hpp | 113 ++++++++++++------ src/admb-code/iscam.tpl | 4 +- 3 files changed, 92 insertions(+), 52 deletions(-) diff --git a/scripts/Template.ctl b/scripts/Template.ctl index bb28e38d..f3ff2643 100644 --- a/scripts/Template.ctl +++ b/scripts/Template.ctl @@ -100,20 +100,23 @@ ## ———————————————————————————————————————————————————————————————————————————————————— ## ## ## -## ------------------------------------------------------------------------- ## -## PRIORS FOR SURVEY Q ## -## Prior type: ## -## 0 - uninformative prior ## -## 1 - normal prior density for log(q) ## -## 2 - random walk in q ## -## ------------------------------------------------------------------------- ## -1 # -number of surveys (nits) -0 # -prior type (see legend above) -0 # -prior log(mean) -0 # -prior sd -## ------------------------------------------------------------------------- ## +## ———————————————————————————————————————————————————————————————————————————————————— ## +## ABUNDANCE OBSERVATION MODELS +## ———————————————————————————————————————————————————————————————————————————————————— ## +## QTYPE: +## 0 = FIXED SURVEY Q (specify log(mean) for prior log(mean)) +## 1 = CONSTANT Q (use MLE for q and optional informative prior) +## 2 = RANDOM WALK Q (use prior mean & sd for penalized random walk) +## ———————————————————————————————————————————————————————————————————————————————————— ## +1 # -number of surveys (n_it_nobs) + 2 # -QTYPE (see legend above) + 0 # -prior log(mean) + 0 # -prior sd (set to 0 for uniformative prior) + 1 # -Estimation phase +## ———————————————————————————————————————————————————————————————————————————————————— ## ## + ## ------------------------------------------------------------------------- ## ## OTHER MISCELANEOUS CONTROLS ## ## ------------------------------------------------------------------------- ## diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index d6946a35..fa78f89a 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -80,12 +80,13 @@ namespace acl ivector m_jmax; /// index for ragged end columns. DATA m_O; /// observed composition object + DVAR m_P; /// predicted composition object public: // constructors negLogLikelihood(){} - negLogLikelihood(const DATA& _O) - :m_O(_O) + negLogLikelihood(const DATA& _O, const DVAR& _P) + :m_O(_O), m_P(_P) { r1 = m_O.rowmin(); r2 = m_O.rowmax(); @@ -97,12 +98,14 @@ namespace acl // pure virtual methods virtual ~negLogLikelihood(){}; - virtual const DVAR nloglike(const DVAR & _P) const = 0; - virtual const DVAR residual(const DVAR & _P) const = 0; + virtual const dvariable nloglike() const = 0; + virtual const DVAR residual() const = 0; // virtual methods inline const DATA tailCompression(const DATA &_M) const; + // inline + // const DVAR tailCompression(const DVAR &_M) const; // methods inline @@ -111,29 +114,33 @@ namespace acl // getters & setters void set_O(DATA & O) { this -> m_O = O; } DATA get_O() const{ return m_O; } - - }; - template - inline - const DATA acl::negLogLikelihood::tailCompression(const DATA &_M) const - { - DATA R; - DATA M; - R.allocate(r1,r2,m_jmin,m_jmax); - R.initialize(); - // fill ragged array R - for(int i = r1; i <= r2; i++ ) + DVAR get_P() const {return m_P; } + void set_P(DVAR &P) {this->m_P=P;} + + template + inline + const T compress(const T& _M) const { - M(i) /= sum(M(i)); - R(i)(m_jmin(i),m_jmax(i)) = M(i)(m_jmin(i),m_jmax(i)); + T R; + T M = _M; + R.allocate(r1,r2,m_jmin,m_jmax); + R.initialize(); + // fill ragged array R + for(int i = r1; i <= r2; i++ ) + { + M(i) /= sum(M(i)); + R(i)(m_jmin(i),m_jmax(i)) = M(i)(m_jmin(i),m_jmax(i)); - // add cumulative sum to tails. - R(i)(m_jmin(i)) = sum(M(i)(c1,m_jmin(i))); - R(i)(m_jmax(i)) = sum(M(i)(m_jmax(i),c2)); + // add cumulative sum to tails. + R(i)(m_jmin(i)) = sum(M(i)(c1,m_jmin(i))); + R(i)(m_jmax(i)) = sum(M(i)(m_jmax(i),c2)); + } + + return R; } - return R; - } + + }; /** * @brief Get column indicies for each row to compress matrix into ragged array. @@ -154,28 +161,35 @@ namespace acl { dvector o = m_O(i); - m_jmin(i) = c1+1; // index for min column + m_jmin(i) = c1; // index for min column m_jmax(i) = c2; // index for max column dvector cumsum = o/sum(o); - for(int j = c1+1; j <= c2; j++ ) + for(int j = c1; j < c2; j++ ) { - cumsum(j) += cumsum(j-1); cumsum(j) <= pmin ? m_jmin(i)++ : NULL; - j != c2 ? 1.0 - cumsum(j) < pmin ? m_jmax(i)-- : NULL : NULL; + j != c2 ? 1.0 - cumsum(j) <= pmin ? m_jmax(i)-- : NULL : NULL; + + cumsum(j+1) += cumsum(j); } } - } // |--------------------------------------------------------------------------------| // | MULTIVARIATE LOGISTIC NEGATIVE LOGLIKELIHOOD derived class | // |--------------------------------------------------------------------------------| - template + template inline - const DVAR dmvlogistic(const DATA& O, const DVAR& P) + const T dmvlogistic(const DATA& O, const DVAR& P) { - DVAR nll; + T nll; + // nll.allocate(P); + // int i; + // for( i = O.rowmin(); i<= O.rowmax(); i++) + // { + // nll(i) = log(O(i)) - log(P(i)); + // nll(i) -= mean(nll(i)); + // } return nll; } @@ -184,6 +198,13 @@ namespace acl const DVAR dmvlogisticResidual(const DATA& O, const DVAR& P) { DVAR nu; + nu.allocate(P); + int i; + for( i = O.rowmin(); i<= O.rowmax(); i++) + { + nu(i) = log(O(i)) - log(P(i)); + nu(i) -= mean(nu(i)); + } return nu; } @@ -196,26 +217,40 @@ namespace acl class multivariteLogistic: public negLogLikelihood { private: - DVAR m_P; /// predicted composition object. + DVAR m_rO; /// ragged observed composition object. + DVAR m_rP; /// ragged predicted composition object. + DVAR m_nu; /// logistic residuals. public: // constructor + // todo: add eps value to constructor. multivariteLogistic(const DATA &_O, const DVAR &_P) - :negLogLikelihood(_O),m_P(_P) {} + :negLogLikelihood(_O,_P) + { + DATA tmp = this->compress(this->get_O()); + cout<compress(_P); + // m_nu = acl::dmvlogisticResidual(m_rO,m_rP); + } - const DVAR nloglike(const DVAR &P) const + + const dvariable nloglike() const { - //DATA rO = tailCompression(this->get_O); - return acl::dmvlogistic(this->get_O(), this->get_P()); + // DATA rO = this->compress(this->get_O()); + // DVAR rP = this->compress(this->get_P()); + // m_nu = acl::dmvlogisticResidual(rO,rP); + // return acl::dmvlogistic(rO, rP); + return 0; } - const DVAR residual(const DVAR &P) const + const DVAR residual() const { return acl::dmvlogisticResidual(this->get_O(), this->get_P()); } - DVAR get_P() const {return m_P; } - void set_P(DVAR P) {this->m_P=P;} + + DVAR get_nu() const {return m_nu; } + void set_nu(DVAR &R){this->m_nu=R;} }; diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index 4154e837..590ad0ef 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -3182,7 +3182,9 @@ FUNCTION calcObjectiveFunction //if( iyr < syr ) iaa++; } ptr_AgeCompLike = new acl::multivariteLogistic(O,P); - //dvar_matrix ell = ptr_AgeCompLike->nloglike(P); + dvariable ell = ptr_AgeCompLike->nloglike(); + COUT(ell); + //dmatrix O = trans(trans(d3_A_obs(k)).sub(n_A_sage(k),n_A_nage(k))).sub(iaa,naa); //dvar_matrix P = trans(trans(A_hat(k)).sub(n_A_sage(k),n_A_nage(k))).sub(iaa,naa); From 9fc0b0243f83df2d9266ddd7565e66d801197738 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Tue, 31 Mar 2015 22:01:37 -0700 Subject: [PATCH 14/19] More progress on the multivariateLogistic derived class --- .../include/compositionLikelihoods.hpp | 54 ++++++++++--------- src/admb-code/iscam.tpl | 2 +- 2 files changed, 31 insertions(+), 25 deletions(-) diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index fa78f89a..e8df0747 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -178,19 +178,20 @@ namespace acl // |--------------------------------------------------------------------------------| // | MULTIVARIATE LOGISTIC NEGATIVE LOGLIKELIHOOD derived class | // |--------------------------------------------------------------------------------| - template + /** + * @brief multivariate logistic negative log likelihood + * @details The negative log likeihood for the multivariate + * logistic distribution based on the MLE of the variance. + * + * @param NU Matrix of residuals + * @return the negative loglikelihood + */ + template inline - const T dmvlogistic(const DATA& O, const DVAR& P) + const T dmvlogistic(const DVAR& NU) { - T nll; - // nll.allocate(P); - // int i; - // for( i = O.rowmin(); i<= O.rowmax(); i++) - // { - // nll(i) = log(O(i)) - log(P(i)); - // nll(i) -= mean(nll(i)); - // } - return nll; + T var = 1.0/size_count(NU)*norm2(NU); + return size_count(NU) * log(var); } template @@ -213,34 +214,35 @@ namespace acl * @author Steve Martell * @param T usually a dmatrix */ - template + template class multivariteLogistic: public negLogLikelihood { private: - DVAR m_rO; /// ragged observed composition object. + DATA m_rO; /// ragged observed composition object. DVAR m_rP; /// ragged predicted composition object. DVAR m_nu; /// logistic residuals. public: // constructor // todo: add eps value to constructor. - multivariteLogistic(const DATA &_O, const DVAR &_P) + multivariteLogistic(const DATA &_O, const DVAR &_P, const double eps=1.e-8) :negLogLikelihood(_O,_P) { + // tail compression DATA tmp = this->compress(this->get_O()); - cout<compress(_P); - // m_nu = acl::dmvlogisticResidual(m_rO,m_rP); + set_rO(tmp); + DVAR vmp = this->compress(this->get_P()); + set_rP(vmp); + + // residuals + DVAR tnu = acl::dmvlogisticResidual(this->get_rO(),this->get_rP()); + set_nu(tnu); } - const dvariable nloglike() const + const T nloglike() const { - // DATA rO = this->compress(this->get_O()); - // DVAR rP = this->compress(this->get_P()); - // m_nu = acl::dmvlogisticResidual(rO,rP); - // return acl::dmvlogistic(rO, rP); - return 0; + return acl::dmvlogistic(this->get_nu()); } const DVAR residual() const @@ -248,8 +250,12 @@ namespace acl return acl::dmvlogisticResidual(this->get_O(), this->get_P()); } - + DATA get_rO() const {return m_rO; } + DVAR get_rP() const {return m_rP; } DVAR get_nu() const {return m_nu; } + + void set_rO(DATA &X){this->m_rO=X;} + void set_rP(DVAR &X){this->m_rP=X;} void set_nu(DVAR &R){this->m_nu=R;} }; diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index 590ad0ef..564d8948 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -3181,7 +3181,7 @@ FUNCTION calcObjectiveFunction //if( iyr <= nyr ) naa++; //if( iyr < syr ) iaa++; } - ptr_AgeCompLike = new acl::multivariteLogistic(O,P); + ptr_AgeCompLike = new acl::multivariteLogistic(O,P); dvariable ell = ptr_AgeCompLike->nloglike(); COUT(ell); From c1843f361cf48f200a93fc9bf346717d523afddd Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Tue, 31 Mar 2015 22:03:11 -0700 Subject: [PATCH 15/19] More progress on the multivariateLogistic derived class --- src/admb-code/include/compositionLikelihoods.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index e8df0747..2a881be9 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -230,9 +230,10 @@ namespace acl { // tail compression DATA tmp = this->compress(this->get_O()); - set_rO(tmp); + set_rO(tmp+eps); DVAR vmp = this->compress(this->get_P()); - set_rP(vmp); + set_rP(vmp+eps); + // residuals DVAR tnu = acl::dmvlogisticResidual(this->get_rO(),this->get_rP()); From 689b828a0babf3237320e7fd98cf025922955d82 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Tue, 31 Mar 2015 22:04:45 -0700 Subject: [PATCH 16/19] Now have an actual likelihood from the derived class multivariateLogistic --- src/admb-code/include/compositionLikelihoods.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index 2a881be9..a486b7e2 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -229,10 +229,10 @@ namespace acl :negLogLikelihood(_O,_P) { // tail compression - DATA tmp = this->compress(this->get_O()); - set_rO(tmp+eps); - DVAR vmp = this->compress(this->get_P()); - set_rP(vmp+eps); + DATA tmp = this->compress(this->get_O()) + eps; + set_rO(tmp); + DVAR vmp = this->compress(this->get_P()) + eps; + set_rP(vmp); // residuals From f2c99bcc78308c34bee82398739e27c33c84b8b3 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Tue, 31 Mar 2015 22:36:20 -0700 Subject: [PATCH 17/19] removed old code --- .../include/compositionLikelihoods.hpp | 49 +------------------ src/admb-code/iscam.tpl | 3 +- 2 files changed, 3 insertions(+), 49 deletions(-) diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index a486b7e2..b229f820 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -16,53 +16,6 @@ namespace acl { - // template - // inline - // const T tailCompression(const T &_M, const double pmin=0.0) - // { - // int r1,r2,c1,c2; - // r1 = _M.rowmin(); - // r2 = _M.rowmax(); - // c1 = _M.colmin(); - // c2 = _M.colmax(); - - // ivector jmin(r1,r2); - // ivector jmax(r1,r2); - // for(int i = r1; i <= r2; i++ ) - // { - // auto o = _M(i); - - // jmin(i) = c1+1; // index for min column - // jmax(i) = c2; // index for max column - // auto cumsum = o/sum(o); - // for(int j = c1+1; j <= c2; j++ ) - // { - // cumsum(j) += cumsum(j-1); - // cumsum(j) <= pmin ? jmin(i)++ : NULL; - // j != c2 ? 1.0 - cumsum(j) < pmin ? jmax(i)-- : NULL : NULL; - // } - // } - - // // Now compress the matrix. - // T M_; - // T M = _M; - // M_.allocate(r1,r2,jmin,jmax); - // M_.initialize(); - - // // fill ragged array M_ - // for(int i = r1; i <= r2; i++ ) - // { - // M(i) /= sum(M(i)); - // M_(i)(jmin(i),jmax(i)) = M(i)(jmin(i),jmax(i)); - - // // add cumulative sum to tails. - // M_(i)(jmin(i)) = sum(M(i)(c1,jmin(i))); - // M_(i)(jmax(i)) = sum(M(i)(jmax(i),c2)); - // } - // return M_; - - // } - /** * @brief Base class for composition likelihoods. * @details Virtual methods for negative loglikelihood and standardized residuals. @@ -248,7 +201,7 @@ namespace acl const DVAR residual() const { - return acl::dmvlogisticResidual(this->get_O(), this->get_P()); + return this->get_nu(); } DATA get_rO() const {return m_rO; } diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index 564d8948..64afce75 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -3183,7 +3183,7 @@ FUNCTION calcObjectiveFunction } ptr_AgeCompLike = new acl::multivariteLogistic(O,P); dvariable ell = ptr_AgeCompLike->nloglike(); - COUT(ell); + //COUT(ell); //dmatrix O = trans(trans(d3_A_obs(k)).sub(n_A_sage(k),n_A_nage(k))).sub(iaa,naa); @@ -3198,6 +3198,7 @@ FUNCTION calcObjectiveFunction { case 1: nlvec(4,k) = dmvlogistic(O,P,nu,age_tau2(k),dMinP(k)); + cout<<"like: "< Date: Wed, 1 Apr 2015 14:22:17 -0700 Subject: [PATCH 18/19] Have aggregation working for dmvlogistic and abstract base class, results are not the same however --- .../include/compositionLikelihoods.hpp | 98 ++++++++++++++++--- src/admb-code/iscam.tpl | 2 +- 2 files changed, 84 insertions(+), 16 deletions(-) diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index b229f820..20645429 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -31,9 +31,14 @@ namespace acl int c1,c2; /// index for rist and last column. ivector m_jmin; /// index for ragged start columns. ivector m_jmax; /// index for ragged end columns. + imatrix m_jagg; /// index for aggregation among bins. DATA m_O; /// observed composition object + DATA m_rO; /// observed ragged composition object + DVAR m_P; /// predicted composition object + DVAR m_rP; /// predicted ragged composition object + public: // constructors @@ -57,13 +62,14 @@ namespace acl // virtual methods inline const DATA tailCompression(const DATA &_M) const; - // inline - // const DVAR tailCompression(const DVAR &_M) const; - + // methods inline void getRaggedVectors(double pmin = 0.0); + inline + void aggregate(const double pmin = 0.0); + // getters & setters void set_O(DATA & O) { this -> m_O = O; } DATA get_O() const{ return m_O; } @@ -71,6 +77,13 @@ namespace acl DVAR get_P() const {return m_P; } void set_P(DVAR &P) {this->m_P=P;} + DATA get_rO() const {return m_rO; } + void set_rO(DATA &X){this->m_rO=X;} + + DVAR get_rP() const {return m_rP; } + void set_rP(DVAR &X){this->m_rP=X;} + + template inline const T compress(const T& _M) const @@ -95,6 +108,44 @@ namespace acl }; + template + inline + void acl::negLogLikelihood::aggregate(const double pmin) + { + m_jmin.allocate(r1,r2); + m_jmax.allocate(r1,r2); + ivector n(r1,r2); + n.initialize(); + + // get number of observations > pmin each year. + for(int i = r1; i <= r2; i++ ) + { + dvector oo = m_O(i)/sum(m_O(i)); + for(int j = c1; j <= c2; j++ ) + { + if( oo(j) > pmin ) n(i)++; + } + } + + m_rO.allocate(r1,r2,1,n); + m_rP.allocate(r1,r2,1,n); + m_rO.initialize(); + m_rP.initialize(); + + for(int i = r1; i <= r2; i++ ) + { + int k = 1; + dvector oo = m_O(i)/sum(m_O(i)); + dvar_vector pp = m_P(i)/sum(m_P(i)); + for(int j = c1; j <= c2; j++ ) + { + m_rO(i)(k) += oo(j); + m_rP(i)(k) += pp(j); + if( oo(j)>pmin && k pmin ? n(i)++ : NULL; } } + } + + + + + + + + + + + + + + // |--------------------------------------------------------------------------------| // | MULTIVARIATE LOGISTIC NEGATIVE LOGLIKELIHOOD derived class | // |--------------------------------------------------------------------------------| + /** * @brief multivariate logistic negative log likelihood * @details The negative log likeihood for the multivariate @@ -178,19 +247,22 @@ namespace acl public: // constructor // todo: add eps value to constructor. - multivariteLogistic(const DATA &_O, const DVAR &_P, const double eps=1.e-8) + multivariteLogistic(const DATA &_O, const DVAR &_P, const double eps=1.e-3) :negLogLikelihood(_O,_P) { // tail compression DATA tmp = this->compress(this->get_O()) + eps; - set_rO(tmp); + this->set_rO(tmp); DVAR vmp = this->compress(this->get_P()) + eps; - set_rP(vmp); + this->set_rP(vmp); + // residuals + this->aggregate(eps); DVAR tnu = acl::dmvlogisticResidual(this->get_rO(),this->get_rP()); set_nu(tnu); + } @@ -204,12 +276,8 @@ namespace acl return this->get_nu(); } - DATA get_rO() const {return m_rO; } - DVAR get_rP() const {return m_rP; } DVAR get_nu() const {return m_nu; } - void set_rO(DATA &X){this->m_rO=X;} - void set_rP(DVAR &X){this->m_rP=X;} void set_nu(DVAR &R){this->m_nu=R;} }; diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index 64afce75..04566ac8 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -3181,7 +3181,7 @@ FUNCTION calcObjectiveFunction //if( iyr <= nyr ) naa++; //if( iyr < syr ) iaa++; } - ptr_AgeCompLike = new acl::multivariteLogistic(O,P); + ptr_AgeCompLike = new acl::multivariteLogistic(O,P,dMinP(k)); dvariable ell = ptr_AgeCompLike->nloglike(); //COUT(ell); From b8a9d8af818ee53b641172cc6b3ad5ffc2c38cf0 Mon Sep 17 00:00:00 2001 From: Steve Martell Date: Wed, 1 Apr 2015 14:27:02 -0700 Subject: [PATCH 19/19] Now have multivariateLogistic producing the same results with aggregation method --- src/admb-code/include/compositionLikelihoods.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp index 20645429..185a7e82 100644 --- a/src/admb-code/include/compositionLikelihoods.hpp +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -212,8 +212,9 @@ namespace acl inline const T dmvlogistic(const DVAR& NU) { - T var = 1.0/size_count(NU)*norm2(NU); - return size_count(NU) * log(var); + int n = size_count(NU) - (NU.rowmax()-NU.rowmin()+1); + T var = 1.0/n * norm2(NU); + return n * log(var); } template