sshist_2d.pro

; Author: Shigenobu Hirose at JAMSTEC
; based on original paper
; Shimazaki and Shinomoto, Neural Computation 19, 1503-1527, 2007
;
function sshist_2d, data1, data2, x1=x1, x2=x2, cost=cost, nbin=nbin

  COMPILE_OPT idl2

  nbin_min = 2
  nbin_max = 200
  
  ntrial = nbin_max - nbin_min + 1
  
  nbin   = INDGEN(ntrial) + nbin_min
  
  delta1 = FLTARR(ntrial)
  delta2 = FLTARR(ntrial)
  cost   = FLTARR(ntrial)
  
  for n = 0, ntrial-1  do begin
     delta1[n] = (MAX(data1) - MIN(data1)) / (nbin[n] - 1) * (1. - (MACHAR()).EPS)
     delta2[n] = (MAX(data2) - MIN(data2)) / (nbin[n] - 1) * (1. - (MACHAR()).EPS)
     
     k = HIST_2D(data1, data2, bin1=delta1[n], bin2=delta2[n], min1=MIN(data1), min2=MIN(data2))
     
     kmean = MEAN(k)
     kvari = MEAN((k - kmean)^2)
     cost[n] = (2. * kmean - kvari) / (delta1[n] * delta2[n])^2
  endfor
  
  n = (WHERE(cost eq MIN(cost)))[0]
  k = HIST_2D(data1, data2, bin1=delta1[n], bin2=delta2[n], min1=MIN(data1), min2=MIN(data2))
  
  x1 = FINDGEN(nbin[n]) * delta1[n] + MIN(data1)
  x2 = FINDGEN(nbin[n]) * delta2[n] + MIN(data2)

  return, k
  
end