0
$\begingroup$

I define a Random Variable drawn from a custom probability density.

u[xi_, yi_, zi_, τ_, δt_, θ0_, k_] := 
 RandomVariate[ProbabilityDistribution[(Sqrt[δt/τ]*(Cosh[(u*δt)/τ] +
(zi*Cos[θ0] + Sin[θ0]*(xi*Cos[k*δt] + yi*Sin[k*δt]))*Sinh[(u*δt)/τ]))/
(E^(((1 + u^2)*δt)/(2*τ))*Sqrt[2*Pi]), {u, -Infinity, Infinity}]]

I want to use this random variable inside the delayed function xfun:

xfun[xi_, yi_, zi_, θ0_, k_, δt_, τ_] := (xi - xi*Cos[2*θ0] - 4*Cos[2*k*Pi*δt]*Sin[θ0] + 
  2*xi*Cos[4*k*Pi*δt]*Sin[θ0]^2 + 2*yi*Sin[4*k*Pi*δt]*Sin[θ0]^2 +
  4*E^((2*δt*u[xi, yi, zi, τ, δt, θ0, k])/τ)*Cos[2*k*Pi*δt]*Sin[θ0]*
  (1 + zi*Cos[θ0] + xi*Cos[2*k*Pi*δt]*Sin[θ0] + yi*Sin[2*k*Pi*δt]*Sin[θ0]) +
  2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0] + 2*E^((δt*u[xi, yi, zi, τ, δt, θ0, k])/τ)* 
  (3*xi + xi*Cos[2*θ0] - 2*(xi*Cos[4*k*Pi*δt] + yi*Sin[4*k*Pi*δt])*
  Sin[θ0]^2 - 2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0]))/(-4*(-1 + zi*Cos[θ0] + 
  xi*Cos[2*k*Pi*δt]*Sin[θ0] + yi*Sin[2*k*Pi*δt]*Sin[θ0]) +
  4*E^((2*δt*u[xi, yi, zi, τ, δt, θ0, k])/τ)*(1 + zi*Cos[θ0] + 
  xi*Cos[2*k*Pi*δt]*Sin[θ0] + yi*Sin[2*k*Pi*δt]*Sin[θ0]))

Now checking the special case (specified by the arguments of xfun) I find that I will sometimes get values greater than 1 for my function

xfun[1, 0, 0, Pi/2, 1, 0.01, 1]

This should be impossible since I can check, replacing the function u[,,,,] with the symbol u in the RHS of the function xfun, that the maximum value xfun should take for any u in this special case should be 1.

Maximize[(xi - xi*Cos[2*θ0] - 4*Cos[2*k*Pi*δt]*Sin[θ0] +
 2*xi*Cos[4*k*Pi*δt]*Sin[θ0]^2 + 
  2*yi*Sin[4*k*Pi*δt]*Sin[θ0]^2 + 4*E^((2*δt*u)/τ)*Cos[2*k*Pi*δt]*Sin[θ0]*
  (1 + zi*Cos[θ0] + xi*Cos[2*k*Pi*δt]*Sin[θ0] + yi*Sin[2*k*Pi*δt]*Sin[θ0]) + 
  2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0] + 2*E^((δt*u)/τ)*(3*xi + xi*Cos[2*θ0] - 
  2*(xi*Cos[4*k*Pi*δt] + yi*Sin[4*k*Pi*δt])*Sin[θ0]^2 -  2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0]))/
  (-4*(-1 + zi*Cos[θ0] + xi*Cos[2*k*Pi*δt]*Sin[θ0] +  yi*Sin[2*k*Pi*δt]*Sin[θ0]) +    
  4*E^((2*δt*u)/τ)*(1 + zi*Cos[θ0] + xi*Cos[2*k*Pi*δt]*Sin[θ0] + 
  yi*Sin[2*k*Pi*δt]*Sin[θ0])) /. {θ0 -> Pi/2, k -> 1, δt -> 0.01, τ -> 1, 
   xi -> 1, yi -> 0, zi -> 0}, u]

Does anybody know what the source of the problem could be? Is there something about the behavior of the RandomVariate function that I am not understanding? Thanks!

$\endgroup$

2 Answers 2

2
$\begingroup$

In the definition of xfun you called u three times. Each time generated a different random variable. Use Module to structure xfun to call u only once and use that value all three times.

Clear["Global`*"]

u[xi_, yi_, zi_, τ_, δt_, θ0_, k_] := 
 RandomVariate[
  ProbabilityDistribution[(Sqrt[δt/τ]*(Cosh[(u*δt)/τ] \
+ (zi*Cos[θ0] + 
           Sin[θ0]*(xi*Cos[k*δt] + yi*Sin[k*δt]))*
         Sinh[(u*δt)/τ]))/(E^(((1 + u^2)*δt)/(2*τ))*
      Sqrt[2*Pi]), {u, -Infinity, Infinity}]]

Modified definition of xfun

xfun[xi_, yi_, zi_, θ0_, k_, δt_, τ_] := 
 Module[{urv = u[xi, yi, zi, τ, δt, θ0, k]},
  (xi - xi*Cos[2*θ0] - 4*Cos[2*k*Pi*δt]*Sin[θ0] + 
     2*xi*Cos[4*k*Pi*δt]*Sin[θ0]^2 + 
     2*yi*Sin[4*k*Pi*δt]*Sin[θ0]^2 + 
     4*E^((2*δt*urv)/τ)*Cos[2*k*Pi*δt]*
      Sin[θ0]*(1 + zi*Cos[θ0] + 
        xi*Cos[2*k*Pi*δt]*Sin[θ0] + 
        yi*Sin[2*k*Pi*δt]*Sin[θ0]) + 
     2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0] + 
     2*E^((δt*urv)/τ)*(3*xi + xi*Cos[2*θ0] - 
        2*(xi*Cos[4*k*Pi*δt] + yi*Sin[4*k*Pi*δt])*
         Sin[θ0]^2 - 
        2*zi*Cos[2*k*Pi*δt]*Sin[2*θ0]))/(-4*(-1 + 
        zi*Cos[θ0] + xi*Cos[2*k*Pi*δt]*Sin[θ0] + 
        yi*Sin[2*k*Pi*δt]*Sin[θ0]) + 
     4*E^((2*δt*urv)/τ)*(1 + zi*Cos[θ0] + 
        xi*Cos[2*k*Pi*δt]*Sin[θ0] + 
        yi*Sin[2*k*Pi*δt]*Sin[θ0]))]

Testing,

SeedRandom[1234];

Max@Table[xfun[1, 0, 0, Pi/2, 1, 0.01, 1], {50}]

(* 1. *)
$\endgroup$
0
$\begingroup$

This is not an answer but just an observation. Occasionally there are errors such as the following (with just a slight modification of your function u):

u2[xi_, yi_, zi_, τ_, δt_, θ0_, k_, n_] := 
 RandomVariate[ProbabilityDistribution[(Sqrt[δt/τ]*(Cosh[(u*δt)/τ] + 
 (zi*Cos[θ0] + Sin[θ0]*(xi*Cos[k*δt] + yi*Sin[k*δt]))*
         Sinh[(u*δt)/τ]))/(E^(((1 + u^2)*δt)/(2*τ))*Sqrt[2*Pi]), {u, -Infinity, Infinity}], n]

x = u2[1, 0, 0, 1, 0.01, π/2, 1, 1000000];

Error messages

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.