Hi there, I'm relatively inexperienced in mata and I'm struggeling to optimize a function.
In short, I want to model hypothetical microdata (n=1,000) that underlie a given average for a number of years (1990-2020). I'm studying coverage of a medication program, so that starts at 0, increases, and eventually gets to a maximum. I want to optimize 3 parameters (b0_1, b0_2, b0_3 below), one related to the starting year, one related to the year the maximum is reached, and one related to the number of years between starting and reaching maximum coverage.
i have attached my code (you can run this all the way)
for now I can work with this and fill in values for b0_1, b0_2 and b0_3 by hand. But how do I optimize these three parameters (and minimize f by doing so)?
many thanks
Joost
In short, I want to model hypothetical microdata (n=1,000) that underlie a given average for a number of years (1990-2020). I'm studying coverage of a medication program, so that starts at 0, increases, and eventually gets to a maximum. I want to optimize 3 parameters (b0_1, b0_2, b0_3 below), one related to the starting year, one related to the year the maximum is reached, and one related to the number of years between starting and reaching maximum coverage.
i have attached my code (you can run this all the way)
for now I can work with this and fill in values for b0_1, b0_2 and b0_3 by hand. But how do I optimize these three parameters (and minimize f by doing so)?
many thanks
Joost
Code:
* start here
set more off
clear
set obs 1000
* make random set for prevalence figures (P)
gen uniform=runiform()
* generate the logit of P
gen logitP=invnormal(uniform)
gen uniform1=runiform()
gen uniform2=runiform()
gen uniform3=runiform()
mata
mata clear
logitP=st_data(.,"logitP")
uniform1=st_data(.,"uniform1")
uniform2=st_data(.,"uniform2")
uniform3=st_data(.,"uniform3")
* choose slope parameters
b1=0.3
b2=-0.3
b3=-0.1
logitcovmax=b1*logitP+1*invnormal(uniform1)
yrsuntilmaxcov=exp(b2*logitP+0.5*invnormal(uniform2))
yrsuntilPCT=exp(b3*logitP+0.5*invnormal(uniform3))
x1990=J(rows(logitcovmax),1,1990)
x1991=J(rows(logitcovmax),1,1991)
x1992=J(rows(logitcovmax),1,1992)
x1993=J(rows(logitcovmax),1,1993)
x1994=J(rows(logitcovmax),1,1994)
x1995=J(rows(logitcovmax),1,1995)
x1996=J(rows(logitcovmax),1,1996)
x1997=J(rows(logitcovmax),1,1997)
x1998=J(rows(logitcovmax),1,1998)
x1999=J(rows(logitcovmax),1,1999)
x2000=J(rows(logitcovmax),1,2000)
x2001=J(rows(logitcovmax),1,2001)
x2002=J(rows(logitcovmax),1,2002)
x2003=J(rows(logitcovmax),1,2003)
x2004=J(rows(logitcovmax),1,2004)
x2005=J(rows(logitcovmax),1,2005)
x2006=J(rows(logitcovmax),1,2006)
x2007=J(rows(logitcovmax),1,2007)
x2008=J(rows(logitcovmax),1,2008)
x2009=J(rows(logitcovmax),1,2009)
x2010=J(rows(logitcovmax),1,2010)
x2011=J(rows(logitcovmax),1,2011)
x2012=J(rows(logitcovmax),1,2012)
x2013=J(rows(logitcovmax),1,2013)
x2014=J(rows(logitcovmax),1,2014)
x2015=J(rows(logitcovmax),1,2015)
x2016=J(rows(logitcovmax),1,2016)
x2017=J(rows(logitcovmax),1,2017)
x2018=J(rows(logitcovmax),1,2018)
x2019=J(rows(logitcovmax),1,2019)
x2020=J(rows(logitcovmax),1,2020)
* set intercept parameters (FOR NOW, these I want to optimize)
b0_1=2.5
b0_2=0.5
b0_3=1.5
for(i=1; i<=rows(logitcovmax); i++)
{
if (x1990[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x1990[i]=0
if (x1990[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x1990[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x1990[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x1990[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x1990[i]=(x1990[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x1991[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x1991[i]=0
if (x1991[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x1991[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x1991[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x1991[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x1991[i]=(x1991[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x1992[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x1992[i]=0
if (x1992[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x1992[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x1992[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x1992[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x1992[i]=(x1992[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x1993[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x1993[i]=0
if (x1993[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x1993[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x1993[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x1993[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x1993[i]=(x1993[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x1994[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x1994[i]=0
if (x1994[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x1994[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x1994[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x1994[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x1994[i]=(x1994[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x1995[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x1995[i]=0
if (x1995[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x1995[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x1995[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x1995[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x1995[i]=(x1995[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x1996[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x1996[i]=0
if (x1996[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x1996[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x1996[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x1996[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x1996[i]=(x1996[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x1997[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x1997[i]=0
if (x1997[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x1997[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x1997[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x1997[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x1997[i]=(x1997[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x1998[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x1998[i]=0
if (x1998[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x1998[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x1998[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x1998[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x1998[i]=(x1998[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x1999[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x1999[i]=0
if (x1999[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x1999[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x1999[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x1999[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x1999[i]=(x1999[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2000[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2000[i]=0
if (x2000[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2000[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2000[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2000[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2000[i]=(x2000[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2001[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2001[i]=0
if (x2001[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2001[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2001[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2001[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2001[i]=(x2001[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2002[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2002[i]=0
if (x2002[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2002[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2002[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2002[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2002[i]=(x2002[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2003[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2003[i]=0
if (x2003[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2003[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2003[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2003[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2003[i]=(x2003[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2004[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2004[i]=0
if (x2004[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2004[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2004[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2004[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2004[i]=(x2004[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2005[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2005[i]=0
if (x2005[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2005[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2005[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2005[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2005[i]=(x2005[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2006[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2006[i]=0
if (x2006[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2006[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2006[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2006[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2006[i]=(x2006[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2007[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2007[i]=0
if (x2007[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2007[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2007[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2007[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2007[i]=(x2007[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2008[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2008[i]=0
if (x2008[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2008[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2008[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2008[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2008[i]=(x2008[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2009[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2009[i]=0
if (x2009[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2009[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2009[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2009[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2009[i]=(x2009[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2010[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2010[i]=0
if (x2010[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2010[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2010[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2010[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2010[i]=(x2010[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2011[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2011[i]=0
if (x2011[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2011[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2011[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2011[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2011[i]=(x2011[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2012[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2012[i]=0
if (x2012[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2012[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2012[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2012[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2012[i]=(x2012[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2013[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2013[i]=0
if (x2013[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2013[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2013[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2013[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2013[i]=(x2013[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2014[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2014[i]=0
if (x2014[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2014[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2014[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2014[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2014[i]=(x2014[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2015[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2015[i]=0
if (x2015[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2015[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2015[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2015[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2015[i]=(x2015[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2016[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2016[i]=0
if (x2016[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2016[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2016[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2016[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2016[i]=(x2016[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2017[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2017[i]=0
if (x2017[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2017[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2017[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2017[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2017[i]=(x2017[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2018[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2018[i]=0
if (x2018[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2018[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2018[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2018[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2018[i]=(x2018[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2019[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2019[i]=0
if (x2019[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2019[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2019[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2019[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2019[i]=(x2019[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
if (x2020[i]<(1990+yrsuntilPCT[i]*exp(b0_3))) x2020[i]=0
if (x2020[i]>(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2))) x2020[i]=0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i])))))
if (x2020[i]<(1990+yrsuntilPCT[i]*exp(b0_3)+yrsuntilmaxcov[i]*exp(b0_2)) & x2020[i]>(1990+yrsuntilPCT[i]*exp(b0_3))) x2020[i]=(x2020[i]-(1990+yrsuntilPCT[i]*exp(b0_3)))*((0.95/(1+(exp(-(b0_1+b1*logitP[i]+invnormal(uniform1[i]))))))/(yrsuntilmaxcov[i]*exp(b0_2)))
}
f=(mean(x1990))^2+(mean(x1991))^2+(mean(x1992))^2+(mean(x1993)-0.1)^2+(mean(x1994)-0.1)^2+(mean(x1995)-0.15)^2+(mean(x1996)-0.2)^2+(mean(x1997)-0.4)^2+(mean(x1998)-0.6)^2+(mean(x1999)-0.7)^2 +(mean(x2000)-0.7)^2+(mean(x2001)-0.75)^2+(mean(x2002)-0.8)^2+(mean(x2003)-0.8)^2+(mean(x2004)-0.8)^2+(mean(x2005)-0.8)^2+(mean(x2006)-0.8)^2+(mean(x2007)-0.8)^2+(mean(x2008)-0.8)^2+(mean(x2009)-0.8)^2+(mean(x2010)-0.8)^2+(mean(x2011)-0.8)^2+(mean(x2012)-0.8)^2+(mean(x2013)-0.8)^2+(mean(x2014)-0.8)^2+(mean(x2015)-0.8)^2+(mean(x2016)-0.8)^2+(mean(x2017)-0.8)^2+(mean(x2018)-0.8)^2+(mean(x2019)-0.8)^2+(mean(x2020)-0.8)^2
f
end
