(*
   Burton Rosenberg
   6:43 PM 10/5/00
*)

(* include packages *)
<< Statistics`ContinuousDistributions`

(* define functions *)

rn[ m_, sd_ ] := Random[NormalDistribution[m, sd]]

myClip[x_, s_, d_] := If[ (x > s + d), 
                          (s + d), 
                          If[  (x < s - d), 
                               s - d, 
                               x]
                        ]

rw2[n_, c_, m_, sd_] := Module[ 
                        {    s1 = 0,
                             sl1 = {0},
                             s2 = 0,
                             sl2 = {0}
                        },
                        Do[ s1 = s1 + rn[ m, sd ] ;
                            s2 = myClip[ s1, s2, c] ;
                            sl1 = Append[sl1, s1] ;
                            sl2 = Append[sl2, s2] ,
                            {n}] ;
                        { sl1, sl2 }
                      ]

lastXInRange[ r_, rstar_ ] := Module[ { n = 0,
                                        x = 0
                                      },
                                    Do[  If[ ( r[[i]]==rstar[[i]] ),
                                             n = i ; x = r[[i]]
                                         ]
                                         ,{ i, 1, Length[r]}
                                    ] ;
                                    { x, n }
                               ]

muEstimate[ r_, rstar_ ] := Module[ { x, n },
                                    {x,n} = lastXInRange[r, rstar] ;
                                    If[ (n>1),
                                         x/(n-1)
                                    ]
                            ]

uv[ r_, rstar_] := Module[ { u = {},
                             v = {},
                             lookForUncollared = True
                           },
                       Do[ 
                            If[ (lookForUncollared),
                                     (* looking for uncollared *)
                                     If[ ( r[[i]]==rstar[[i]] ),
                                           v = Append[ v, i ] ; 
                                           lookForUncollared = False
                                     ],
                                     (* looking for next collared *)
                                     If[ ( r[[i]]!=rstar[[i]] ),
                                           u = Append[ u, i ] ;
                                           lookForUncollared = True
                                     ]
                            ]
                            , {i, 1, Length[r]}
                       ] ;
                       { v, u }  
                   ]

sumOfWaits[ r_, rstar_ ] := Module[ { v, u, wTotal=0 },
                                    { v,u } = uv[r,rstar] ;
                                    Do[ wTotal += u[[i]] - v[[i]]
                                    , { i, 1, Length[u] }
                                    ] ;
                                    { wTotal, Length[u] }
                                ]



(*
  now solve: 
  { wTotal, m } = sumOfWaits[ r, rstar] ;
                  1 - ( m/wTotal ) = Integral[ { x, -collar, +collar },
                                          (1/(Sqrt[2 Pi] sigma)) 
                                             Exp[ -(x-muEstimate)^2/(2 sigma^2) ]
                                     ]
  for sigma
*)


myInt[ limit_, sigma_ ] := Integrate[
                             ( 1 / (sigma Sqrt[ 2 Pi ]) ) Exp[ - t^2 / ( 2 sigma^2 ) ], 
                             { t, 0, limit }
                           ]



ro[l_] := Module[ { ltemp = {} }, 
                  Do[ ltemp = Append[ltemp, {i, Exp[l[[i]]]}], 
                      {i, 1, Length[l]}]; 
                  ltemp
                ]

(* sample run *)
{ s, sStar} = rw2[100, .1, 0, .1]
{ v, u } = uv[ s, sStar ]

myShow[l1_,l2_] := Show[
                       Graphics[{RGBColor[1, 0, 0], 
                                Line[ ro[l1] ], 
                                RGBColor[0, 0, 1], 
                                Line[ ro[l2] ]}
                       ],
                       {Axes -> Automatic}]

myShow[s,sStar]






