/**********************************************************************/


data sample (keep=id time conc weight)            
     parms (keep=time conc cmax tmax              
            rename=(time=lasttime conc=lastconc));
   infile cards eof=out;                          
   retain cmax tmax time conc;                    
   input id time conc weight;                     
   output sample;                                 
   if conc>cmax then do;                          
      cmax=conc;                                  
      tmax=time;                                  
   end;                                           
return;                                           

out:                        
   output parms;            
   cards;                   
1 0 0 0                     
1 0.5 0 0                   
1 1 1.15 0                  
1 1.5 1.98 0                
1 2 2.82 0                  
1 8 0.54 1                  
run;                        
                            
data sample;                
   set sample;              
   if _n_=1 then set parms; 
run;                        
                            
proc print;                 
run; 


/**********************************************************************/


proc model data=compart;                              
   endogenous conc;                                   
   exogenous time;                                    
   parms a1 95.0 a2 120.0 theta1 1.5 theta2 0.12;     
   conc=(a1*exp(-theta1*time))+(a2*exp(-theta2*time));
   fit conc start=(a1 90 95 100 a2 110 120 130        
       theta1 1 1.5 2 theta2 .10 .12 .14);            
run;                                                  
quit;


/**********************************************************************/
                                                                                                                          

    /* Implement parameter BOUNDS using penalty functions */
data nl1ka;                                                 
   input patient time conc @@;                              
   dose=51500;                                              
   wt=1;                                                    
   cards;                                                   
2 .5 6.68 2 1 7.5 2 2 6.64 2 4 6.16 2 8 5.23 2 12 4.24      
2 18 3.52 2 24 2.59 2 36 1.26 2 48 .77 2 60 .45 2 72 .21    
;
                                                            
proc model data=nl1ka ;                                     
   endogenous conc;                                         
   exogenous time;                                          
                                                            
      /* Define regression function                       */
   conc=a*exp(-k*time)-a*exp(-ka*time);                     
                                                            
      /* Deny infeasible set so we don't step out         */
   if a le 0 or ka le 0 or ka ge 30 or k le 0 then stop;    
                                                            
      /* Add penalty functions to deflect                 */
      /* search from boundary                             */
   penalty=   p1 / (a)                            /* a>0  */
            + p2 / (ka)                           /* ka>0 */
            + p3 / (30-ka)                       /* ka<30 */
            + p4 / (k) ;                         /* k>0   */
   resid.conc=sqrt(resid.conc ** 2 + scale * penalty );     
                                                            
      /* Set initial parameter guess                      */
   parms a 5 ka 25 k .01;                                   
                                                            
      /* Set initial penalty weights                      */
   control p1 1 p2 1 p3 1 p4 1 scale 1e-5;                  
                                                            
      /* Get initial estimates and print results at       */
      /* each iteration                                   */
   fit conc / method=gauss iter=100 itprint;                
   run;                                                     


/**********************************************************************/