#!/usr/bin/perl
################################################################
#
#  EM for simple linear interpolation between two fixed models
# 
#
#  To run, just execute
#
#     ./interpolation_em.prl
#
#  You can edit the values under "Models and Initialization"
#  to play with different models for p1 and p2, different
#  initial values for the lambdas, or different 
#  proportions in observed data.
#
#  Under program parameters set $verbose = 1 for verbose
#  output showing intermediate values in the updates,
#  and modify $num_iterations as you see fit. (This version
#  iterates a fixed number of times rather than seeing if
#  the lambdas have converged.)
#
#
#  Author: Philip Resnik
#  Date:   7 October 2005
#
################################################################


###
### Program parameters
###

# Set to nonzero for verbose output
my $verbose = 1;

# Number of iterations
my $num_iterations = 100;



###
### Models and initialization
###

# Initial lambdas (weights) for states 1 and 2
# Note: $lambda[0] has no meaning since there's no state 0.
my @lambda = ("ignore", 0.5, 0.5);

# Observables
my @Y = ("j", "c", "d");

# Observed proportion
my %ptilde = (
  "j" => .1,
  "c" => .6,
  "d" => .3	
);


# Fixed distribution for state (model) 1
my %p1 = ( 
  "j" => .1,
  "c" => .6,
  "d" => .3
);

# Fixed distribution for state (model) 2
my %p2 = ( 
  "j" => .2,
  "c" => .2,
  "d" => .6
);


###
### EM
###

foreach (my $iteration = 0; $iteration < $num_iterations; ++$iteration)
{
  print "\nIteration = $iteration\n" if ($verbose > 0);

  $C1 = 0;
  foreach my $y (@Y) 
  {
      my $p1_contribution += ($lambda[1] * $p1{$y});
      my $p2_contribution += ($lambda[2] * $p2{$y});
      my $denominator = $p1_contribution + $p2_contribution;

      print "C1: ptilde($y) = $ptilde{$y}, lambda1 = $lambda[1], p1($y) = $p1{$y}, p2($y) = $p2{$y}, denominator = $denominator\n"
	if ($verbose > 0);

      my $to_add = $ptilde{$y} * $p1_contribution / $denominator;

      print "C1: for y=$y, adding $ptilde{$y} * $p1_contribution / ($p1_contribution + $p2_contribution) = $to_add\n"
	if ($verbose > 0);

      $C1 += $to_add;
  }

  print "\n" if ($verbose > 0);

  $C2 = 0;
  foreach my $y (@Y) 
  {
      my $p1_contribution += ($lambda[1] * $p1{$y});
      my $p2_contribution += ($lambda[2] * $p2{$y});
      my $denominator = $p1_contribution + $p2_contribution;
      my $to_add = $ptilde{$y} * $p2_contribution / $denominator;

      print "C2: ptilde($y) = $ptilde{$y}, lambda2 = $lambda[2], p1($y) = $p1{$y}, p2($y) = $p2{$y}\n"
	if ($verbose > 0);

      print "C2: for y=$y, adding $ptilde{$y} * $p2_contribution / ($p1_contribution + $p2_contribution) = $to_add\n"
	if ($verbose > 0);

      $C2 += $to_add;
  }


  print "At end of E step, C1=$C1, C2=$C2\n" if ($verbose > 0);

  $newlambda1 = $C1 / ($C1 + $C2);
  $newlambda2 = $C2 / ($C1 + $C2);

  printf "New lambdas at end of iteration %3d:  %8.7f  %8.7f\n", $iteration, $newlambda1, $newlambda2;

  @lambda = ("ignore", $newlambda1, $newlambda2);
}