/* ----------------------------------------------------------------------------

Copyright (C) 2005, 2006.

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-----------------------------------------------------------------------------*/

#include "libscl.h"
#include "libgsm.h"
#include "gsmusr.h"
using namespace std;  
using namespace scl;
using namespace libgsm;
using namespace gsm;

den_val 
one_dim_tinker_toy_sci_mod::prior (const realmat& parm, const realmat& func)
{
  return den_val(true,0.0);
}

bool 
one_dim_tinker_toy_sci_mod::gen_sim(realmat& sim, realmat& func)
{
  if (theta[1] <= 0.0) return false;

  INT_32BIT seed = 770116;

  sim.resize(1,len_sim);
  func.resize(3,1);
  REAL mean = 0;
  REAL var = 0;
  
  for (INTEGER i=1; i<=len_sim; ++i) {
    sim[i] = theta[1] + theta[1]*scl::unsk(seed);
    mean += sim[i];
  }

  mean /= REAL(len_sim);

  for (INTEGER i=1; i<=len_sim; ++i) {
    var += pow(sim[i]-mean,2);
  }

  var /= REAL(len_sim-1);

  func[1] = mean;
  func[2] = var;
  func[3] = mean/sqrt(var);

  return true;
}
    
den_val 
one_dim_tinker_toy_stat_mod::prior
  (const realmat& parm, const realmat& func)
{
  return den_val(true,0.0);
}
  
bool 
one_dim_tinker_toy_stat_mod::gen_sim(realmat& sim, realmat& func)
{
  INT_32BIT seed = 770116;

  sim.resize(1,len_sim);
  func.resize(3,1);
  REAL mean = 0;
  REAL var = 0;

  REAL mu = eta[1];
  REAL sd = eta[2];

  for (INTEGER i=1; i<=len_sim; ++i) {
    sim[i] = mu + sd*scl::unsk(seed);
    mean += sim[i];
  }

  mean /= REAL(len_sim);

  for (INTEGER i=1; i<=len_sim; ++i) {
    var += pow(sim[i]-mean,2);
  }

  var /= REAL(len_sim-1);

  func[1] = mean;
  func[2] = var;
  func[3] = mean/sqrt(var);

  return true;
}

REAL
one_dim_tinker_toy_stat_mod::penalty(realmat& dpenalty)
{
  dpenalty.resize(1,len_parm());
  smomax smooth;
  dpenalty[1] = 0.0;
  dpenalty[2] = smooth.max1(-eta[2]); 
  return smooth.max0(-eta[2]);
}

den_val 
one_dim_tinker_toy_stat_mod::loglikelihood()
{
  REAL mu = eta[1];
  REAL sd = eta[2];

  const REAL half_log_two_pi = 0.91893853320467266954;

  INTEGER n = data_ptr->size();

  REAL sumsq = 0.0;
  for (INTEGER i=1; i<=n; ++i) { 
    REAL x = (*data_ptr)[i]-mu;
    sumsq += pow(x,2); 
  }

  REAL log_top = (-0.5)*sumsq/pow(sd,2);

  REAL log_bot = REAL(n)*(half_log_two_pi + log(sd));

  REAL log_den = log_top - log_bot;

  return den_val(true, log_den);
}

den_val 
one_dim_tinker_toy_stat_mod::loglikelihood(realmat& dlogl)
{
  dlogl.resize(1,2);

  REAL mu = eta[1];
  REAL sd = eta[2];

  const REAL half_log_two_pi = 0.91893853320467266954;

  INTEGER n = data_ptr->size();

  REAL sum = 0.0;
  REAL sumsq = 0.0;
  for (INTEGER i=1; i<=n; ++i) { 
    REAL x = (*data_ptr)[i]-mu;
    sumsq += pow(x,2); 
    sum += x;
  }

  REAL log_top = (-0.5)*sumsq/pow(sd,2);

  REAL log_bot = REAL(n)*(half_log_two_pi + log(sd));

  REAL log_den = log_top - log_bot;

  dlogl[1] = sum/pow(sd,2);

  dlogl[2] = sumsq/pow(sd,3) - REAL(n)/sd;

  return den_val(true, log_den);
}

INTEGER 
one_dim_tinker_toy_stat_mod::mle(realmat& stat_parm_mle, realmat& stat_parm_V)
{
  stat_parm_mle.resize(2,1);
  stat_parm_V.resize(2,2);

  INTEGER n = data_ptr->size();

  REAL sum = 0.0;
  for (INTEGER i=1; i<=n; ++i) { 
    sum += (*data_ptr)[i];
  }
  REAL mu = sum/REAL(n);

  REAL sumsq = 0.0;
  for (INTEGER i=1; i<=n; ++i) { 
    REAL x = (*data_ptr)[i]-mu;
    sumsq += pow(x,2); 
  }
  REAL var = sumsq/REAL(n);
  REAL sd = sqrt(var);

  stat_parm_mle[1] = mu;
  stat_parm_mle[2] = sd;

  // Using distribution theory

  stat_parm_V(1,1) = var/REAL(n);
  stat_parm_V(2,1) = 0.0;
  stat_parm_V(1,2) = 0.0;
  stat_parm_V(2,2) = 0.5*var/REAL(n);

  return 1;

  /* 
  // An alternative approach that relies on stat_mod_parm being set to the 
  // previous mle by sci_mcmc::draw prior to calling sci_mcmc::make_sci_val 

  realmat parm;
  get_parm(parm);
  parm -= stat_parm_mle;
  stat_parm_V = parm*T(parm);

  return 1;
  */ 
}

