# likelihood for Fister's twin criminals data (see notes on proportions)

freq=6652325*p  

# where p is the vector of null (central) hypergeometric frequences..

## psi > 1 pushes noncentral distribution to the right, psi < 1 pushes it to left

# could also use p itself without the scaling and the 6652325 cancels out anyway

L = function(psi) freq[11] * psi^10 / (sum( freq * psi^(0:12)))

pmf = function(psi) freq * psi^(0:12) / (sum( freq * psi^(0:12))) 



# transcript of in-class session...

> p=dhyper(0:12, 12,18, 13 )
> p
 [1] 7.154318e-05 1.860123e-03 1.753830e-02 8.038387e-02 2.009597e-01 2.893819e-01
 [7] 2.455362e-01 1.227681e-01 3.541387e-02 5.621250e-03 4.497000e-04 1.533068e-05
[13] 1.503008e-07
> sum(p)
[1] 1
> plot(0:12,p)
> dhyper(0:12, 13,17, 12 )
 [1] 7.154318e-05 1.860123e-03 1.753830e-02 8.038387e-02 2.009597e-01 2.893819e-01
 [7] 2.455362e-01 1.227681e-01 3.541387e-02 5.621250e-03 4.497000e-04 1.533068e-05
[13] 1.503008e-07
> sum(p[11:13])
[1] 0.0004651809
> 2150/6652325
[1] 0.0003231953
> p
 [1] 7.154318e-05 1.860123e-03 1.753830e-02 8.038387e-02 2.009597e-01 2.893819e-01
 [7] 2.455362e-01 1.227681e-01 3.541387e-02 5.621250e-03 4.497000e-04 1.533068e-05
[13] 1.503008e-07
> p[11:13]
[1] 4.497000e-04 1.533068e-05 1.503008e-07
> 1/2150
[1] 0.0004651163
> p
 [1] 7.154318e-05 1.860123e-03 1.753830e-02 8.038387e-02 2.009597e-01 2.893819e-01
 [7] 2.455362e-01 1.227681e-01 3.541387e-02 5.621250e-03 4.497000e-04 1.533068e-05
[13] 1.503008e-07
> 6652325*p
 [1] 4.759285e+02 1.237414e+04 1.166705e+05 5.347396e+05 1.336849e+06 1.925063e+06
 [7] 1.633386e+06 8.166932e+05 2.355846e+05 3.739438e+04 2.991550e+03 1.019847e+02
[13] 9.998497e-01

# likelihhod (noncentral hypergeometric) based on an observed a of 10


> freq=6652325*p


> L = function(psi) freq[11] * psi^10 / (sum( freq * psi^(0:12)))


> L(1)
[1] 0.0004497

> L(2)
[1] 0.008070642

> L(30)
[1] 0.3522977

> L(22)
[1] 0.3744265

> L(21)
[1] 0.3745935

> L(21.3)
[1] 0.3746362

# just about at top of the L fn. !!



# density (non-central Hypergeometric)

> pmf = function(psi) freq * psi^(0:12) / (sum( freq * psi^(0:12)))
> pmf(1)
 [1] 7.154318e-05 1.860123e-03 1.753830e-02 8.038387e-02 2.009597e-01 2.893819e-01
 [7] 2.455362e-01 1.227681e-01 3.541387e-02 5.621250e-03 4.497000e-04 1.533068e-05
[13] 1.503008e-07
> pmf(21.5)
 [1] 2.823910e-15 1.578566e-12 3.199978e-10 3.153312e-08 1.694905e-06 5.247427e-05
 [7] 9.572578e-04 1.029052e-02 6.382102e-02 2.178019e-01 3.746193e-01 2.745789e-01
[13] 5.787692e-02
> plot(0:12,pmf(21.5))
> plot(0:12,pmf(1.5))

# balancing equations to get value of psi so E(a | psi) = observed a

> E.a = sum(pmf(21.5)*(0:12)) ; E.a
[1] 10.00992

> E.a = sum(pmf(21.3)*(0:12)) ; E.a
[1] 9.999727
>