/* 

  Secretary problem in Picat.

  From Siegrist "Probability, Mathematical Statistics, and Stochastic Processes"
  Section "12.9: The Secretary Problem"

  Cf
  - ppl_sultans_dowry.pi 

  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

import ppl_distributions, ppl_utils.
import util.
% import ordset.

main => go.


/*
  """
  With n candidates, let X denote the
  number (arrival order) of the best candidate, and let S denote the event
  of success for strategy k (we select the best candidate).
  Xn is uniformly distributed on 1..n.
  Next we will compute the conditional probability of success given the arrival order of
  the best candidate.
  For n in N+, and k 2..n
                    0 if j in 1.. k-1
  P(S(n,k) | Xn=j)  (k-1)/(j-1) j in k..n
  """
*/
secretary_p(N,K,J) = Res =>
  if J < K then
    Res = 0
  else
    Res = (K-1) / (J-1)
  end.

pn(N,K) = Res =>
 if K == 1 then
   Res = 1 / N
 else
   Res = ((K - 1) / N) * sum([ 1/(J-1) : J in K..N])
 end.

pnk(N,K) = Res =>
   Res = (K - 1) / sum([ (1/N) * N / (J-1) : J in K..N]).

/*

  N = 6

  secretary_p(N,K,J):
  [k = 2,j = 1,p = 0]
  [k = 2,j = 2,p = 1.0]
  [k = 2,j = 3,p = 0.5]
  [k = 2,j = 4,p = 0.333333]
  [k = 2,j = 5,p = 0.25]
  [k = 3,j = 1,p = 0]
  [k = 3,j = 2,p = 0]
  [k = 3,j = 3,p = 1.0]
  [k = 3,j = 4,p = 0.666667]
  [k = 3,j = 5,p = 0.5]
  [k = 4,j = 1,p = 0]
  [k = 4,j = 2,p = 0]
  [k = 4,j = 3,p = 0]
  [k = 4,j = 4,p = 1.0]
  [k = 4,j = 5,p = 0.75]

*/
go ?=>
  N = 6,
  println("secretary_p(N,K,J):"),
  foreach(K in 2..N-1)
    foreach(J in 1..N)
      println([k=K,j=J,p=secretary_p(N,K,J)])
    end
  end,
  nl.
go => true.

/*
  N = 10

  pn(N,K):
  [k = 1,pn = 0.1]
  [k = 2,pn = 0.282897]
  [k = 3,pn = 0.365794]
  [k = 4,pn = 0.39869]
  [k = 5,pn = 0.398254]
  [k = 6,pn = 0.372817]
  [k = 7,pn = 0.327381]
  [k = 8,pn = 0.265278]
  [k = 9,pn = 0.188889]
  [k = 10,pn = 0.1]

  pnk(N,K):
  [k = 2,pnk = 0.353486]
  [k = 3,pnk = 1.09351]
  [k = 4,pnk = 2.25739]
  [k = 5,pnk = 4.01754]
  [k = 6,pnk = 6.70569]
  [k = 7,pnk = 10.9964]
  [k = 8,pnk = 18.4712]
  [k = 9,pnk = 33.8824]
  [k = 10,pnk = 81.0]



*/
go2 ?=>
  N = 10,
  println("pn(N,K):"),
  foreach(K in 1..N)
    println([k=K,pn=pn(N,K)])
  end,
  nl,
  println("pnk(N,K):"),
  foreach(K in 2..N)
    println([k=K,pnk=pnk(N,K)])
  end,
  nl.
go2 => true.
  

/*
  N = 100

  PN:
  pn = [0.01,0.0517738,0.0835476,0.110321,0.133762,0.154702,0.173643,0.190916,0.206762,0.221357,0.234841,0.247325,0.2589,0.269642,0.279614,0.288872,0.297464,0.30543,0.312808,0.319631,0.325928,0.331724,0.337044,0.34191,0.346341,0.350355,0.353969,0.357199,0.360058,0.36256,0.364717,0.366541,0.368042,0.369231,0.370117,0.370709,0.371015,0.371043,0.370801,0.370295,0.369534,0.368522,0.367267,0.365773,0.364047,0.362093,0.359918,0.357524,0.354919,0.352104,0.349086,0.345868,0.342453,0.338847,0.335051,0.331071,0.326909,0.322568,0.318051,0.313363,0.308504,0.303479,0.29829,0.29294,0.287431,0.281766,0.275947,0.269977,0.263857,0.25759,0.251179,0.244624,0.237929,0.231094,0.224123,0.217016,0.209777,0.202405,0.194904,0.187275,0.179519,0.171638,0.163633,0.155507,0.14726,0.138894,0.13041,0.12181,0.113096,0.104267,0.0953262,0.0862743,0.0771125,0.067842,0.0584639,0.0489795,0.0393898,0.0296959,0.019899,0.01]
  [argmax = 38,ix = 0.371043]

  PNK:
  pnk = [0.193148,0.478769,0.815799,1.19616,1.61601,2.07322,2.56657,3.09535,3.65925,4.2582,4.89235,5.56199,6.26758,7.00966,7.78891,8.60609,9.46206,10.3578,11.2943,12.2727,13.2942,14.3601,15.4719,16.631,17.8391,19.0977,20.4088,21.7743,23.1962,24.6767,26.2181,27.8229,29.4937,31.2334,33.0448,34.9312,36.896,38.9428,41.0753,43.2978,45.6146,48.0305,50.5505,53.18,55.9248,58.7912,61.786,64.9163,68.19,71.6156,75.2021,78.9596,82.8988,87.0314,91.3702,95.929,100.723,105.769,111.085,116.692,122.611,128.868,135.488,142.504,149.947,157.856,166.274,175.246,184.828,195.08,206.071,217.881,230.599,244.33,259.197,275.34,292.927,312.154,333.254,356.509,382.259,410.919,443.003,479.153,520.182,567.134,621.376,684.731,759.684,849.714,959.846,1097.62,1274.87,1511.36,1842.61,2339.69,3168.45,4826.38,9801.0]

*/
go3 ?=>
  N = 100,
  println("PN:"),
  PN = [pn(N,K) :  K in 1..N],
  println(pn=PN),
  PNArgMax = argmax(PN).first,
  
  println([argmax=PNArgMax,ix=PN[PNArgMax]]),
  nl,
  println("PNK:"),
  PNK = [pnk(N,K) :  K in 2..N],
  println(pnk=PNK),
  nl.
go3 => true.